model.md

---
title: Implementing a Model in DuMu^x^
---

# What is a DuMu^x^ model

## What is a DuMu^x^ model
A DuMu^x^ model is an implementation of a discretized **mathematical model**, generally given by partial differential equations.

## What is a DuMu^x^ model
Mathematical model (PDE):
$$
\begin{aligned}
\frac{\partial S(u)}{\partial t} + \nabla \cdot \mathbf{F}(u) = q(u), \quad \forall (t,\mathbf{x}) \in (0,T] \times \Omega
\end{aligned}
$$

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

## What is a DuMu^x^ model
Discrete model using finite volumes:
$$
\begin{aligned}
\small
|B| \frac{S_B(\mathbf{u}^{n+1}_h) - S_B(\mathbf{u}^{n}_h)}{\Delta t}  + \sum_{\sigma \in \Sigma_B} F_{B,\sigma}(\mathbf{u}^{n+1}_h) = |B|  q_B(\mathbf{u}^{n+1}_h), \quad \forall t_{n+1}\leq T, \; \forall B
\end{aligned}
$$

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

Where to implement these terms in DuMu^x^?

`LocalResidual`

# `LocalResidual`

## `LocalResidual`

Implements terms of a PDE in the functions

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


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

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

with

- $c:$ concentration
- $D:$ constant diffusion coefficient
- $\Omega:$ spatial domain
- $T:$ end time

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

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

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

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

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

with

- $c_B^n:$ concentration at time $t_n$ and control volume $B$
- $c^n_h:$ global discrete solution at time $t_n$, interpolated using __basis functions__
- $\mathbf{n}:$ unit outer normal vector
- $\sigma:$ sub control volume face (scvf)

## `LocalResidual`
The local residual of the diffusion model:
```cpp
template<class TypeTag>
class DiffusionModelLocalResidual
: public DiscretizationDefaultLocalOperator<TypeTag>
{
    ...
}
```
<span style="font-size: 0.4em; position: relative; top: -38px; color: gray;">File: `dumux/examples/diffusion/model.hh`</span>

Inherits from the `DiscretizationDefaultLocalOperator`, 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;
 }
```
<span style="font-size: 0.4em; position: relative; top: -38px; color: gray;">File: `dumux/examples/diffusion/model.hh`</span>

## Flux term

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

with

- $c^n_h:$ global discrete solution at time $t_n$, interpolated using __basis functions__
- $\mathbf{n}:$ unit outer normal vector
- $\sigma:$ sub control volume face (scvf)

## Flux term
<font size=8>
```cpp
NumEqVector computeFlux(const Problem& problem,
                        const Element& element,
                        const FVElementGeometry& fvGeometry,
                        const ElementVolumeVariables& elemVolVars,
                        const SubControlVolumeFace& scvf,
                        const ElementFluxVariablesCache& elemFluxVarsCache) const
{
    ...
}
```
</font>
<span style="font-size: 0.4em; position: relative; top: -60px; color: gray;">File: `dumux/examples/diffusion/model.hh`</span>

## 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())
        );
    }
    ...
}
```
<span style="font-size: 0.4em; position: relative; top: -38px; color: gray;">File: `dumux/examples/diffusion/model.hh`</span>

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

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

    return flux;
}
```
<span style="font-size: 0.4em; position: relative; top: -38px; color: gray;">File: `dumux/examples/diffusion/model.hh`</span>

## `LocalResidual`

A `LocalResidual` implements the discretized mathematical model.

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

# Model properties

## Model properties

All properties are defined within the namespace `Dumux::Properties`.

```cpp
namespace Dumux::Properties {
    //define all properties
} // end namespace Dumux::Properties
```
<span style="font-size: 0.4em; position: relative; top: -38px; color: gray;">File: `dumux/examples/diffusion/model.hh`</span>

## Model type tag

The property type tag is an empty struct with the respective name, e.g. `DiffusionModel`.
All properties related to the desired model can be attributed using the property type tag.

```cpp
namespace Dumux::Properties::TTag {
//! The diffusion model tag that we can specialize properties for
struct DiffusionModel {};
} // end namespace Dumux::Properties::TTag
```
<span style="font-size: 0.4em; position: relative; top: -38px; color: gray;">File: `dumux/examples/diffusion/model.hh`</span>

## 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>; };
```
<span style="font-size: 0.4em; position: relative; top: -38px; color: gray;">File: `dumux/examples/diffusion/model.hh`</span>

## 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; }
    };
};
```
<span style="font-size: 0.4em; position: relative; top: -38px; color: gray;">File: `dumux/examples/diffusion/model.hh`</span>

## 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** for denoising of an MRI image

<img src="img/exercise_model_mri_denoise.gif" alt="denoising" width="300"/>

## Tasks

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

* Go to [Model exercise](https://git.iws.uni-stuttgart.de/dumux-repositories/dumux-course/tree/master/exercises/exercise-model#exercise-model-dumux-course)