From c515a261161d80a219f1405d4f16f011e4c5ebab Mon Sep 17 00:00:00 2001
From: Timo Koch <timokoch@math.uio.no>
Date: Fri, 24 Mar 2023 17:42:11 +0100
Subject: [PATCH] [example][cahnhilliard] Further expand and improve example
description
---
examples/cahn_hilliard/.doc_config | 16 +-
examples/cahn_hilliard/CMakeLists.txt | 1 +
examples/cahn_hilliard/README.md | 131 ++++--
examples/cahn_hilliard/doc/_intro.md | 127 ++++--
.../doc/{mainfile_intro.md => _main.md} | 0
.../doc/{modelheader_intro.md => _model.md} | 0
.../doc/{mainfile.md => main.md} | 374 +++++++----------
.../doc/{modelheader.md => model.md} | 311 +++++++--------
examples/cahn_hilliard/main.cc | 376 +++++++-----------
examples/cahn_hilliard/model.hh | 307 +++++++-------
examples/diffusion/doc/main.md | 4 +-
examples/diffusion/main.cc | 4 +-
12 files changed, 757 insertions(+), 894 deletions(-)
rename examples/cahn_hilliard/doc/{mainfile_intro.md => _main.md} (100%)
rename examples/cahn_hilliard/doc/{modelheader_intro.md => _model.md} (100%)
rename examples/cahn_hilliard/doc/{mainfile.md => main.md} (53%)
rename examples/cahn_hilliard/doc/{modelheader.md => model.md} (50%)
diff --git a/examples/cahn_hilliard/.doc_config b/examples/cahn_hilliard/.doc_config
index d1a87cc8a3..6df89661cd 100644
--- a/examples/cahn_hilliard/.doc_config
+++ b/examples/cahn_hilliard/.doc_config
@@ -3,21 +3,21 @@
"doc/_intro.md"
],
- "doc/mainfile.md" : [
- "doc/mainfile_intro.md",
- "main.cc"
+ "doc/model.md" : [
+ "doc/_model.md",
+ "model.hh"
],
- "doc/modelheader.md" : [
- "doc/modelheader_intro.md",
- "model.hh"
+ "doc/main.md" : [
+ "doc/_main.md",
+ "main.cc"
],
"navigation" : {
"mainpage" : "README.md",
"subpages" : [
- "doc/modelheader.md",
- "doc/mainfile.md"
+ "doc/model.md",
+ "doc/main.md"
],
"subtitles" : [
"Model implementation",
diff --git a/examples/cahn_hilliard/CMakeLists.txt b/examples/cahn_hilliard/CMakeLists.txt
index 73113b9863..f69a567fda 100644
--- a/examples/cahn_hilliard/CMakeLists.txt
+++ b/examples/cahn_hilliard/CMakeLists.txt
@@ -14,6 +14,7 @@ dumux_add_test(NAME example_cahn_hilliard
dumux_add_test(NAME example_cahn_hilliard_parallel
LABELS example
TARGET example_cahn_hilliard
+ CMAKE_GUARD MPI_FOUND
COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
CMD_ARGS --script fuzzy
--files ${CMAKE_SOURCE_DIR}/test/references/example_cahn_hilliard_parallel_p0-reference.vtu
diff --git a/examples/cahn_hilliard/README.md b/examples/cahn_hilliard/README.md
index a192a3346b..11552d3711 100644
--- a/examples/cahn_hilliard/README.md
+++ b/examples/cahn_hilliard/README.md
@@ -1,70 +1,127 @@
<!-- Important: This file has been automatically generated by generate_example_docs.py. Do not edit this file directly! -->
-# Cahn-Hilliard model
+# Cahn-Hilliard equation
-A random initial distribution of two phases separating according to the Cahn-Hilliard model.
+In this example, we create an application solving the Cahn-Hilliard equation which describes
+a phase separation process.
+This example is implemented in three files, [`model.hh`](model.hh), [`params.input`](params.input), [`main.cc`](main.cc).
+The [`model.hh`](model.hh) header implements a simple self-contained Cahn-Hilliard equation model
+to be used with a control-volume finite element discretization like the Box method.
+The [`main.cc`](main.cc) source file contains the main program and [`params.input`](params.input) contains
+runtime parameters. The executable is configured in [`CMakeLists.txt`](CMakeLists.txt) and created with CMake.
+The source code will be discussed in Part 1 and Part 2 after an introduction to the equations and the problem description.
__The main points illustrated in this example are__
-* A base setup for solving a nonlinear partial differential equation
-__Table of contents__. This description is structured as follows:
-
-[[_TOC_]]
+* setting up a model with two coupled nonlinear equations
+* creating a vector-valued random initial solution in parallel
+* solving a time-dependent nonlinear problem
+* using a Newton solver instance
+* enabling partial reassembly for the Newton method
+* using adaptive time stepping
-__Result__. The result will look like below.
-You see the concentration variable $c$ capturing the phase distribution after 1 second.
+__Results__. The result will look like below.
<figure>
<center>
- <img src="img/results_phase_distribution.png" alt="Concentration field modeling phase distribution" width="50%"/>
+ <img src="img/animation.gif" alt="Cahn Hilliard result" width="400px"/>
+ <figcaption> <b> Fig.1 </b> - Evolution of an initially random concentration field according to the Cahn-Hilliard equation.</figcaption>
</center>
</figure>
-## Problem set-up
+__Table of contents__. This description is structured as follows:
-A square two-dimensional domain is initialized with a concentration `c`, randomly perturbed by
-noise following a uniform distribution between 0.41 and 0.43.
-With time the concentration field evolves towards attaining mostly values near to 0 or 1 while
-conserving the total concentration, modeling the separation of two immiscible fluids.
+[TOC]
-## Model description
+## Equation and problem description
-The Cahn-Hilliard model uses a pair of second order nonlinear partial differential equations,
-for a concentration field $c$ and a chemical potential $\mu$. The former is conserved and
-evolves with a flux proportional to $\nabla \mu$, while the latter depends on the laplacian of
-the concentration $\nabla^2 c$ and a nonlinear function of $c$, a derivative of a free energy
-functional.
+The Cahn-Hilliard equation is a forth order partial differential equation (PDE) and reads
```math
-\frac{\partial c}{\partial t} = M \nabla^2 \mu \\
-\mu = E (4 c^3 - 6 c^2 +2 c) - \gamma \nabla^2 c
+\frac{\partial c}{\partial t} - \nabla\cdot \left( M \nabla\left[ \frac{\partial f}{\partial c} - \gamma \nabla^2 c \right]\right) = 0 \quad \text{in} \; \Omega \times (0, T]
```
-Here $M$ denotes the mobility of the concentration, while the energy scale $E$ and surface tension
-$\gamma$ balance the two contributions of the free energy functional.
+with the concentration $c(x,t)$, the mobility coefficient $M$, and surface tension $\gamma$.
+The domain $\Omega \subset \mathbb{R}^2$ is initialized with a concentration field
+$c(x,t=0) = 0.42 + \zeta$, randomly perturbed by
+noise $\zeta$ following a uniform distribution $\zeta \sim U(-0.02, 0.02)$.
+With time the concentration field evolves towards attaining mostly values near to $0$ or $1$ while
+conserving the total concentration. The model describes the separation of two immiscible fluids.
-## Implementation of a simple nonlinear PDE
+The fourth order PDE cannot be solved by a standard finite volume scheme. We therefore
+split the equation in two second order PDEs
-For this example the C++ code is contained in two files, `main.cc` and `model.hh`.
-The `model.hh` header contains the `ModelTraits` and properties for the model TypeTag
-`CahnHilliardModel`, as well as volumevariables `CahnHilliardModelVolumeVariables`
-and the basic local residual `CahnHilliardModelLocalResidual`.
-The residual's storage term is extended in the problem definition `CahnHilliardTestProblem`
-in the file `main.cc`, which also contains more specific properties of the problem setup (for
-type tag `CahnHilliardTest`) as well as the actual simulation loop usually found in a mainfile.
+```math
+\begin{aligned}
+\frac{\partial c}{\partial t} - \nabla\cdot \left( M \nabla \mu \right) &= 0\\
+- \gamma \nabla^2 c &= \mu - \frac{\partial f}{\partial c}.
+\end{aligned}
+```
+
+were $\mu$ is called chemical potential.
+Using the free energy functional $f = E c^2(1-c)^2$, where $E$ is a scaling parameter,
+we obtain
+
+```math
+\begin{aligned}
+\frac{\partial c}{\partial t} - \nabla\cdot \left( M \nabla \mu \right) &= 0\\
+- \gamma \nabla^2 c &= \mu - E (4 c^3 - 6 c^2 +2 c).
+\end{aligned}
+```
+
+The concentration field is conserved and
+evolves with a flux proportional to $\nabla \mu$, while the latter depends on the Laplacian of
+the concentration $\nabla^2 c$ and the free energy functional which is a nonlinear function of $c$. We therefore have a system of two coupled nonlinear PDEs. We will use homogeneous Neumann
+boundary conditions everywhere on $\partial \Omega$. This means that the total concentration in $\Omega$ stays constant.
+
+For the implementation, it is useful to write the conservation equation in the following abstract
+form in vector notation,
-More details are given in [main.cc](doc/mainfile.md) and [model.hh](doc/modelheader.md).
+```math
+\frac{\partial \boldsymbol{S}(\boldsymbol{U})}{\partial t} + \nabla\cdot \boldsymbol{F}(\boldsymbol{U}) = \boldsymbol{Q}(\boldsymbol{U})
+```
+
+where for the Cahn-Hillard equation, we have
+
+```math
+\boldsymbol{U} := \begin{bmatrix} c \\ \mu \end{bmatrix}, \quad
+\boldsymbol{S}(\boldsymbol{U}) := \begin{bmatrix} c \\ 0 \end{bmatrix}, \quad
+\boldsymbol{F}(\boldsymbol{U}) := \begin{bmatrix} -M \nabla \mu \\ -\gamma \nabla c \end{bmatrix}, \quad
+\boldsymbol{Q}(\boldsymbol{U}) := \begin{bmatrix} 0 \\ \mu - E (4 c^3 - 6 c^2 +2 c) \end{bmatrix}.
+```
-Additionally the folder contains `params.input` containing runtime parameters and the build
-system file `CMakeLists.txt` defining conditions for the accompanying test.
+### Model parameters
+
+* $\Omega = [0,1]\times[0,1]$ (unit square)
+* End time $T = 1$, initial time step size $\Delta t = 0.0025$, and maximum time step size $\Delta t = 0.01$
+* $100 \times 100$ cells (structured Cartesian grid)
+* $M = 0.0001$
+* $\gamma = 0.01$
+* $f = E c^2(1-c)^2$, with $E = 100$
+* $c^0_{h,B} = 0.42 + \zeta$, where $\zeta \sim U(-0.02,0.02)$
+* $\mu^0_{h,B} = 0.0$
+* homogeneous Neumann boundary conditions
+
+### Running the example in parallel
+
+See [Diffusion equation example](https://git.iws.uni-stuttgart.de/dumux-repositories/dumux/-/blob/master/examples/diffusion/README.md).
+
+# Implementation
+
+For this example the C++ code is contained in two files, [`model.hh`](model.hh) and [`main.cc`](main.cc). The [`model.hh`](model.hh) header contains the `ModelTraits` and properties for the
+model type tag `CahnHilliardModel`, as well as the volume variables class
+`CahnHilliardModelVolumeVariables` and the local residual class `CahnHilliardModelLocalResidual`.
+The source term is extended in the problem definition `CahnHilliardTestProblem`
+in the file `main.cc`, which also contains more specific properties of the problem setup (for
+type tag `CahnHilliardTest`) as well as the actual simulation loop usually found in the main function.
## Part 1: Model implementation
-| [:arrow_right: Click to continue with part 1 of the documentation](doc/modelheader.md) |
+| [:arrow_right: Click to continue with part 1 of the documentation](doc/model.md) |
|---:|
## Part 2: Main program
-| [:arrow_right: Click to continue with part 2 of the documentation](doc/mainfile.md) |
+| [:arrow_right: Click to continue with part 2 of the documentation](doc/main.md) |
|---:|
diff --git a/examples/cahn_hilliard/doc/_intro.md b/examples/cahn_hilliard/doc/_intro.md
index 0c63dab293..42f94f6c00 100644
--- a/examples/cahn_hilliard/doc/_intro.md
+++ b/examples/cahn_hilliard/doc/_intro.md
@@ -1,57 +1,114 @@
-# Cahn-Hilliard model
+# Cahn-Hilliard equation
-A random initial distribution of two phases separating according to the Cahn-Hilliard model.
+In this example, we create an application solving the Cahn-Hilliard equation which describes
+a phase separation process.
+This example is implemented in three files, [`model.hh`](model.hh), [`params.input`](params.input), [`main.cc`](main.cc).
+The [`model.hh`](model.hh) header implements a simple self-contained Cahn-Hilliard equation model
+to be used with a control-volume finite element discretization like the Box method.
+The [`main.cc`](main.cc) source file contains the main program and [`params.input`](params.input) contains
+runtime parameters. The executable is configured in [`CMakeLists.txt`](CMakeLists.txt) and created with CMake.
+The source code will be discussed in Part 1 and Part 2 after an introduction to the equations and the problem description.
__The main points illustrated in this example are__
-* A base setup for solving a nonlinear partial differential equation
-__Table of contents__. This description is structured as follows:
-
-[[_TOC_]]
+* setting up a model with two coupled nonlinear equations
+* creating a vector-valued random initial solution in parallel
+* solving a time-dependent nonlinear problem
+* using a Newton solver instance
+* enabling partial reassembly for the Newton method
+* using adaptive time stepping
-__Result__. The result will look like below.
-You see the concentration variable $c$ capturing the phase distribution after 1 second.
+__Results__. The result will look like below.
<figure>
<center>
- <img src="img/results_phase_distribution.png" alt="Concentration field modeling phase distribution" width="50%"/>
+ <img src="img/animation.gif" alt="Cahn Hilliard result" width="400px"/>
+ <figcaption> <b> Fig.1 </b> - Evolution of an initially random concentration field according to the Cahn-Hilliard equation.</figcaption>
</center>
</figure>
-## Problem set-up
+__Table of contents__. This description is structured as follows:
-A square two-dimensional domain is initialized with a concentration `c`, randomly perturbed by
-noise following a uniform distribution between 0.41 and 0.43.
-With time the concentration field evolves towards attaining mostly values near to 0 or 1 while
-conserving the total concentration, modeling the separation of two immiscible fluids.
+[TOC]
-## Model description
+## Equation and problem description
-The Cahn-Hilliard model uses a pair of second order nonlinear partial differential equations,
-for a concentration field $c$ and a chemical potential $\mu$. The former is conserved and
-evolves with a flux proportional to $\nabla \mu$, while the latter depends on the laplacian of
-the concentration $\nabla^2 c$ and a nonlinear function of $c$, a derivative of a free energy
-functional.
+The Cahn-Hilliard equation is a forth order partial differential equation (PDE) and reads
```math
-\frac{\partial c}{\partial t} = M \nabla^2 \mu \\
-\mu = E (4 c^3 - 6 c^2 +2 c) - \gamma \nabla^2 c
+\frac{\partial c}{\partial t} - \nabla\cdot \left( M \nabla\left[ \frac{\partial f}{\partial c} - \gamma \nabla^2 c \right]\right) = 0 \quad \text{in} \; \Omega \times (0, T]
```
-Here $M$ denotes the mobility of the concentration, while the energy scale $E$ and surface tension
-$\gamma$ balance the two contributions of the free energy functional.
+with the concentration $c(x,t)$, the mobility coefficient $M$, and surface tension $\gamma$.
+The domain $\Omega \subset \mathbb{R}^2$ is initialized with a concentration field
+$c(x,t=0) = 0.42 + \zeta$, randomly perturbed by
+noise $\zeta$ following a uniform distribution $\zeta \sim U(-0.02, 0.02)$.
+With time the concentration field evolves towards attaining mostly values near to $0$ or $1$ while
+conserving the total concentration. The model describes the separation of two immiscible fluids.
-## Implementation of a simple nonlinear PDE
+The fourth order PDE cannot be solved by a standard finite volume scheme. We therefore
+split the equation in two second order PDEs
-For this example the C++ code is contained in two files, `main.cc` and `model.hh`.
-The `model.hh` header contains the `ModelTraits` and properties for the model TypeTag
-`CahnHilliardModel`, as well as volumevariables `CahnHilliardModelVolumeVariables`
-and the basic local residual `CahnHilliardModelLocalResidual`.
-The residual's storage term is extended in the problem definition `CahnHilliardTestProblem`
-in the file `main.cc`, which also contains more specific properties of the problem setup (for
-type tag `CahnHilliardTest`) as well as the actual simulation loop usually found in a mainfile.
+```math
+\begin{aligned}
+\frac{\partial c}{\partial t} - \nabla\cdot \left( M \nabla \mu \right) &= 0\\
+- \gamma \nabla^2 c &= \mu - \frac{\partial f}{\partial c}.
+\end{aligned}
+```
+
+were $\mu$ is called chemical potential.
+Using the free energy functional $f = E c^2(1-c)^2$, where $E$ is a scaling parameter,
+we obtain
+
+```math
+\begin{aligned}
+\frac{\partial c}{\partial t} - \nabla\cdot \left( M \nabla \mu \right) &= 0\\
+- \gamma \nabla^2 c &= \mu - E (4 c^3 - 6 c^2 +2 c).
+\end{aligned}
+```
+
+The concentration field is conserved and
+evolves with a flux proportional to $\nabla \mu$, while the latter depends on the Laplacian of
+the concentration $\nabla^2 c$ and the free energy functional which is a nonlinear function of $c$. We therefore have a system of two coupled nonlinear PDEs. We will use homogeneous Neumann
+boundary conditions everywhere on $\partial \Omega$. This means that the total concentration in $\Omega$ stays constant.
+
+For the implementation, it is useful to write the conservation equation in the following abstract
+form in vector notation,
-More details are given in [main.cc](doc/mainfile.md) and [model.hh](doc/modelheader.md).
+```math
+\frac{\partial \boldsymbol{S}(\boldsymbol{U})}{\partial t} + \nabla\cdot \boldsymbol{F}(\boldsymbol{U}) = \boldsymbol{Q}(\boldsymbol{U})
+```
+
+where for the Cahn-Hillard equation, we have
+
+```math
+\boldsymbol{U} := \begin{bmatrix} c \\ \mu \end{bmatrix}, \quad
+\boldsymbol{S}(\boldsymbol{U}) := \begin{bmatrix} c \\ 0 \end{bmatrix}, \quad
+\boldsymbol{F}(\boldsymbol{U}) := \begin{bmatrix} -M \nabla \mu \\ -\gamma \nabla c \end{bmatrix}, \quad
+\boldsymbol{Q}(\boldsymbol{U}) := \begin{bmatrix} 0 \\ \mu - E (4 c^3 - 6 c^2 +2 c) \end{bmatrix}.
+```
-Additionally the folder contains `params.input` containing runtime parameters and the build
-system file `CMakeLists.txt` defining conditions for the accompanying test.
+### Model parameters
+
+* $\Omega = [0,1]\times[0,1]$ (unit square)
+* End time $T = 1$, initial time step size $\Delta t = 0.0025$, and maximum time step size $\Delta t = 0.01$
+* $100 \times 100$ cells (structured Cartesian grid)
+* $M = 0.0001$
+* $\gamma = 0.01$
+* $f = E c^2(1-c)^2$, with $E = 100$
+* $c^0_{h,B} = 0.42 + \zeta$, where $\zeta \sim U(-0.02,0.02)$
+* $\mu^0_{h,B} = 0.0$
+* homogeneous Neumann boundary conditions
+
+### Running the example in parallel
+
+See [Diffusion equation example](https://git.iws.uni-stuttgart.de/dumux-repositories/dumux/-/blob/master/examples/diffusion/README.md).
+
+# Implementation
+
+For this example the C++ code is contained in two files, [`model.hh`](model.hh) and [`main.cc`](main.cc). The [`model.hh`](model.hh) header contains the `ModelTraits` and properties for the
+model type tag `CahnHilliardModel`, as well as the volume variables class
+`CahnHilliardModelVolumeVariables` and the local residual class `CahnHilliardModelLocalResidual`.
+The source term is extended in the problem definition `CahnHilliardTestProblem`
+in the file `main.cc`, which also contains more specific properties of the problem setup (for
+type tag `CahnHilliardTest`) as well as the actual simulation loop usually found in the main function.
diff --git a/examples/cahn_hilliard/doc/mainfile_intro.md b/examples/cahn_hilliard/doc/_main.md
similarity index 100%
rename from examples/cahn_hilliard/doc/mainfile_intro.md
rename to examples/cahn_hilliard/doc/_main.md
diff --git a/examples/cahn_hilliard/doc/modelheader_intro.md b/examples/cahn_hilliard/doc/_model.md
similarity index 100%
rename from examples/cahn_hilliard/doc/modelheader_intro.md
rename to examples/cahn_hilliard/doc/_model.md
diff --git a/examples/cahn_hilliard/doc/mainfile.md b/examples/cahn_hilliard/doc/main.md
similarity index 53%
rename from examples/cahn_hilliard/doc/mainfile.md
rename to examples/cahn_hilliard/doc/main.md
index ef55ebc606..778b01e761 100644
--- a/examples/cahn_hilliard/doc/mainfile.md
+++ b/examples/cahn_hilliard/doc/main.md
@@ -1,60 +1,78 @@
<!-- Important: This file has been automatically generated by generate_example_docs.py. Do not edit this file directly! -->
-| [:arrow_left: Back to the main documentation](../README.md) | [:arrow_left: Go back to part 1](modelheader.md) |
+| [:arrow_left: Back to the main documentation](../README.md) | [:arrow_left: Go back to part 1](model.md) |
|---|---:|
-# Problem, test properties/traits and main program flow (`main.cc`)
-In this example the file `main.cc` contains the problem class `CahnHilliardTestProblem`,
-properties and traits specific to the test case as well as the main program flow in the form of
-`main` function.
+# Cahn-Hilliard equation main program
+The file [`main.cc`](main.cc) contains the problem class `CahnHilliardTestProblem`,
+properties and traits specific to the test case as well as the `main` function.
+The setup consists of four steps:
+1. Define the problem setting boundary conditions and the diffusion coefficient
+2. Configure the property system reusing the model defined in Part 1
+3. Define a function for setting the random initial condition
+4. The main program defining all steps of the program
-```cpp
-#include <config.h>
-```
+__Table of contents__
-## Problem
-The __problem class__ defines boundary conditions and extends the storage term defined in the
-model's localresidual by the derivative of the free energy.
+[TOC]
+We start in `main.cc` with the necessary header includes:
+<details><summary> Click to show includes</summary>
-<details open>
-<summary><b>Click to hide/show the file documentation</b> (or inspect the [source code](../main.cc))</summary>
-
-### Include headers
+```cpp
+#include <config.h>
+#include <type_traits>
+#include <random>
-```cpp
-// use the property system and runtime parameters
+#include <dumux/common/initialize.hh>
#include <dumux/common/properties.hh>
#include <dumux/common/parameters.hh>
-// common DuMux vector for discretized equations
#include <dumux/common/numeqvector.hh>
-// types of boundary conditions to use
#include <dumux/common/boundarytypes.hh>
-// generic problem for finite volume simulations
#include <dumux/common/fvproblem.hh>
+
+#include <dumux/discretization/box.hh>
+#include <dumux/linear/linearsolvertraits.hh>
+#include <dumux/linear/linearalgebratraits.hh>
+#include <dumux/linear/istlsolvers.hh>
+#include <dumux/nonlinear/newtonsolver.hh>
+#include <dumux/assembly/fvassembler.hh>
+
+#include <dune/grid/yaspgrid.hh>
+#include <dumux/io/vtkoutputmodule.hh>
+#include <dumux/io/grid/gridmanager_yasp.hh>
+
+#include "model.hh"
```
+</details>
+
+## 1. The problem class
-### The problem class `CahnHilliardTestProblem`
-In this abstract problem we inherit from the generic `FVProblem`.
+The problem class implements boundary conditions and the source term.
+It also provides an interface that is used by the local residual (see Part 1) to obtain
+the mobility and surface tension. The values are read from the parameter configuration tree.
+
+<details open>
+<summary><b>Click to hide/show the file documentation</b> (or inspect the [source code](../main.cc))</summary>
```cpp
namespace Dumux {
-
template<class TypeTag>
class CahnHilliardTestProblem : public FVProblem<TypeTag>
{
```
-<details><summary> Click to show alias definitions and local variables</summary>
+<details><summary> Click to show boilerplate code</summary>
```cpp
+
using ParentType = FVProblem<TypeTag>;
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
using FVElementGeometry = typename GridGeometry::LocalView;
@@ -71,6 +89,8 @@ class CahnHilliardTestProblem : public FVProblem<TypeTag>
```
</details>
+In the constructor, we read the parameter constants from the
+parameter tree (which is initialized with the content of `params.input`).
```cpp
public:
@@ -83,13 +103,9 @@ public:
}
```
-
-### Problem source term
-
-Here we implement the derivative of the free energy, setting a source for the equation for
-the chemical potential. The `computeSource` function in the local residual adds the terms
-defined here.
-
+In the `source` function, we implement the derivative of the free energy.
+This demonstrates how parts of the local residual can be split into model specific
+parts and parts that might change from scenario to scenario.
```cpp
template<class ElementVolumeVariables>
@@ -105,13 +121,9 @@ defined here.
}
```
-
-### Boundary conditions
-
For the boundary we choose boundary flux (or Neumann) conditions for all equations and on
every part of the boundary, specifying zero flux everywhere for both equations.
-
```cpp
BoundaryTypes boundaryTypesAtPos(const GlobalPosition& globalPos) const
{
@@ -124,11 +136,9 @@ every part of the boundary, specifying zero flux everywhere for both equations.
{ return { 0.0, 0.0 }; }
```
-
-<details><summary> Click to show coefficients and private variables</summary>
-The problem class offers access to the mobility and surface tension coefficients as read from
-the parameter file by default `params.input`.
-
+The parameters interfaces are used in the local residual (see Part 1).
+We can name this interface however we want as long as we adapt the calling site
+in the `LocalResidual` class in `model.hh`.
```cpp
Scalar mobility() const
@@ -142,136 +152,61 @@ private:
Scalar surfaceTension_;
Scalar energyScale_;
};
-
} // end namespace Dumux
```
-</details>
</details>
-## Test case properties/traits
-Within the `Dumux::Properties` namespace we specialize properties and traits to the considered
-test case by using the test's TypeTag.
+## 2. Property tag and specializations
+
+We create a new tag `DiffusionTest` that inherits all properties
+specialized for the tag `DiffusionModel` (we created this in Part 1)
+and the tag `BoxModel` which provides types relevant for the spatial
+discretization scheme (see [dumux/discretization/box.hh](https://git.iws.uni-stuttgart.de/dumux-repositories/dumux/-/blob/master/dumux/discretization/box.hh)).
<details open>
<summary><b>Click to hide/show the file documentation</b> (or inspect the [source code](../main.cc))</summary>
-### Include headers
-
-
-```cpp
-// Include the grid to be used
-#include <dune/grid/yaspgrid.hh>
-// The header for the box discretization scheme
-#include <dumux/discretization/box.hh>
-// The model header including the model traits and properties
-#include "model.hh"
-```
-
-
-### TypeTag `CahnHilliardTest`
-
-We define a type tag for the test case, allowing us to further specify properties and traits. To
-use those set for the Cahn-Hilliard model we inherit from its type tag.
-
-
```cpp
-namespace Dumux::Properties {
-
-// Inheriting properties of the Cahn-Hilliard model and the box finite volume discretization
-namespace TTag {
-struct CahnHilliardTest { using InheritsFrom = std::tuple<CahnHilliardModel, BoxModel>; };
-} // end namespace TTag
-```
-
-
-### Test properties
-
-We specify a grid to be used in the test, select our problem class and enable caching.
-
-
-```cpp
-// Set the grid type
-template<class TypeTag>
-struct Grid<TypeTag, TTag::CahnHilliardTest>
-{ using type = Dune::YaspGrid<2>; };
-
-// Select the problem class defined above
-template<class TypeTag>
-struct Problem<TypeTag, TTag::CahnHilliardTest>
-{ using type = CahnHilliardTestProblem<TypeTag>; };
+namespace Dumux::Properties::TTag {
+struct CahnHilliardTest
+{
+ using InheritsFrom = std::tuple<CahnHilliardModel, BoxModel>;
-// Enable caching
-template<class TypeTag>
-struct EnableGridVolumeVariablesCache<TypeTag, TTag::CahnHilliardTest>
-{ static constexpr bool value = true; };
+ using Grid = Dune::YaspGrid<2>;
-template<class TypeTag>
-struct EnableGridFluxVariablesCache<TypeTag, TTag::CahnHilliardTest>
-{ static constexpr bool value = true; };
-
-template<class TypeTag>
-struct EnableGridGeometryCache<TypeTag, TTag::CahnHilliardTest>
-{ static constexpr bool value = true; };
+ template<class TypeTag>
+ using Problem = CahnHilliardTestProblem<TypeTag>;
+ using EnableGridGeometryCache = std::true_type;
+ using EnableGridVolumeVariablesCache = std::true_type;
+ using EnableGridFluxVariablesCache = std::true_type;
+};
} // end namespace Dumux::Properties
```
</details>
-## The main program flow
-The main program flow in the `main` function is often the only content of `main.cc`. It sets up
-the simulation framework, initializes runtime parameters, creates the grid and storage vectors
-for the variables, primary and secondary. It specifies and constructs and assembler, which
-assembles the discretized residual and system matrix (Jacobian of the model residual), as well as
-linear and nonlinear solvers that solve the resulting linear system and handle the convergence of
-nonlinear iterations. The time loop controls the time stepping, with adaptive time step size in
-coordination with the nonlinear solver.
+## 3. Creating the initial solution
+
+To initialize with a random field in parallel, where each processor
+creates its own random number sequence, we need to communicate the
+resulting values on the processor border and overlap.
+See [Diffusion equation example](https://git.iws.uni-stuttgart.de/dumux-repositories/dumux/-/blob/master/examples/diffusion/docs/main.md).
+for details. Here in addition, we need to provide a custom scatter operation
+that handles vector types. We only need to communicate the first entry (concentration).
<details open>
<summary><b>Click to hide/show the file documentation</b> (or inspect the [source code](../main.cc))</summary>
-### Include headers
-
-
-```cpp
-// standard header to generate random initial data
-#include <random>
-// common DuMux header for parallelization
-#include <dumux/common/initialize.hh>
-// headers to use the property system and runtime parameters
-#include <dumux/common/properties.hh>
-#include <dumux/common/parameters.hh>
-// module for VTK output, to write out fields of interest
-#include <dumux/io/vtkoutputmodule.hh>
-// gridmanager for the grid used in the test
-#include <dumux/io/grid/gridmanager_yasp.hh>
-// headers for linear and non-linear solvers as well as the assembler
-#include <dumux/linear/linearsolvertraits.hh>
-#include <dumux/linear/linearalgebratraits.hh>
-#include <dumux/linear/istlsolvers.hh>
-#include <dumux/nonlinear/newtonsolver.hh>
-#include <dumux/assembly/fvassembler.hh>
-```
-
-
-### Creating the initial solution
-
-We define a helper struct and function to handle communication for parallel runs.
-For our initial conditions we create a random field of values around a mean of 0.42.
-The random values are created with an offset based on the processor rank for communication
-purposes, which is removed afterwards. For more information see the description of the
-diffusion example
-[examples/diffusion/doc/main.md](https://git.iws.uni-stuttgart.de/dumux-repositories/dumux/-/blob/master/examples/diffusion/doc/main.md)
-
-
```cpp
+// functor for data communication with MPI
struct MinScatter
{
template<class A, class B>
@@ -279,6 +214,7 @@ struct MinScatter
{ a[0] = std::min(a[0], b[0]); }
};
+// create the random initial solution
template<class SolutionVector, class GridGeometry>
SolutionVector createInitialSolution(const GridGeometry& gg)
{
@@ -294,7 +230,8 @@ SolutionVector createInitialSolution(const GridGeometry& gg)
sol[n][1] = 0.0;
}
- // Take the value of the processor with the minimum rank and subtract the rank offset
+ // We take the value of the processor with the minimum rank
+ // and subtract the rank offset
if (gg.gridView().comm().size() > 1)
{
Dumux::VectorCommDataHandle<
@@ -314,143 +251,122 @@ SolutionVector createInitialSolution(const GridGeometry& gg)
```
-### The main function
+</details>
-The main function takes the command line arguments, optionally specifying an input file of
-parameters and/or individual key-value pairs of runtime parameters.
+## 4. The main program
+
+The `main` function sets up the simulation framework, initializes runtime parameters,
+creates the grid and storage vectors
+for the variables, primary and secondary. It specifies and constructs and assembler, which
+assembles the discretized residual and system matrix (Jacobian of the model residual), as well as
+a linear solver. A Newton method is used to solve the nonlinear equations where in each Newton iteration
+the Jacobian and the residual needs to be reassembled and the resulting linear system is solved.
+The time loop controls the time stepping. Here, we let the Newton solver suggest an adaption of
+the time step size based on the convergence history of the nonlinear solver.
+
+
+<details open>
+<summary><b>Click to hide/show the file documentation</b> (or inspect the [source code](../main.cc))</summary>
```cpp
int main(int argc, char** argv)
{
using namespace Dumux;
-
- // define the type tag for this problem
- using TypeTag = Properties::TTag::CahnHilliardTest;
```
-
-We initialize parallelization backends as well as runtime parameters
-
+We initialize MPI and/or multithreading backend. When not running
+in parallel the MPI setup is skipped.
```cpp
- // maybe initialize MPI and/or multithreading backend
Dumux::initialize(argc, argv);
-
- // initialize parameter tree
- Parameters::init(argc, argv);
```
-
-### Grid setup
-
-Set up the grid as well as a grid geometry to access the (sub-)control-volumes and their
-faces
-
+Then we initialize parameter tree.
```cpp
- // initialize the grid
- GridManager<GetPropType<TypeTag, Properties::Grid>> gridManager;
- gridManager.init();
-
- // we compute on the leaf grid view
- const auto& leafGridView = gridManager.grid().leafGridView();
-
- // create the finite volume grid geometry
- using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
- auto gridGeometry = std::make_shared<GridGeometry>(leafGridView);
+ Parameters::init(argc, argv);
```
-
-### Problem setup
-
-We instantiate also the problem class according to the test properties
-
+We create an alias for the type tag for this problem
+and extract type information through properties.
```cpp
- // the problem (initial and boundary conditions)
+ using TypeTag = Properties::TTag::CahnHilliardTest;
+
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using Grid = GetPropType<TypeTag, Properties::Grid>;
+ using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
using Problem = GetPropType<TypeTag, Properties::Problem>;
- auto problem = std::make_shared<Problem>(gridGeometry);
+ using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
+ using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
```
-
-### Applying initial conditions
-
-After writing the initial data to the storage for previous and current time-step, we
-initialize the grid variables, also computing secondary variables.
-
+We initialize the grid, grid geometry, problem, solution vector, and grid variables.
```cpp
- // the solution vector
- using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
- auto sol = createInitialSolution<SolutionVector>(*gridGeometry);
- // copy the vector to store state of previous time step
- auto oldSol = sol;
+ GridManager<Grid> gridManager;
+ gridManager.init();
- // the grid variables
- using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
+ auto gridGeometry = std::make_shared<GridGeometry>(gridManager.grid().leafGridView());
+ auto problem = std::make_shared<Problem>(gridGeometry);
+ auto sol = createInitialSolution<SolutionVector>(*gridGeometry);
auto gridVariables = std::make_shared<GridVariables>(problem, gridGeometry);
gridVariables->init(sol);
```
-
-### Initialize VTK output
-
+We initialize the VTK output module and write out the initial concentration and chemical potential
```cpp
- // initialize the vtk output module
VtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables, sol, problem->name());
vtkWriter.addVolumeVariable([](const auto& vv){ return vv.concentration(); }, "c");
vtkWriter.addVolumeVariable([](const auto& vv){ return vv.chemicalPotential(); }, "mu");
vtkWriter.write(0.0);
```
-
-### Set up time loop
-
+We instantiate time loop using start and end time as well as
+the time step size from the parameter tree (`params.input`)
```cpp
- // get some time loop parameters
- using Scalar = GetPropType<TypeTag, Properties::Scalar>;
- const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
- const auto dt = getParam<Scalar>("TimeLoop.InitialTimeStepSize");
- const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
-
- // instantiate time loop
- auto timeLoop = std::make_shared<CheckPointTimeLoop<Scalar>>(0.0, dt, tEnd);
- timeLoop->setMaxTimeStepSize(maxDt);
+ auto timeLoop = std::make_shared<CheckPointTimeLoop<Scalar>>(
+ getParam<Scalar>("TimeLoop.TStart", 0.0),
+ getParam<Scalar>("TimeLoop.InitialTimeStepSize"),
+ getParam<Scalar>("TimeLoop.TEnd")
+ );
```
+We set the maximum time step size allowed in the adaptive time stepping scheme.
-### Assembler, linear and nonlinear solver
+```cpp
+ timeLoop->setMaxTimeStepSize(getParam<Scalar>("TimeLoop.MaxTimeStepSize"));
+```
+Next, we choose the type of assembler, linear solver and PDE solver
+and construct instances of these classes.
```cpp
- // the assembler with time loop for a transient problem
using Assembler = FVAssembler<TypeTag, DiffMethod::numeric>;
- auto assembler = std::make_shared<Assembler>(problem, gridGeometry, gridVariables, timeLoop, oldSol);
+ using LinearSolver = SSORBiCGSTABIstlSolver<LinearSolverTraits<GridGeometry>, LinearAlgebraTraitsFromAssembler<Assembler>>;
+ using Solver = Dumux::NewtonSolver<Assembler, LinearSolver>;
- // the linear solver
- using LinearSolver = SSORBiCGSTABIstlSolver<
- LinearSolverTraits<GridGeometry>,
- LinearAlgebraTraitsFromAssembler<Assembler>
- >;
+ auto oldSol = sol;
+ auto assembler = std::make_shared<Assembler>(problem, gridGeometry, gridVariables, timeLoop, oldSol);
auto linearSolver = std::make_shared<LinearSolver>(gridGeometry->gridView(), gridGeometry->dofMapper());
-
- // the solver
- using Solver = Dumux::NewtonSolver<Assembler, LinearSolver>;
Solver solver(assembler, linearSolver);
```
-
-### Time loop
+The time loop is where most of the actual computations happen.
+We assemble and solve the linear system of equations, update the solution,
+write the solution to a VTK file and continue until the we reach the
+final simulation time.
```cpp
- // time loop
timeLoop->start(); do
{
- // assemble & solve
+ // Assemble & solve
+ // Passing the time loop enables the solver to repeat the time step
+ // with a reduced time step size if the Newton solver does not converge.
solver.solve(sol, *timeLoop);
// make the new solution the old solution
@@ -460,7 +376,7 @@ initialize the grid variables, also computing secondary variables.
// advance to the time loop to the next step
timeLoop->advanceTimeStep();
- // write VTK output
+ // write VTK output (concentration field and chemical potential)
vtkWriter.write(timeLoop->time());
// report statistics of this time step
@@ -470,15 +386,9 @@ initialize the grid variables, also computing secondary variables.
timeLoop->setTimeStepSize(solver.suggestTimeStepSize(timeLoop->timeStepSize()));
} while (!timeLoop->finished());
-```
-
-
-### Finalize
-
-
-```cpp
- timeLoop->finalize(leafGridView.comm());
+ // print the final report
+ timeLoop->finalize(gridGeometry->gridView().comm());
return 0;
}
```
@@ -487,6 +397,6 @@ initialize the grid variables, also computing secondary variables.
</details>
-| [:arrow_left: Back to the main documentation](../README.md) | [:arrow_left: Go back to part 1](modelheader.md) |
+| [:arrow_left: Back to the main documentation](../README.md) | [:arrow_left: Go back to part 1](model.md) |
|---|---:|
diff --git a/examples/cahn_hilliard/doc/modelheader.md b/examples/cahn_hilliard/doc/model.md
similarity index 50%
rename from examples/cahn_hilliard/doc/modelheader.md
rename to examples/cahn_hilliard/doc/model.md
index 125ea472c8..b9fc2a3d73 100644
--- a/examples/cahn_hilliard/doc/modelheader.md
+++ b/examples/cahn_hilliard/doc/model.md
@@ -1,30 +1,81 @@
<!-- Important: This file has been automatically generated by generate_example_docs.py. Do not edit this file directly! -->
-| [:arrow_left: Back to the main documentation](../README.md) | [:arrow_right: Continue with part 2](mainfile.md) |
+| [:arrow_left: Back to the main documentation](../README.md) | [:arrow_right: Continue with part 2](main.md) |
|---|---:|
-# Volume variables, local residual and model traits (`model.hh`)
-In this example the file `model.hh` contains the classes `CahnHilliardModelVolumeVariables` and
-`CahnHilliardModelLocalResidual` as well as general model traits and properties.
+# Cahn-Hilliard equation model definition
-## VolumeVariables
+In the file `model.hh`, we define the model equations and
+set all default model properties. The setup consist of four steps:
+1. Create a model type tag (used to specialize properties)
+2. Define the volume variables class computing and storing variables for a control volume
+3. Define the local residual class implementing the discrete equation
+4. Specialize important properties of the model such that Dumux knows how to assemble the system matrix
-The volume variables store the local element volume variables, both primary and secondary.
+We implement the classes
+
+* `CahnHilliardModel` (step 1),
+* `CahnHilliardModelVolumeVariables` (step 2) and
+* `CahnHilliardModelLocalResidual` (step 3).
+
+__Table of contents__
+
+[TOC]
+
+We start in `model.hh` with the necessary header includes:
+<details><summary> Click to show includes</summary>
+
+```cpp
+
+#include <dumux/common/math.hh>
+#include <dumux/common/properties.hh>
+#include <dumux/common/numeqvector.hh>
+```
+
+</details>
+
+## 1. Property Tag
+
+The property tag is simply an empty struct with the name `CahnHilliardModel`.
<details open>
<summary><b>Click to hide/show the file documentation</b> (or inspect the [source code](../model.hh))</summary>
-### The `CahnHilliardModelVolumeVariables` class
+```cpp
+namespace Dumux::Properties::TTag {
+struct CahnHilliardModel {};
+} // end namespace Dumux::Properties::TTag
+```
+
+
+</details>
+
+## 2. The volume variables
+
+The volume variables store the control volume variables, both primary and secondary.
+Let's have a look at the class implementation.
+
+
+<details open>
+<summary><b>Click to hide/show the file documentation</b> (or inspect the [source code](../model.hh))</summary>
+
+
+We write out a full volume variable class compatible with the Dumux assembly.
+We have to fulfill the same interface as the [`BasicVolumeVariables`](https://git.iws.uni-stuttgart.de/dumux-repositories/dumux/-/blob/master/dumux/common/volumevariables.hh).
+To shorten the code, we could also inherit from `BasicVolumeVariables`. However,
+we want to show the entire interface here as a basis for the implementation of more complex variables.
+Note that in `update` we could, for example, compute some secondary variable that is a
+possibly nonlinear function of the primary variables. That secondary variable could be exposed
+via an interface and then used below in the `CahnHilliardModelLocalResidual` class implementation.
```cpp
namespace Dumux {
-
template <class Traits>
class CahnHilliardModelVolumeVariables
{
@@ -38,16 +89,13 @@ public:
```
-### Update variables
-
-The `update` function stores the local primary variables of the current solution and
-potentially recomputes secondary variables.
+**Update variables:** The `update` function stores the local primary variables of the current solution and
+potentially computes secondary variables. Secondary variables can be nonlinear functions
+of the primary variables.
```cpp
- /*!
- * \brief Update all quantities for a given control volume
- */
+ //! Update all quantities for a given control volume
template<class ElementSolution, class Problem, class Element, class SubControlVolume>
void update(const ElementSolution& elemSol,
const Problem& problem,
@@ -59,10 +107,9 @@ potentially recomputes secondary variables.
```
-### Access functions
-
-Named and generic functions to access different primary variables
-
+**Access functions:** Named and generic functions to access different primary variables.
+The volumevariables are also expected to return an extrusion factor.
+Here we don't extrude our domain and therefore return $1.0$.
```cpp
Scalar concentration() const
@@ -76,113 +123,53 @@ Named and generic functions to access different primary variables
const PrimaryVariables& priVars() const
{ return priVars_; }
-```
-
-
-### Extrusion factor
-
-The volumevariables are also expected to return the extrusion factor
-
-```cpp
- // for compatibility with more general models
Scalar extrusionFactor() const
{ return 1.0; }
-```
-
-
-### Storage of local variables
-
-```cpp
private:
PrimaryVariables priVars_;
};
-
} // end namespace Dumux
```
</details>
-## LocalResidual
+## 3. The local residual
-The local residual defines the discretized and integrated partial differential equation through
-terms for storage, fluxes and sources, with the residual given as
-d/dt storage + div(fluxes) - sources = 0
-The individual terms can be further adjusted or replaced in the more specific problem.
+The local residual assembles the contribution to the residual for
+all degrees of freedom associated with an element.
+See [examples/diffusion/doc/model.hh](https://git.iws.uni-stuttgart.de/dumux-repositories/dumux/-/blob/master/examples/diffusion/doc/main.md)
+for a more detailed explanation for control-volume finite element method local residuals.
+Let's have a look at the class implementation.
<details open>
<summary><b>Click to hide/show the file documentation</b> (or inspect the [source code](../model.hh))</summary>
-### Include headers
-
-
-```cpp
-#include <dumux/common/math.hh>
-// use the property system
-#include <dumux/common/properties.hh>
-// common DuMux vector for discretized equations
-#include <dumux/common/numeqvector.hh>
-```
-
-
-### The local residual class `CahnHilliardModelLocalResidual` inherits from a base class set in
+The class `CahnHilliardModelLocalResidual` inherits from a base class set in
the model properties, depending on the discretization scheme.
-
+See [examples/diffusion/doc/model.hh](https://git.iws.uni-stuttgart.de/dumux-repositories/dumux/-/blob/master/examples/diffusion/doc/main.md)
+for details on the `BaseLocalResidual`.
```cpp
namespace Dumux {
-
template<class TypeTag>
class CahnHilliardModelLocalResidual
: public GetPropType<TypeTag, Properties::BaseLocalResidual>
{
- // the base local residual is selected depending on the chosen discretization scheme
- using ParentType = GetPropType<TypeTag, Properties::BaseLocalResidual>;
```
-<details><summary> Click to show alias definitions</summary>
-
-```cpp
- using Scalar = GetPropType<TypeTag, Properties::Scalar>;
- using Problem = GetPropType<TypeTag, Properties::Problem>;
-
- using NumEqVector = Dumux::NumEqVector<GetPropType<TypeTag, Properties::PrimaryVariables>>;
- using VolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::VolumeVariables;
-
- using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
- using ElementFluxVariablesCache = typename GetPropType<TypeTag, Properties::GridFluxVariablesCache>::LocalView;
- using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
-
- using SubControlVolume = typename FVElementGeometry::SubControlVolume;
- using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
- using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
- using Element = typename GridView::template Codim<0>::Entity;
-
- using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
- using Indices = typename ModelTraits::Indices;
- static constexpr int dimWorld = GridView::dimensionworld;
-public:
- using ParentType::ParentType;
-```
-
-</details>
-
-### The storage term
-The function `computeStorage` receives the volumevariables at the previous or current time
-step and computes the value of the storage terms.
+**Storage term:** The function `computeStorage` receives the volume variables
+at the previous or current time step and computes the value of the storage terms.
In this case the mass balance equation is a conservation equation of the concentration and
the equation for the chemical potential does not have a storage term.
```cpp
- /*!
- * \brief Evaluate the rate of change of all conserved quantities
- */
NumEqVector computeStorage(const Problem& problem,
const SubControlVolume& scv,
const VolumeVariables& volVars) const
@@ -194,20 +181,19 @@ the equation for the chemical potential does not have a storage term.
}
```
-
-### The flux terms
-`computeFlux` computes the fluxes over a subcontrolvolumeface, including the integration over
+**Flux term:** The function `computeFlux` computes the fluxes
+over a sub control volume faces, including the integration over
the area of the face.
+```math
+\begin{aligned}
+F_{K,\sigma,0} &= -M \sum_{B \in \mathcal{B}_K} \mu_{h,B} \nabla N_B \cdot\boldsymbol{n} \vert \sigma \vert \\
+F_{K,\sigma,1} &= -\gamma \sum_{B \in \mathcal{B}_K} c_{h,B} \nabla N_B \cdot\boldsymbol{n} \vert \sigma \vert
+\end{aligned}
+````
+
```cpp
- /*!
- * \brief Evaluate the fluxes over a face of a sub control volume
- * Here we evaluate the flow rate, F1 = -M∇mu·n A, F2 = -gamma∇c·n A
- *
- * TODO: why is this called flux, if we expect it to be integrated already?
- * computeFluxIntegral?
- */
NumEqVector computeFlux(const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
@@ -215,40 +201,46 @@ the area of the face.
const SubControlVolumeFace& scvf,
const ElementFluxVariablesCache& elemFluxVarsCache) const
{
+ static_assert(DiscretizationMethods::isCVFE<typename GridGeometry::DiscretizationMethod>,
+ "This local residual is hard-coded to control-volume finite element schemes");
+
const auto& fluxVarCache = elemFluxVarsCache[scvf];
Dune::FieldVector<Scalar, dimWorld> gradConcentration(0.0);
Dune::FieldVector<Scalar, dimWorld> gradChemicalPotential(0.0);
for (const auto& scv : scvs(fvGeometry))
{
const auto& volVars = elemVolVars[scv];
- gradConcentration.axpy(volVars.concentration(), fluxVarCache.gradN(scv.indexInElement()));
- gradChemicalPotential.axpy(volVars.chemicalPotential(), fluxVarCache.gradN(scv.indexInElement()));
+ // v.axpy(a, w) means v <- v + a*w
+ gradConcentration.axpy(
+ volVars.concentration(),
+ fluxVarCache.gradN(scv.indexInElement())
+ );
+ gradChemicalPotential.axpy(
+ volVars.chemicalPotential(),
+ fluxVarCache.gradN(scv.indexInElement())
+ );
}
- const auto M = problem.mobility();
- const auto gamma = problem.surfaceTension();
-
NumEqVector flux;
- flux[Indices::massBalanceEqIdx] = -1.0*vtmv(scvf.unitOuterNormal(), M, gradChemicalPotential)*scvf.area();
- flux[Indices::chemicalPotentialEqIdx] = -1.0*vtmv(scvf.unitOuterNormal(), gamma, gradConcentration)*scvf.area();
+ // Compute the flux with `vtmv` (vector transposed times matrix times vector).
+ // The mobility and surface tension coefficients comes from the `problem` (see Part 2 of the example).
+ flux[Indices::massBalanceEqIdx] = -1.0*vtmv(
+ scvf.unitOuterNormal(), problem.mobility(), gradChemicalPotential
+ )*scvf.area();
+ flux[Indices::chemicalPotentialEqIdx] = -1.0*vtmv(
+ scvf.unitOuterNormal(), problem.surfaceTension(), gradConcentration
+ )*scvf.area();
return flux;
}
```
-
-### Source terms
-
-`computeSource` defines the sources terms at a sub control volume.
+**Source term:** The function `computeSource` computes the sources terms for a sub control volume.
We implement a model-specific source term for the chemical potential equation before
deferring further implementation to the problem where we add the derivative of the free
-energy. The default implementation of this function also defers the calculation to the
-problem.
+energy.
```cpp
- /*!
- * \brief Calculate the source term of the equation
- */
NumEqVector computeSource(const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
@@ -256,90 +248,59 @@ problem.
const SubControlVolume &scv) const
{
NumEqVector source(0.0);
-
+ source[Indices::massBalanceEqIdx] = 0.0;
source[Indices::chemicalPotentialEqIdx] = elemVolVars[scv].chemicalPotential();
-
// add contributions from problem (e.g. double well potential)
source += problem.source(element, fvGeometry, elemVolVars, scv);
-
return source;
}
};
-
} // end namespace Dumux
```
</details>
-## Model properties/traits
+## 4. The model properties/traits
-We set some general model traits and properties, using a TypeTag `CahnHilliardModel`.
-For this type tag we define a `ModelTraits` struct and set a number of properties by specifying
-further structs within the `Dumux::Properties` namespace.
+In the `Dumux::Properties` namespace, we specialize properties for
+the created type tag `CahnHilliardModel`.
<details open>
<summary><b>Click to hide/show the file documentation</b> (or inspect the [source code](../model.hh))</summary>
-### Include the header for the property system and common properties and switch to the
-`Properties` namespace.
-
-
```cpp
-#include <dumux/common/properties.hh>
-
namespace Dumux::Properties {
```
-
-### Define the type tag
-
+The type of the local residual is the class defined above.
```cpp
-namespace TTag {
-struct CahnHilliardModel {};
-} // end namespace TTag
+template<class TypeTag>
+struct LocalResidual<TypeTag, TTag::CahnHilliardModel>
+{ using type = CahnHilliardModelLocalResidual<TypeTag>; };
```
-
-### Basic model properties
-
-Define some general properties to be used for this modedl such as a datatype for scalars and the
-default vector for the primary variables. Here we can also use the model traits we specify below.
-
+The default scalar type is double.
+We compute with double precision floating point numbers.
```cpp
-//! Set the default type of scalar values to double
template<class TypeTag>
-struct Scalar<TypeTag, TTag:: CahnHilliardModel >
+struct Scalar<TypeTag, TTag::CahnHilliardModel>
{ using type = double; };
-
-//! Set the default primary variable vector to a vector of size of number of equations
-template<class TypeTag>
-struct PrimaryVariables<TypeTag, TTag:: CahnHilliardModel >
-{
- using type = Dune::FieldVector<
- GetPropType<TypeTag, Properties::Scalar>,
- GetPropType<TypeTag, Properties::ModelTraits>::numEq()
- >;
-};
```
-
-### Model traits
-
-We specify general traits of the implemented model, defining indices (often in `indices.hh`)
-and the number of equations in the model.
-
+The model traits specify some information about our equation system.
+Here we have two equations. The indices allow to access primary variables
+and equations with a named indices.
```cpp
-//! Set the model traits property
template<class TypeTag>
struct ModelTraits<TypeTag, TTag::CahnHilliardModel>
{
- struct Traits
+ struct type
{
struct Indices
{
@@ -352,24 +313,26 @@ struct ModelTraits<TypeTag, TTag::CahnHilliardModel>
static constexpr int numEq() { return 2; }
};
-
- using type = Traits;
};
```
-
-### Further model properties
-
-Define further properties of the model, selecting the local residual and volumevariables defined
-above.
-
+The primary variable vector has entries of type `Scalar` and is
+as large as the number of equations (here 2) but we keep it general.
```cpp
template<class TypeTag>
-struct LocalResidual<TypeTag, TTag::CahnHilliardModel>
-{ using type = CahnHilliardModelLocalResidual<TypeTag>; };
+struct PrimaryVariables<TypeTag, TTag::CahnHilliardModel>
+{
+ using type = Dune::FieldVector<
+ GetPropType<TypeTag, Properties::Scalar>,
+ GetPropType<TypeTag, Properties::ModelTraits>::numEq()
+ >;
+};
+```
+
+Finally, the type of the volume variables is the class defined above.
-//! Set the volume variables property
+```cpp
template<class TypeTag>
struct VolumeVariables<TypeTag, TTag::CahnHilliardModel>
{
@@ -389,6 +352,6 @@ struct VolumeVariables<TypeTag, TTag::CahnHilliardModel>
</details>
-| [:arrow_left: Back to the main documentation](../README.md) | [:arrow_right: Continue with part 2](mainfile.md) |
+| [:arrow_left: Back to the main documentation](../README.md) | [:arrow_right: Continue with part 2](main.md) |
|---|---:|
diff --git a/examples/cahn_hilliard/main.cc b/examples/cahn_hilliard/main.cc
index 4e18dcff99..b91a6d9acb 100644
--- a/examples/cahn_hilliard/main.cc
+++ b/examples/cahn_hilliard/main.cc
@@ -17,43 +17,59 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
-// # Problem, test properties/traits and main program flow (`main.cc`)
-// In this example the file `main.cc` contains the problem class `CahnHilliardTestProblem`,
-// properties and traits specific to the test case as well as the main program flow in the form of
-// `main` function.
+// # Cahn-Hilliard equation main program
//
-#include <config.h>
-// ## Problem
-// The __problem class__ defines boundary conditions and extends the storage term defined in the
-// model's localresidual by the derivative of the free energy.
+// The file [`main.cc`](main.cc) contains the problem class `CahnHilliardTestProblem`,
+// properties and traits specific to the test case as well as the `main` function.
+// The setup consists of four steps:
+// 1. Define the problem setting boundary conditions and the diffusion coefficient
+// 2. Configure the property system reusing the model defined in Part 1
+// 3. Define a function for setting the random initial condition
+// 4. The main program defining all steps of the program
//
-// [[content]]
-// ### Include headers
+// __Table of contents__
//
-// [[codeblock]]
-// use the property system and runtime parameters
+// [TOC]
+//
+// We start in `main.cc` with the necessary header includes:
+// [[details]] includes
+#include <config.h>
+#include <type_traits>
+#include <random>
+
+#include <dumux/common/initialize.hh>
#include <dumux/common/properties.hh>
#include <dumux/common/parameters.hh>
-// common DuMux vector for discretized equations
#include <dumux/common/numeqvector.hh>
-// types of boundary conditions to use
#include <dumux/common/boundarytypes.hh>
-// generic problem for finite volume simulations
#include <dumux/common/fvproblem.hh>
-// [[/codeblock]]
+
+#include <dumux/discretization/box.hh>
+#include <dumux/linear/linearsolvertraits.hh>
+#include <dumux/linear/linearalgebratraits.hh>
+#include <dumux/linear/istlsolvers.hh>
+#include <dumux/nonlinear/newtonsolver.hh>
+#include <dumux/assembly/fvassembler.hh>
+
+#include <dune/grid/yaspgrid.hh>
+#include <dumux/io/vtkoutputmodule.hh>
+#include <dumux/io/grid/gridmanager_yasp.hh>
+
+#include "model.hh"
+// [[/details]]
+
//
-// ### The problem class `CahnHilliardTestProblem`
-// In this abstract problem we inherit from the generic `FVProblem`.
+// ## 1. The problem class
//
-// [[codeblock]]
+// The problem class implements boundary conditions and the source term.
+// It also provides an interface that is used by the local residual (see Part 1) to obtain
+// the mobility and surface tension. The values are read from the parameter configuration tree.
+// [[content]]
namespace Dumux {
-
template<class TypeTag>
class CahnHilliardTestProblem : public FVProblem<TypeTag>
{
-// [[/codeblock]]
-// [[details]] alias definitions and local variables
-// [[codeblock]]
+ // [[details]] boilerplate code
using ParentType = FVProblem<TypeTag>;
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
using FVElementGeometry = typename GridGeometry::LocalView;
@@ -67,9 +83,9 @@ class CahnHilliardTestProblem : public FVProblem<TypeTag>
using NumEqVector = Dumux::NumEqVector<PrimaryVariables>;
using BoundaryTypes = Dumux::BoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>;
using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
-// [[/codeblock]]
-// [[/details]]
-// [[codeblock]]
+ // [[/details]]
+ // In the constructor, we read the parameter constants from the
+ // parameter tree (which is initialized with the content of `params.input`).
public:
CahnHilliardTestProblem(std::shared_ptr<const GridGeometry> gridGeometry)
: ParentType(gridGeometry)
@@ -78,15 +94,10 @@ public:
surfaceTension_ = getParam<Scalar>("Problem.SurfaceTension");
energyScale_ = getParam<Scalar>("Problem.EnergyScale");
}
- // [[/codeblock]]
- //
- // ### Problem source term
- //
- // Here we implement the derivative of the free energy, setting a source for the equation for
- // the chemical potential. The `computeSource` function in the local residual adds the terms
- // defined here.
- //
- // [[codeblock]]
+
+ // In the `source` function, we implement the derivative of the free energy.
+ // This demonstrates how parts of the local residual can be split into model specific
+ // parts and parts that might change from scenario to scenario.
template<class ElementVolumeVariables>
NumEqVector source(const Element &element,
const FVElementGeometry& fvGeometry,
@@ -98,13 +109,9 @@ public:
values[Indices::chemicalPotentialEqIdx] = -energyScale_*2.0*c*(2.0*c*c - 3*c + 1);
return values;
}
- // [[/codeblock]]
- //
- // ### Boundary conditions
- //
+
// For the boundary we choose boundary flux (or Neumann) conditions for all equations and on
// every part of the boundary, specifying zero flux everywhere for both equations.
- //
// [[codeblock]]
BoundaryTypes boundaryTypesAtPos(const GlobalPosition& globalPos) const
{
@@ -116,11 +123,10 @@ public:
NumEqVector neumannAtPos(const GlobalPosition& globalPos) const
{ return { 0.0, 0.0 }; }
// [[/codeblock]]
- //
- // [[details]] coefficients and private variables
- // The problem class offers access to the mobility and surface tension coefficients as read from
- // the parameter file by default `params.input`.
- //
+
+ // The parameters interfaces are used in the local residual (see Part 1).
+ // We can name this interface however we want as long as we adapt the calling site
+ // in the `LocalResidual` class in `model.hh`.
// [[codeblock]]
Scalar mobility() const
{ return mobility_; }
@@ -133,118 +139,48 @@ private:
Scalar surfaceTension_;
Scalar energyScale_;
};
-
} // end namespace Dumux
// [[/codeblock]]
-// [[/details]]
// [[/content]]
-//
-// ## Test case properties/traits
-// Within the `Dumux::Properties` namespace we specialize properties and traits to the considered
-// test case by using the test's TypeTag.
-//
-// [[content]]
-//
-// ### Include headers
-//
-// [[codeblock]]
-// Include the grid to be used
-#include <dune/grid/yaspgrid.hh>
-// The header for the box discretization scheme
-#include <dumux/discretization/box.hh>
-// The model header including the model traits and properties
-#include "model.hh"
-// [[/codeblock]]
-//
-// ### TypeTag `CahnHilliardTest`
-//
-// We define a type tag for the test case, allowing us to further specify properties and traits. To
-// use those set for the Cahn-Hilliard model we inherit from its type tag.
-//
-// [[codeblock]]
-namespace Dumux::Properties {
-// Inheriting properties of the Cahn-Hilliard model and the box finite volume discretization
-namespace TTag {
-struct CahnHilliardTest { using InheritsFrom = std::tuple<CahnHilliardModel, BoxModel>; };
-} // end namespace TTag
-// [[/codeblock]]
//
-// ### Test properties
+// ## 2. Property tag and specializations
//
-// We specify a grid to be used in the test, select our problem class and enable caching.
+// We create a new tag `DiffusionTest` that inherits all properties
+// specialized for the tag `DiffusionModel` (we created this in Part 1)
+// and the tag `BoxModel` which provides types relevant for the spatial
+// discretization scheme (see [dumux/discretization/box.hh](https://git.iws.uni-stuttgart.de/dumux-repositories/dumux/-/blob/master/dumux/discretization/box.hh)).
//
-// [[codeblock]]
-// Set the grid type
-template<class TypeTag>
-struct Grid<TypeTag, TTag::CahnHilliardTest>
-{ using type = Dune::YaspGrid<2>; };
-
-// Select the problem class defined above
-template<class TypeTag>
-struct Problem<TypeTag, TTag::CahnHilliardTest>
-{ using type = CahnHilliardTestProblem<TypeTag>; };
-
-// Enable caching
-template<class TypeTag>
-struct EnableGridVolumeVariablesCache<TypeTag, TTag::CahnHilliardTest>
-{ static constexpr bool value = true; };
+// [[content]]
+namespace Dumux::Properties::TTag {
+struct CahnHilliardTest
+{
+ using InheritsFrom = std::tuple<CahnHilliardModel, BoxModel>;
-template<class TypeTag>
-struct EnableGridFluxVariablesCache<TypeTag, TTag::CahnHilliardTest>
-{ static constexpr bool value = true; };
+ using Grid = Dune::YaspGrid<2>;
-template<class TypeTag>
-struct EnableGridGeometryCache<TypeTag, TTag::CahnHilliardTest>
-{ static constexpr bool value = true; };
+ template<class TypeTag>
+ using Problem = CahnHilliardTestProblem<TypeTag>;
+ using EnableGridGeometryCache = std::true_type;
+ using EnableGridVolumeVariablesCache = std::true_type;
+ using EnableGridFluxVariablesCache = std::true_type;
+};
} // end namespace Dumux::Properties
-// [[/codeblock]]
// [[/content]]
//
-// ## The main program flow
-// The main program flow in the `main` function is often the only content of `main.cc`. It sets up
-// the simulation framework, initializes runtime parameters, creates the grid and storage vectors
-// for the variables, primary and secondary. It specifies and constructs and assembler, which
-// assembles the discretized residual and system matrix (Jacobian of the model residual), as well as
-// linear and nonlinear solvers that solve the resulting linear system and handle the convergence of
-// nonlinear iterations. The time loop controls the time stepping, with adaptive time step size in
-// coordination with the nonlinear solver.
-//
-// [[content]]
-//
-// ### Include headers
-//
-// [[codeblock]]
-// standard header to generate random initial data
-#include <random>
-// common DuMux header for parallelization
-#include <dumux/common/initialize.hh>
-// headers to use the property system and runtime parameters
-#include <dumux/common/properties.hh>
-#include <dumux/common/parameters.hh>
-// module for VTK output, to write out fields of interest
-#include <dumux/io/vtkoutputmodule.hh>
-// gridmanager for the grid used in the test
-#include <dumux/io/grid/gridmanager_yasp.hh>
-// headers for linear and non-linear solvers as well as the assembler
-#include <dumux/linear/linearsolvertraits.hh>
-#include <dumux/linear/linearalgebratraits.hh>
-#include <dumux/linear/istlsolvers.hh>
-#include <dumux/nonlinear/newtonsolver.hh>
-#include <dumux/assembly/fvassembler.hh>
-// [[/codeblock]]
-//
-// ### Creating the initial solution
+// ## 3. Creating the initial solution
//
-// We define a helper struct and function to handle communication for parallel runs.
-// For our initial conditions we create a random field of values around a mean of 0.42.
-// The random values are created with an offset based on the processor rank for communication
-// purposes, which is removed afterwards. For more information see the description of the
-// diffusion example
-// [examples/diffusion/doc/main.md](https://git.iws.uni-stuttgart.de/dumux-repositories/dumux/-/blob/master/examples/diffusion/doc/main.md)
+// To initialize with a random field in parallel, where each processor
+// creates its own random number sequence, we need to communicate the
+// resulting values on the processor border and overlap.
+// See [Diffusion equation example](https://git.iws.uni-stuttgart.de/dumux-repositories/dumux/-/blob/master/examples/diffusion/docs/main.md).
+// for details. Here in addition, we need to provide a custom scatter operation
+// that handles vector types. We only need to communicate the first entry (concentration).
//
+// [[content]]
// [[codeblock]]
+// functor for data communication with MPI
struct MinScatter
{
template<class A, class B>
@@ -252,6 +188,7 @@ struct MinScatter
{ a[0] = std::min(a[0], b[0]); }
};
+// create the random initial solution
template<class SolutionVector, class GridGeometry>
SolutionVector createInitialSolution(const GridGeometry& gg)
{
@@ -267,7 +204,8 @@ SolutionVector createInitialSolution(const GridGeometry& gg)
sol[n][1] = 0.0;
}
- // Take the value of the processor with the minimum rank and subtract the rank offset
+ // We take the value of the processor with the minimum rank
+ // and subtract the rank offset
if (gg.gridView().comm().size() > 1)
{
Dumux::VectorCommDataHandle<
@@ -285,127 +223,91 @@ SolutionVector createInitialSolution(const GridGeometry& gg)
return sol;
}
// [[/codeblock]]
+// [[/content]]
//
-// ### The main function
+// ## 4. The main program
//
-// The main function takes the command line arguments, optionally specifying an input file of
-// parameters and/or individual key-value pairs of runtime parameters.
+// The `main` function sets up the simulation framework, initializes runtime parameters,
+// creates the grid and storage vectors
+// for the variables, primary and secondary. It specifies and constructs and assembler, which
+// assembles the discretized residual and system matrix (Jacobian of the model residual), as well as
+// a linear solver. A Newton method is used to solve the nonlinear equations where in each Newton iteration
+// the Jacobian and the residual needs to be reassembled and the resulting linear system is solved.
+// The time loop controls the time stepping. Here, we let the Newton solver suggest an adaption of
+// the time step size based on the convergence history of the nonlinear solver.
//
-// [[codeblock]]
+// [[content]]
int main(int argc, char** argv)
{
using namespace Dumux;
- // define the type tag for this problem
- using TypeTag = Properties::TTag::CahnHilliardTest;
- // [[/codeblock]]
- //
- // We initialize parallelization backends as well as runtime parameters
- //
- // [[codeblock]]
- // maybe initialize MPI and/or multithreading backend
+ // We initialize MPI and/or multithreading backend. When not running
+ // in parallel the MPI setup is skipped.
Dumux::initialize(argc, argv);
- // initialize parameter tree
+ // Then we initialize parameter tree.
Parameters::init(argc, argv);
- // [[/codeblock]]
- //
- // ### Grid setup
- //
- // Set up the grid as well as a grid geometry to access the (sub-)control-volumes and their
- // faces
- //
- // [[codeblock]]
- // initialize the grid
- GridManager<GetPropType<TypeTag, Properties::Grid>> gridManager;
- gridManager.init();
- // we compute on the leaf grid view
- const auto& leafGridView = gridManager.grid().leafGridView();
+ // We create an alias for the type tag for this problem
+ // and extract type information through properties.
+ using TypeTag = Properties::TTag::CahnHilliardTest;
- // create the finite volume grid geometry
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using Grid = GetPropType<TypeTag, Properties::Grid>;
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
- auto gridGeometry = std::make_shared<GridGeometry>(leafGridView);
- // [[/codeblock]]
- //
- // ### Problem setup
- //
- // We instantiate also the problem class according to the test properties
- //
- // [[codeblock]]
- // the problem (initial and boundary conditions)
using Problem = GetPropType<TypeTag, Properties::Problem>;
- auto problem = std::make_shared<Problem>(gridGeometry);
- // [[/codeblock]]
- //
- // ### Applying initial conditions
- //
- // After writing the initial data to the storage for previous and current time-step, we
- // initialize the grid variables, also computing secondary variables.
- //
- // [[codeblock]]
- // the solution vector
using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
- auto sol = createInitialSolution<SolutionVector>(*gridGeometry);
- // copy the vector to store state of previous time step
- auto oldSol = sol;
-
- // the grid variables
using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
+
+ // We initialize the grid, grid geometry, problem, solution vector, and grid variables.
+ GridManager<Grid> gridManager;
+ gridManager.init();
+
+ auto gridGeometry = std::make_shared<GridGeometry>(gridManager.grid().leafGridView());
+ auto problem = std::make_shared<Problem>(gridGeometry);
+ auto sol = createInitialSolution<SolutionVector>(*gridGeometry);
auto gridVariables = std::make_shared<GridVariables>(problem, gridGeometry);
gridVariables->init(sol);
- // [[/codeblock]]
- //
- // ### Initialize VTK output
- //
- // [[codeblock]]
- // initialize the vtk output module
+
+ // We initialize the VTK output module and write out the initial concentration and chemical potential
VtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables, sol, problem->name());
vtkWriter.addVolumeVariable([](const auto& vv){ return vv.concentration(); }, "c");
vtkWriter.addVolumeVariable([](const auto& vv){ return vv.chemicalPotential(); }, "mu");
vtkWriter.write(0.0);
- // [[/codeblock]]
- //
- // ### Set up time loop
- //
- // [[codeblock]]
- // get some time loop parameters
- using Scalar = GetPropType<TypeTag, Properties::Scalar>;
- const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
- const auto dt = getParam<Scalar>("TimeLoop.InitialTimeStepSize");
- const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
- // instantiate time loop
- auto timeLoop = std::make_shared<CheckPointTimeLoop<Scalar>>(0.0, dt, tEnd);
- timeLoop->setMaxTimeStepSize(maxDt);
- // [[/codeblock]]
- //
- // ### Assembler, linear and nonlinear solver
- //
- // [[codeblock]]
- // the assembler with time loop for a transient problem
- using Assembler = FVAssembler<TypeTag, DiffMethod::numeric>;
- auto assembler = std::make_shared<Assembler>(problem, gridGeometry, gridVariables, timeLoop, oldSol);
+ // We instantiate time loop using start and end time as well as
+ // the time step size from the parameter tree (`params.input`)
+ auto timeLoop = std::make_shared<CheckPointTimeLoop<Scalar>>(
+ getParam<Scalar>("TimeLoop.TStart", 0.0),
+ getParam<Scalar>("TimeLoop.InitialTimeStepSize"),
+ getParam<Scalar>("TimeLoop.TEnd")
+ );
- // the linear solver
- using LinearSolver = SSORBiCGSTABIstlSolver<
- LinearSolverTraits<GridGeometry>,
- LinearAlgebraTraitsFromAssembler<Assembler>
- >;
- auto linearSolver = std::make_shared<LinearSolver>(gridGeometry->gridView(), gridGeometry->dofMapper());
+ // We set the maximum time step size allowed in the adaptive time stepping scheme.
+ timeLoop->setMaxTimeStepSize(getParam<Scalar>("TimeLoop.MaxTimeStepSize"));
- // the solver
+ // Next, we choose the type of assembler, linear solver and PDE solver
+ // and construct instances of these classes.
+ using Assembler = FVAssembler<TypeTag, DiffMethod::numeric>;
+ using LinearSolver = SSORBiCGSTABIstlSolver<LinearSolverTraits<GridGeometry>, LinearAlgebraTraitsFromAssembler<Assembler>>;
using Solver = Dumux::NewtonSolver<Assembler, LinearSolver>;
+
+ auto oldSol = sol;
+ auto assembler = std::make_shared<Assembler>(problem, gridGeometry, gridVariables, timeLoop, oldSol);
+ auto linearSolver = std::make_shared<LinearSolver>(gridGeometry->gridView(), gridGeometry->dofMapper());
Solver solver(assembler, linearSolver);
- // [[/codeblock]]
- //
- // ### Time loop
+
+ // The time loop is where most of the actual computations happen.
+ // We assemble and solve the linear system of equations, update the solution,
+ // write the solution to a VTK file and continue until the we reach the
+ // final simulation time.
//
// [[codeblock]]
- // time loop
timeLoop->start(); do
{
- // assemble & solve
+ // Assemble & solve
+ // Passing the time loop enables the solver to repeat the time step
+ // with a reduced time step size if the Newton solver does not converge.
solver.solve(sol, *timeLoop);
// make the new solution the old solution
@@ -415,7 +317,7 @@ int main(int argc, char** argv)
// advance to the time loop to the next step
timeLoop->advanceTimeStep();
- // write VTK output
+ // write VTK output (concentration field and chemical potential)
vtkWriter.write(timeLoop->time());
// report statistics of this time step
@@ -425,13 +327,9 @@ int main(int argc, char** argv)
timeLoop->setTimeStepSize(solver.suggestTimeStepSize(timeLoop->timeStepSize()));
} while (!timeLoop->finished());
- // [[/codeblock]]
- //
- // ### Finalize
- //
- // [[codeblock]]
- timeLoop->finalize(leafGridView.comm());
+ // print the final report
+ timeLoop->finalize(gridGeometry->gridView().comm());
return 0;
}
// [[/codeblock]]
diff --git a/examples/cahn_hilliard/model.hh b/examples/cahn_hilliard/model.hh
index 9f69f9dada..f77c3d4f3f 100644
--- a/examples/cahn_hilliard/model.hh
+++ b/examples/cahn_hilliard/model.hh
@@ -20,21 +20,61 @@
#ifndef DUMUX_EXAMPLES_CAHN_HILLIARD_MODEL_HH
#define DUMUX_EXAMPLES_CAHN_HILLIARD_MODEL_HH
-// # Volume variables, local residual and model traits (`model.hh`)
-// In this example the file `model.hh` contains the classes `CahnHilliardModelVolumeVariables` and
-// `CahnHilliardModelLocalResidual` as well as general model traits and properties.
+// # Cahn-Hilliard equation model definition
//
-// ## VolumeVariables
+// In the file `model.hh`, we define the model equations and
+// set all default model properties. The setup consist of four steps:
+// 1. Create a model type tag (used to specialize properties)
+// 2. Define the volume variables class computing and storing variables for a control volume
+// 3. Define the local residual class implementing the discrete equation
+// 4. Specialize important properties of the model such that Dumux knows how to assemble the system matrix
//
-// The volume variables store the local element volume variables, both primary and secondary.
+// We implement the classes
+//
+// * `CahnHilliardModel` (step 1),
+// * `CahnHilliardModelVolumeVariables` (step 2) and
+// * `CahnHilliardModelLocalResidual` (step 3).
+//
+// __Table of contents__
+//
+// [TOC]
+//
+// We start in `model.hh` with the necessary header includes:
+// [[details]] includes
+#include <dumux/common/math.hh>
+#include <dumux/common/properties.hh>
+#include <dumux/common/numeqvector.hh>
+// [[/details]]
+//
+// ## 1. Property Tag
+//
+// The property tag is simply an empty struct with the name `CahnHilliardModel`.
+//
+// [[content]]
+// [[codeblock]]
+namespace Dumux::Properties::TTag {
+struct CahnHilliardModel {};
+} // end namespace Dumux::Properties::TTag
+// [[/codeblock]]
+// [[/content]]
+//
+// ## 2. The volume variables
+//
+// The volume variables store the control volume variables, both primary and secondary.
+// Let's have a look at the class implementation.
//
// [[content]]
//
-// ### The `CahnHilliardModelVolumeVariables` class
+// We write out a full volume variable class compatible with the Dumux assembly.
+// We have to fulfill the same interface as the [`BasicVolumeVariables`](https://git.iws.uni-stuttgart.de/dumux-repositories/dumux/-/blob/master/dumux/common/volumevariables.hh).
+// To shorten the code, we could also inherit from `BasicVolumeVariables`. However,
+// we want to show the entire interface here as a basis for the implementation of more complex variables.
+// Note that in `update` we could, for example, compute some secondary variable that is a
+// possibly nonlinear function of the primary variables. That secondary variable could be exposed
+// via an interface and then used below in the `CahnHilliardModelLocalResidual` class implementation.
//
// [[codeblock]]
namespace Dumux {
-
template <class Traits>
class CahnHilliardModelVolumeVariables
{
@@ -45,17 +85,14 @@ public:
using PrimaryVariables = typename Traits::PrimaryVariables;
//! export the indices type
using Indices = typename Traits::ModelTraits::Indices;
- // [[/codeblock]]
- //
- // ### Update variables
+// [[/codeblock]]
//
- // The `update` function stores the local primary variables of the current solution and
- // potentially recomputes secondary variables.
+ // **Update variables:** The `update` function stores the local primary variables of the current solution and
+ // potentially computes secondary variables. Secondary variables can be nonlinear functions
+ // of the primary variables.
//
// [[codeblock]]
- /*!
- * \brief Update all quantities for a given control volume
- */
+ //! Update all quantities for a given control volume
template<class ElementSolution, class Problem, class Element, class SubControlVolume>
void update(const ElementSolution& elemSol,
const Problem& problem,
@@ -66,10 +103,9 @@ public:
}
// [[/codeblock]]
//
- // ### Access functions
- //
- // Named and generic functions to access different primary variables
- //
+ // **Access functions:** Named and generic functions to access different primary variables.
+ // The volumevariables are also expected to return an extrusion factor.
+ // Here we don't extrude our domain and therefore return $1.0$.
// [[codeblock]]
Scalar concentration() const
{ return priVars_[Indices::concentrationIdx]; }
@@ -82,77 +118,53 @@ public:
const PrimaryVariables& priVars() const
{ return priVars_; }
- // [[/codeblock]]
- //
- // ### Extrusion factor
- //
- // The volumevariables are also expected to return the extrusion factor
- //
- // [[codeblock]]
- // for compatibility with more general models
+
Scalar extrusionFactor() const
{ return 1.0; }
- // [[/codeblock]]
- //
- // ### Storage of local variables
- //
- // [[codeblock]]
+
private:
PrimaryVariables priVars_;
};
-
} // end namespace Dumux
// [[/codeblock]]
// [[/content]]
//
-// ## LocalResidual
+// ## 3. The local residual
//
-// The local residual defines the discretized and integrated partial differential equation through
-// terms for storage, fluxes and sources, with the residual given as
-// d/dt storage + div(fluxes) - sources = 0
-// The individual terms can be further adjusted or replaced in the more specific problem.
+// The local residual assembles the contribution to the residual for
+// all degrees of freedom associated with an element.
+// See [examples/diffusion/doc/model.hh](https://git.iws.uni-stuttgart.de/dumux-repositories/dumux/-/blob/master/examples/diffusion/doc/main.md)
+// for a more detailed explanation for control-volume finite element method local residuals.
+// Let's have a look at the class implementation.
//
// [[content]]
//
-// ### Include headers
-//
-// [[codeblock]]
-#include <dumux/common/math.hh>
-// use the property system
-#include <dumux/common/properties.hh>
-// common DuMux vector for discretized equations
-#include <dumux/common/numeqvector.hh>
-// [[/codeblock]]
-//
-// ### The local residual class `CahnHilliardModelLocalResidual` inherits from a base class set in
+// The class `CahnHilliardModelLocalResidual` inherits from a base class set in
// the model properties, depending on the discretization scheme.
-//
-// [[codeblock]]
+// See [examples/diffusion/doc/model.hh](https://git.iws.uni-stuttgart.de/dumux-repositories/dumux/-/blob/master/examples/diffusion/doc/main.md)
+// for details on the `BaseLocalResidual`.
namespace Dumux {
-
template<class TypeTag>
class CahnHilliardModelLocalResidual
: public GetPropType<TypeTag, Properties::BaseLocalResidual>
{
+ // [[exclude]]
// the base local residual is selected depending on the chosen discretization scheme
using ParentType = GetPropType<TypeTag, Properties::BaseLocalResidual>;
- // [[/codeblock]]
- //
- // [[details]] alias definitions
- // [[codeblock]]
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using Problem = GetPropType<TypeTag, Properties::Problem>;
-
using NumEqVector = Dumux::NumEqVector<GetPropType<TypeTag, Properties::PrimaryVariables>>;
- using VolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::VolumeVariables;
- using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
- using ElementFluxVariablesCache = typename GetPropType<TypeTag, Properties::GridFluxVariablesCache>::LocalView;
- using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
+ using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
+ using VolumeVariables = typename GridVariables::GridVolumeVariables::VolumeVariables;
+ using ElementVolumeVariables = typename GridVariables::GridVolumeVariables::LocalView;
+ using ElementFluxVariablesCache = typename GridVariables::GridFluxVariablesCache::LocalView;
- using SubControlVolume = typename FVElementGeometry::SubControlVolume;
- using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
- using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
+ using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
+ using FVElementGeometry = typename GridGeometry::LocalView;
+ using SubControlVolume = typename GridGeometry::SubControlVolume;
+ using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
+ using GridView = typename GridGeometry::GridView;
using Element = typename GridView::template Codim<0>::Entity;
using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
@@ -160,19 +172,14 @@ class CahnHilliardModelLocalResidual
static constexpr int dimWorld = GridView::dimensionworld;
public:
using ParentType::ParentType;
- // [[/codeblock]]
- // [[/details]]
+ // [[/exclude]]
//
- // ### The storage term
- // The function `computeStorage` receives the volumevariables at the previous or current time
- // step and computes the value of the storage terms.
+ // **Storage term:** The function `computeStorage` receives the volume variables
+ // at the previous or current time step and computes the value of the storage terms.
// In this case the mass balance equation is a conservation equation of the concentration and
// the equation for the chemical potential does not have a storage term.
//
// [[codeblock]]
- /*!
- * \brief Evaluate the rate of change of all conserved quantities
- */
NumEqVector computeStorage(const Problem& problem,
const SubControlVolume& scv,
const VolumeVariables& volVars) const
@@ -183,19 +190,19 @@ public:
return storage;
}
// [[/codeblock]]
- //
- // ### The flux terms
- // `computeFlux` computes the fluxes over a subcontrolvolumeface, including the integration over
+
+ // **Flux term:** The function `computeFlux` computes the fluxes
+ // over a sub control volume faces, including the integration over
// the area of the face.
//
+ // ```math
+ // \begin{aligned}
+ // F_{K,\sigma,0} &= -M \sum_{B \in \mathcal{B}_K} \mu_{h,B} \nabla N_B \cdot\boldsymbol{n} \vert \sigma \vert \\
+ // F_{K,\sigma,1} &= -\gamma \sum_{B \in \mathcal{B}_K} c_{h,B} \nabla N_B \cdot\boldsymbol{n} \vert \sigma \vert
+ // \end{aligned}
+ // ````
+ //
// [[codeblock]]
- /*!
- * \brief Evaluate the fluxes over a face of a sub control volume
- * Here we evaluate the flow rate, F1 = -M∇mu·n A, F2 = -gamma∇c·n A
- *
- * TODO: why is this called flux, if we expect it to be integrated already?
- * computeFluxIntegral?
- */
NumEqVector computeFlux(const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
@@ -203,38 +210,45 @@ public:
const SubControlVolumeFace& scvf,
const ElementFluxVariablesCache& elemFluxVarsCache) const
{
+ static_assert(DiscretizationMethods::isCVFE<typename GridGeometry::DiscretizationMethod>,
+ "This local residual is hard-coded to control-volume finite element schemes");
+
const auto& fluxVarCache = elemFluxVarsCache[scvf];
Dune::FieldVector<Scalar, dimWorld> gradConcentration(0.0);
Dune::FieldVector<Scalar, dimWorld> gradChemicalPotential(0.0);
for (const auto& scv : scvs(fvGeometry))
{
const auto& volVars = elemVolVars[scv];
- gradConcentration.axpy(volVars.concentration(), fluxVarCache.gradN(scv.indexInElement()));
- gradChemicalPotential.axpy(volVars.chemicalPotential(), fluxVarCache.gradN(scv.indexInElement()));
+ // v.axpy(a, w) means v <- v + a*w
+ gradConcentration.axpy(
+ volVars.concentration(),
+ fluxVarCache.gradN(scv.indexInElement())
+ );
+ gradChemicalPotential.axpy(
+ volVars.chemicalPotential(),
+ fluxVarCache.gradN(scv.indexInElement())
+ );
}
- const auto M = problem.mobility();
- const auto gamma = problem.surfaceTension();
-
NumEqVector flux;
- flux[Indices::massBalanceEqIdx] = -1.0*vtmv(scvf.unitOuterNormal(), M, gradChemicalPotential)*scvf.area();
- flux[Indices::chemicalPotentialEqIdx] = -1.0*vtmv(scvf.unitOuterNormal(), gamma, gradConcentration)*scvf.area();
+ // Compute the flux with `vtmv` (vector transposed times matrix times vector).
+ // The mobility and surface tension coefficients comes from the `problem` (see Part 2 of the example).
+ flux[Indices::massBalanceEqIdx] = -1.0*vtmv(
+ scvf.unitOuterNormal(), problem.mobility(), gradChemicalPotential
+ )*scvf.area();
+ flux[Indices::chemicalPotentialEqIdx] = -1.0*vtmv(
+ scvf.unitOuterNormal(), problem.surfaceTension(), gradConcentration
+ )*scvf.area();
return flux;
}
// [[/codeblock]]
- //
- // ### Source terms
- //
- // `computeSource` defines the sources terms at a sub control volume.
+
+ // **Source term:** The function `computeSource` computes the sources terms for a sub control volume.
// We implement a model-specific source term for the chemical potential equation before
// deferring further implementation to the problem where we add the derivative of the free
- // energy. The default implementation of this function also defers the calculation to the
- // problem.
+ // energy.
//
// [[codeblock]]
- /*!
- * \brief Calculate the source term of the equation
- */
NumEqVector computeSource(const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
@@ -242,78 +256,43 @@ public:
const SubControlVolume &scv) const
{
NumEqVector source(0.0);
-
+ source[Indices::massBalanceEqIdx] = 0.0;
source[Indices::chemicalPotentialEqIdx] = elemVolVars[scv].chemicalPotential();
-
// add contributions from problem (e.g. double well potential)
source += problem.source(element, fvGeometry, elemVolVars, scv);
-
return source;
}
};
-
} // end namespace Dumux
// [[/codeblock]]
// [[/content]]
//
-// ## Model properties/traits
+// ## 4. The model properties/traits
//
-// We set some general model traits and properties, using a TypeTag `CahnHilliardModel`.
-// For this type tag we define a `ModelTraits` struct and set a number of properties by specifying
-// further structs within the `Dumux::Properties` namespace.
+// In the `Dumux::Properties` namespace, we specialize properties for
+// the created type tag `CahnHilliardModel`.
//
// [[content]]
-//
-// ### Include the header for the property system and common properties and switch to the
-// `Properties` namespace.
-//
-// [[codeblock]]
-#include <dumux/common/properties.hh>
-
namespace Dumux::Properties {
-// [[/codeblock]]
-//
-// ### Define the type tag
-//
-// [[codeblock]]
-namespace TTag {
-struct CahnHilliardModel {};
-} // end namespace TTag
-// [[/codeblock]]
-//
-// ### Basic model properties
-//
-// Define some general properties to be used for this modedl such as a datatype for scalars and the
-// default vector for the primary variables. Here we can also use the model traits we specify below.
-//
-// [[codeblock]]
-//! Set the default type of scalar values to double
+
+// The type of the local residual is the class defined above.
template<class TypeTag>
-struct Scalar<TypeTag, TTag:: CahnHilliardModel >
-{ using type = double; };
+struct LocalResidual<TypeTag, TTag::CahnHilliardModel>
+{ using type = CahnHilliardModelLocalResidual<TypeTag>; };
-//! Set the default primary variable vector to a vector of size of number of equations
+// The default scalar type is double.
+// We compute with double precision floating point numbers.
template<class TypeTag>
-struct PrimaryVariables<TypeTag, TTag:: CahnHilliardModel >
-{
- using type = Dune::FieldVector<
- GetPropType<TypeTag, Properties::Scalar>,
- GetPropType<TypeTag, Properties::ModelTraits>::numEq()
- >;
-};
-// [[/codeblock]]
-//
-// ### Model traits
-//
-// We specify general traits of the implemented model, defining indices (often in `indices.hh`)
-// and the number of equations in the model.
-//
-// [[codeblock]]
-//! Set the model traits property
+struct Scalar<TypeTag, TTag::CahnHilliardModel>
+{ using type = double; };
+
+// The model traits specify some information about our equation system.
+// Here we have two equations. The indices allow to access primary variables
+// and equations with a named indices.
template<class TypeTag>
struct ModelTraits<TypeTag, TTag::CahnHilliardModel>
{
- struct Traits
+ struct type
{
struct Indices
{
@@ -326,22 +305,20 @@ struct ModelTraits<TypeTag, TTag::CahnHilliardModel>
static constexpr int numEq() { return 2; }
};
-
- using type = Traits;
};
-// [[/codeblock]]
-//
-// ### Further model properties
-//
-// Define further properties of the model, selecting the local residual and volumevariables defined
-// above.
-//
-// [[codeblock]]
+
+// The primary variable vector has entries of type `Scalar` and is
+// as large as the number of equations (here 2) but we keep it general.
template<class TypeTag>
-struct LocalResidual<TypeTag, TTag::CahnHilliardModel>
-{ using type = CahnHilliardModelLocalResidual<TypeTag>; };
+struct PrimaryVariables<TypeTag, TTag::CahnHilliardModel>
+{
+ using type = Dune::FieldVector<
+ GetPropType<TypeTag, Properties::Scalar>,
+ GetPropType<TypeTag, Properties::ModelTraits>::numEq()
+ >;
+};
-//! Set the volume variables property
+// Finally, the type of the volume variables is the class defined above.
template<class TypeTag>
struct VolumeVariables<TypeTag, TTag::CahnHilliardModel>
{
@@ -355,7 +332,7 @@ struct VolumeVariables<TypeTag, TTag::CahnHilliardModel>
};
} // end namespace Dumux::Properties
-// [[/codeblock]]
// [[/content]]
-
+// [[exclude]]
#endif
+// [[/exclude]]
diff --git a/examples/diffusion/doc/main.md b/examples/diffusion/doc/main.md
index 28d3bd539a..28791805ad 100644
--- a/examples/diffusion/doc/main.md
+++ b/examples/diffusion/doc/main.md
@@ -248,7 +248,7 @@ SolutionVector createInitialSolution(const GridGeometry& gg)
for (int n = 0; n < sol.size(); ++n)
sol[n] = dis(gen) + rank;
- // We, take the value of the processor with the minimum rank
+ // We take the value of the processor with the minimum rank
// and subtract the rank offset
if (gg.gridView().comm().size() > 1)
{
@@ -341,7 +341,7 @@ the problem for the boundary conditions, a solution vector and a grid variables
gridVariables->init(sol);
```
-Ww initialize the VTK output module and write out the initial concentration field
+We initialize the VTK output module and write out the initial concentration field
```cpp
VtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables, sol, problem->name());
diff --git a/examples/diffusion/main.cc b/examples/diffusion/main.cc
index fb80e86507..a3585e7eda 100644
--- a/examples/diffusion/main.cc
+++ b/examples/diffusion/main.cc
@@ -204,7 +204,7 @@ SolutionVector createInitialSolution(const GridGeometry& gg)
for (int n = 0; n < sol.size(); ++n)
sol[n] = dis(gen) + rank;
- // We, take the value of the processor with the minimum rank
+ // We take the value of the processor with the minimum rank
// and subtract the rank offset
if (gg.gridView().comm().size() > 1)
{
@@ -274,7 +274,7 @@ int main(int argc, char** argv)
auto gridVariables = std::make_shared<GridVariables>(problem, gridGeometry);
gridVariables->init(sol);
- // Ww initialize the VTK output module and write out the initial concentration field
+ // We initialize the VTK output module and write out the initial concentration field
VtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables, sol, problem->name());
vtkWriter.addVolumeVariable([](const auto& vv){ return vv.priVar(0); }, "c");
vtkWriter.write(0.0);
--
GitLab