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 **mathematical model**, generally given by partial differential equations, by using **discretization** schemes.

## 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 calculation
* $F_{B,\sigma}$: flux over sub control volume face
* $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_h: &\:\text{global discrete solution, interpolation using \textbf{basis functions}} \\
    \mathbf{n}: &\:\text{unit outer normal vector} \\
    \sigma: &\:\text{sub control volume face} \\
\end{aligned}$


## Local residual
The local residual of the diffusion model:
```cpp
namespace Dumux {
template<class TypeTag>
class DiffusionModelLocalResidual
: public GetPropType<TypeTag, Properties::BaseLocalResidual>
{
```

## Local residual
```cpp
namespace Dumux {
template<class TypeTag>
class DiffusionModelLocalResidual
: public GetPropType<TypeTag, Properties::BaseLocalResidual>
{
```
Our class inherits from the class `BaseLocalResidual`, which has to be 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_h: &\:\text{discrete solution, interpolated using basis functions} \\
    \sigma: &\:\text{sub control volume face} \\
\end{aligned}$

## Flux term
Implementation takes place in the dedicated function `computeFlux`:
```cpp
 NumEqVector computeFlux(const Problem& problem,
                         const Element& element,
                         const FVElementGeometry& fvGeometry,
                         const ElementVolumeVariables& elemVolVars,
                         const SubControlVolumeFace& scvf,
                         const ElementFluxVariablesCache& elemFluxVarsCache) const
 {
    ...
```

## Flux term
Implementation takes place in the dedicated function `computeFlux`:
```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
Implementation takes place in the dedicated function `computeFlux`:
```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
By specializing properties for the type tag `DiffusionModel`, every other class that knows about the type tag (this will be for example the assembler or the problem), can extract the type information that we specify here.

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