model.md

---
title: Models in DuMuX
---

# Implementing a new DuMu<sup>X</sup> model

## What is a DuMu<sup>X</sup> model
A DuMu<sup>X</sup> model is an implementation of a discretized **mathematical model**, generally given by partial differential equations.

## What is a DuMu<sup>X</sup> model
Mathematical model (PDE):
$\begin{equation*}
\small
\frac{\partial S(u)}{\partial t} + \nabla \cdot \mathbf{F}(u) = q, \quad \forall (t,\mathbf{x}) \in (0,T] \times \Omega
\end{equation*}$

<img src=img/model_domain_discretization.png width="100%">

Discrete model, e.g. using finite volumes:
$\begin{equation*}
\small
|B| \frac{S_h(\mathbf{u}^{n+1}_h) - S_h(\mathbf{u}^{n}_h)}{\Delta t}  + \sum_{\sigma \in \Sigma_B} F_{B,\sigma}(\mathbf{u}^{n+1}_h) = \int_{B}  q^{n+1} \, dx, \quad \forall t_{n+1}\leq T, \; \forall B
\end{equation*}$

## What is a DuMu<sup>X</sup> model
Discrete model, e.g. using finite volumes:
$\begin{equation*}
\small
|B| \frac{S_h(\mathbf{u}^{n+1}_h) - S_h(\mathbf{u}^{n}_h)}{\Delta t}  + \sum_{\sigma \in \Sigma_B} F_{B,\sigma}(\mathbf{u}^{n+1}_h) = \int_{B}  q \, dx, \quad \forall t_{n+1}\leq T, \; \forall B
\end{equation*}$

* $S_h$: storage term
* $F_{B,\sigma}$: flux term over sub control volume face (scvf)
* $q$ source term

How to implement these terms in DuMu<sup>X</sup>?

**local residual**

# Local Residual

## Local Residual

Implements these terms within

* `computeStorage(...)`
* `computeFlux(...)`
* `computeSource(...)`


## Example: Diffusion equation
Mathematical model (PDE):

$\begin{equation}
\frac{\partial c}{\partial t} - \nabla \cdot (D \nabla c) = 0 \:\: \mathrm{in}\; \Omega \times (0,T] \\
\end{equation}$

with

$\begin{aligned}
    c: &\:\text{concentration} \\
    D: &\:\text{constant diffusion coefficient} \\
    \Omega: &\:\text{spatial domain} \\
    T: &\:\text{end time}
\end{aligned}$

## Example: Diffusion equation
Discrete model using the Box discretization:

$\begin{equation}
\vert B \vert \frac{c_B^{n+1}-c_B^{n}}{\Delta t} - \sum_{\sigma \in \Sigma_{B}} \left[ D \nabla c_h^{n+1} \cdot \boldsymbol{n}_{B,\sigma} \vert \sigma \vert \right] = 0,
\end{equation}$

with

$\begin{aligned}
    c_B^n: &\:\text{concentration at time $t_n$ and control volume $B$} \\
    c^n_h: &\:\text{global discrete solution at time $t_n$, interpolated using \textbf{basis functions}} \\
    \mathbf{n}: &\:\text{unit outer normal vector} \\
    \sigma: &\:\text{sub control volume face (scvf)} \\
\end{aligned}$

## Example: Diffusion equation
Discrete model using the Box discretization:

$\begin{equation}
\vert B \vert \frac{c_B^{n+1}-c_B^{n}}{\Delta t} - \sum_{\sigma \in \Sigma_{B}} \left[ D \nabla c_h^{n+1} \cdot \boldsymbol{n}_{B,\sigma} \vert \sigma \vert \right] = 0,
\end{equation}$

<img src=img/box_scv_scvf.png width="90%">

## Local residual
The local residual of the diffusion model:
```cpp
template<class TypeTag>
class DiffusionModelLocalResidual
: public GetPropType<TypeTag, Properties::BaseLocalResidual>
{
    ...
}
```
Inherits from the `BaseLocalResidual`, which is chosen depending on the discretization scheme, here *Box scheme*.

## Storage term
```cpp
NumEqVector computeStorage(const Problem& problem,
						   const SubControlVolume& scv,
						   const VolumeVariables& volVars) const
 {
    NumEqVector storage;
    storage[Indices::massBalanceEqIdx]
        = volVars.priVar(Indices::concentrationIdx);
    return storage;
 }
```

## Flux term

$\begin{equation}
F_{B,\sigma} = -D \nabla c_h^{n+1} \cdot \boldsymbol{n}_{B,\sigma} \vert \sigma \vert
\end{equation}$

with

$\begin{aligned}
    c^n_h: &\:\text{global discrete solution at time $t_n$, interpolated using \textbf{basis functions}} \\
    \mathbf{n}: &\:\text{unit outer normal vector} \\
    \sigma: &\:\text{sub control volume face (scvf)} \\
\end{aligned}$

## Flux term
```cpp
NumEqVector computeFlux(const Problem& problem,
                        const Element& element,
                        const FVElementGeometry& fvGeometry,
                        const ElementVolumeVariables& elemVolVars,
                        const SubControlVolumeFace& scvf,
                        const ElementFluxVariablesCache& elemFluxVarsCache) const
{
    ...
}
```

## Flux term
```cpp
NumEqVector computeFlux(...) const
{
    // Compute ∇c
    const auto& fluxVarCache = elemFluxVarsCache[scvf];
    Dune::FieldVector<Scalar, dimWorld> gradConcentration(0.0);
    for (const auto& scv : scvs(fvGeometry))
    {
        const auto& volVars = elemVolVars[scv];
        gradConcentration.axpy(
            volVars.priVar(Indices::concentrationIdx),
            fluxVarCache.gradN(scv.indexInElement())
        );
    }
    ...
}
```

## Flux term
```cpp
NumEqVector computeFlux(...) const
{
    ...
    NumEqVector flux;

    // Compute the flux
    flux[Indices::massBalanceEqIdx] = -1.0*vtmv(
        scvf.unitOuterNormal(), problem.diffusionCoefficient(), gradConcentration
    )*scvf.area();

    return flux;
}
```

## Local Residual

A **local residual** implements the discretized mathematical model.

For its implementation different model-specific properties have to be set

# Model properties

## Model type tag

The property tag is an empty struct with the respective name, e.g. `DiffusionModel`

```cpp
namespace Dumux::Properties::TTag {
//! The diffusion model tag that we can specialize properties for
struct DiffusionModel {};
} // end namespace Dumux::Properties::TTag
```
## Model properties
We can set nodel properties for the `DiffusionModel` type tag

All properties are defined within the namespace `Dumux::Properties`
```cpp
namespace Dumux::Properties {
    //define all properties
} // end namespace Dumux::Properties
```

## Defining model properties
The type of the local residual is the class `DiffusionModelLocalResidual` defined from earlier
```cpp
template<class TypeTag>
struct LocalResidual<TypeTag, TTag::DiffusionModel>
{ using type = DiffusionModelLocalResidual<TypeTag>; };
```

## Defining model properties
The model traits specify information about the model:
```cpp
template<class TypeTag>
struct ModelTraits<TypeTag, TTag::DiffusionModel>
{
    struct type
    {
        struct Indices
        {
		    static constexpr int concentrationIdx = 0;
		    static constexpr int massBalanceEqIdx = 0;
        };
        static constexpr int numEq() { return 1; }
    };
};
```
## Defining model properties
Further model specific properties can be set accordingly by using the model property tag,
i.e. `TTag::DiffusionModel`

# Exercise: Model

## Exercise: Model
Implementation of a **nonlinear diffusion model**

<figure>
    <center>
        <img src="../exercises/extradoc/exercisemodel_mri_denoise.gif" alt="denoising"/>
        <figcaption> Denosing of MRI image using nonlinear diffusion model.</figcaption>
    </center>
</figure>

## Tasks

- Implement local residual
- Set model properties
- Use model in test case
- Customize volume variables