diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
index fbf20b2458099cde5fd04af2561380804f40aedb..aacf28405d7293db3555bf49083ee4b425b9947c 100644
--- a/examples/CMakeLists.txt
+++ b/examples/CMakeLists.txt
@@ -1,6 +1,7 @@
add_subdirectory(2pinfiltration)
add_subdirectory(1ptracer)
add_subdirectory(biomineralization)
+add_subdirectory(cahn_hilliard)
add_subdirectory(diffusion)
add_subdirectory(embedded_network_1d3d)
add_subdirectory(shallowwaterfriction)
diff --git a/examples/README.md b/examples/README.md
index beb935118d8db8a5d64f7bcedeb2a8c74f17f499..029c5781ca6ac0dfa9839919b7ebad34b02770e8 100644
--- a/examples/README.md
+++ b/examples/README.md
@@ -22,12 +22,35 @@ __Model equations:__ A diffusion equation model fully developed and contained wi
__Discretization method:__ Vertex-centered finite volumes / control-volume finite elements (Lagrange, P1) (`BoxModel`)
</td>
-<td width="35%"><a href="diffusion/README.md">
+<td width="25%"><a href="diffusion/README.md">
<figure><img src="diffusion/img/diffusion.gif" alt="diffusion"/></figure></td>
</a></td>
</tr></table>
-### [:open_file_folder: Example 2: One-phase flow and tracer transport](1ptracer/README.md)
+### [:open_file_folder: Example 2: Cahn-Hilliard equation](cahn_hilliard/README.md)
+
+<table><tr><td>
+
+In this example we simulate the separation of two immiscible phases using the Cahn-Hilliard equation.
+
+You learn how to
+
+* setup a new nonlinear model equation (Cahn-Hilliard equation)
+* setup a model with two governing equations
+* implement a custom volume variables class to name variables
+* generate a randomly distributed vector-valued initial field (with MPI parallelism)
+* solve a time-dependent diffusion problem in parallel
+
+__Model equations:__ A Cahn-Hilliard equation model fully developed and contained within the example<br />
+__Discretization method:__ Vertex-centered finite volumes / control-volume finite elements (Lagrange, P1) (`BoxModel`)
+
+</td>
+<td width="25%"><a href="cahn_hilliard/README.md">
+<figure><img src="cahn_hilliard/img/animation.gif" alt="phase distribution" width="400px"/></figure></td>
+</a></td>
+</tr></table>
+
+### [:open_file_folder: Example 3: One-phase flow and tracer transport](1ptracer/README.md)
<table><tr><td>
@@ -51,7 +74,7 @@ __Discretization method:__ Cell-centered finite volumes with two-point flux appr
</a></td>
</tr></table>
-### [:open_file_folder: Example 3: Two-phase flow with infiltration and adaptive grid](2pinfiltration/README.md)
+### [:open_file_folder: Example 4: Two-phase flow with infiltration and adaptive grid](2pinfiltration/README.md)
<table><tr><td>
@@ -75,7 +98,7 @@ __Discretization method:__ Cell-centered finite volumes with two-point flux appr
</a></td>
</tr></table>
-### [:open_file_folder: Example 4: Shallow water model](shallowwaterfriction/README.md)
+### [:open_file_folder: Example 5: Shallow water model](shallowwaterfriction/README.md)
<table><tr><td>
@@ -95,7 +118,7 @@ __Discretization method:__ Cell-centered finite volumes with Riemann solver (`CC
</a></td>
</tr></table>
-### [:open_file_folder: Example 5: Freeflow channel](freeflowchannel/README.md)
+### [:open_file_folder: Example 6: Freeflow channel](freeflowchannel/README.md)
<table><tr><td>
@@ -115,7 +138,7 @@ __Discretization method:__ Finite volumes with staggered grid arrangement (`Stag
</a></td>
</tr></table>
-### [:open_file_folder: Example 6: One-phase flow with rotation-symmetric solution](1protationsymmetry/README.md)
+### [:open_file_folder: Example 7: One-phase flow with rotation-symmetric solution](1protationsymmetry/README.md)
<table><tr><td>
@@ -136,7 +159,7 @@ __Discretization method:__ Vertex-centered finite volumes / control-volume finit
</a></td>
</tr></table>
-### [:open_file_folder: Example 7: Biomineralization](biomineralization/README.md)
+### [:open_file_folder: Example 8: Biomineralization](biomineralization/README.md)
<table><tr><td>
@@ -161,7 +184,7 @@ __Discretization method:__ Vertex-centered finite volumes / control-volume finit
</a></td>
</tr></table>
-### [:open_file_folder: Example 8: Lid-driven cavity](liddrivencavity/README.md)
+### [:open_file_folder: Example 9: Lid-driven cavity](liddrivencavity/README.md)
<table><tr><td>
@@ -182,7 +205,7 @@ __Discretization method:__ Finite volumes with staggered grid arrangement (`Stag
</a></td>
</tr></table>
-### [:open_file_folder: Example 9: Permeability estimation using a pore-network model](porenetwork_upscaling/README.md)
+### [:open_file_folder: Example 10: Permeability estimation using a pore-network model](porenetwork_upscaling/README.md)
<table><tr><td>
@@ -203,7 +226,7 @@ __Discretization method:__ Pore-network (`PoreNetworkModel`)
</a></td>
</tr></table>
-### [:open_file_folder: Example 10: Embedded network 1D-3D model (tissue perfusion)](embedded_network_1d3d/README.md)
+### [:open_file_folder: Example 11: Embedded network 1D-3D model (tissue perfusion)](embedded_network_1d3d/README.md)
<table><tr><td>
diff --git a/examples/cahn_hilliard/.doc_config b/examples/cahn_hilliard/.doc_config
new file mode 100644
index 0000000000000000000000000000000000000000..6df89661cd20ce69ae526ae0f44fd3baf4a89b5e
--- /dev/null
+++ b/examples/cahn_hilliard/.doc_config
@@ -0,0 +1,27 @@
+{
+ "README.md" : [
+ "doc/_intro.md"
+ ],
+
+ "doc/model.md" : [
+ "doc/_model.md",
+ "model.hh"
+ ],
+
+ "doc/main.md" : [
+ "doc/_main.md",
+ "main.cc"
+ ],
+
+ "navigation" : {
+ "mainpage" : "README.md",
+ "subpages" : [
+ "doc/model.md",
+ "doc/main.md"
+ ],
+ "subtitles" : [
+ "Model implementation",
+ "Main program"
+ ]
+ }
+}
diff --git a/examples/cahn_hilliard/CMakeLists.txt b/examples/cahn_hilliard/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..f69a567fda9d580add84414e8939f4a7bec35582
--- /dev/null
+++ b/examples/cahn_hilliard/CMakeLists.txt
@@ -0,0 +1,27 @@
+dune_symlink_to_source_files(FILES "params.input")
+
+# A simple test to automatically check the example compiles and produces the same result
+dumux_add_test(NAME example_cahn_hilliard
+ LABELS example
+ SOURCES main.cc
+ COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+ CMD_ARGS --script fuzzy
+ --files ${CMAKE_SOURCE_DIR}/test/references/example_cahn_hilliard-reference.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/cahn_hilliard-00102.vtu
+ --command "${CMAKE_CURRENT_BINARY_DIR}/example_cahn_hilliard params.input
+ -Grid.Cells \" 10 10\" -Problem.Mobility 0.01 -Problem.SurfaceTension 0.5")
+
+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
+ ${CMAKE_CURRENT_BINARY_DIR}/s0002-p0000-cahn_hilliard_parallel-00102.vtu
+ ${CMAKE_SOURCE_DIR}/test/references/example_cahn_hilliard_parallel_p1-reference.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/s0002-p0001-cahn_hilliard_parallel-00102.vtu
+ --command "${MPIEXEC} -np 2
+ ${CMAKE_CURRENT_BINARY_DIR}/example_cahn_hilliard params.input
+ -Problem.Name cahn_hilliard_parallel
+ -Grid.Cells \" 10 10\" -Problem.Mobility 0.01 -Problem.SurfaceTension 0.5")
diff --git a/examples/cahn_hilliard/README.md b/examples/cahn_hilliard/README.md
new file mode 100644
index 0000000000000000000000000000000000000000..11552d371170e8184917058c461b43e2a353c6d5
--- /dev/null
+++ b/examples/cahn_hilliard/README.md
@@ -0,0 +1,127 @@
+<!-- Important: This file has been automatically generated by generate_example_docs.py. Do not edit this file directly! -->
+
+# Cahn-Hilliard equation
+
+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__
+
+* 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
+
+__Results__. The result will look like below.
+
+<figure>
+ <center>
+ <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>
+
+__Table of contents__. This description is structured as follows:
+
+[TOC]
+
+## Equation and problem description
+
+The Cahn-Hilliard equation is a forth order partial differential equation (PDE) and reads
+
+```math
+\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]
+```
+
+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.
+
+The fourth order PDE cannot be solved by a standard finite volume scheme. We therefore
+split the equation in two second order PDEs
+
+```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,
+
+```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}.
+```
+
+### 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/model.md) |
+|---:|
+
+
+## Part 2: Main program
+
+| [: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
new file mode 100644
index 0000000000000000000000000000000000000000..42f94f6c00ef5c09d5439112abc854ce91748f5a
--- /dev/null
+++ b/examples/cahn_hilliard/doc/_intro.md
@@ -0,0 +1,114 @@
+# Cahn-Hilliard equation
+
+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__
+
+* 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
+
+__Results__. The result will look like below.
+
+<figure>
+ <center>
+ <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>
+
+__Table of contents__. This description is structured as follows:
+
+[TOC]
+
+## Equation and problem description
+
+The Cahn-Hilliard equation is a forth order partial differential equation (PDE) and reads
+
+```math
+\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]
+```
+
+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.
+
+The fourth order PDE cannot be solved by a standard finite volume scheme. We therefore
+split the equation in two second order PDEs
+
+```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,
+
+```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}.
+```
+
+### 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/_main.md b/examples/cahn_hilliard/doc/_main.md
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/examples/cahn_hilliard/doc/_model.md b/examples/cahn_hilliard/doc/_model.md
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/examples/cahn_hilliard/doc/main.md b/examples/cahn_hilliard/doc/main.md
new file mode 100644
index 0000000000000000000000000000000000000000..778b01e7616ba460060b8ff8eb9f0db9282f88f4
--- /dev/null
+++ b/examples/cahn_hilliard/doc/main.md
@@ -0,0 +1,402 @@
+<!-- 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](model.md) |
+|---|---:|
+
+
+
+# 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
+
+__Table of contents__
+
+[TOC]
+
+We start in `main.cc` with the necessary header includes:
+<details><summary> Click to show includes</summary>
+
+```cpp
+
+#include <config.h>
+#include <type_traits>
+#include <random>
+
+#include <dumux/common/initialize.hh>
+#include <dumux/common/properties.hh>
+#include <dumux/common/parameters.hh>
+#include <dumux/common/numeqvector.hh>
+#include <dumux/common/boundarytypes.hh>
+#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 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 boilerplate code</summary>
+
+```cpp
+
+ using ParentType = FVProblem<TypeTag>;
+ using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
+ using FVElementGeometry = typename GridGeometry::LocalView;
+ using SubControlVolume = typename GridGeometry::SubControlVolume;
+ using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
+ using NumEqVector = Dumux::NumEqVector<PrimaryVariables>;
+ using BoundaryTypes = Dumux::BoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>;
+ using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
+```
+
+</details>
+In the constructor, we read the parameter constants from the
+parameter tree (which is initialized with the content of `params.input`).
+
+```cpp
+public:
+ CahnHilliardTestProblem(std::shared_ptr<const GridGeometry> gridGeometry)
+ : ParentType(gridGeometry)
+ {
+ mobility_ = getParam<Scalar>("Problem.Mobility");
+ surfaceTension_ = getParam<Scalar>("Problem.SurfaceTension");
+ energyScale_ = getParam<Scalar>("Problem.EnergyScale");
+ }
+```
+
+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>
+ NumEqVector source(const Element &element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolume &scv) const
+ {
+ NumEqVector values(0.0);
+ const auto& c = elemVolVars[scv].concentration();
+ values[Indices::chemicalPotentialEqIdx] = -energyScale_*2.0*c*(2.0*c*c - 3*c + 1);
+ return values;
+ }
+```
+
+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
+ {
+ BoundaryTypes values;
+ values.setAllNeumann();
+ return values;
+ }
+
+ NumEqVector neumannAtPos(const GlobalPosition& globalPos) const
+ { return { 0.0, 0.0 }; }
+```
+
+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
+ { return mobility_; }
+
+ Scalar surfaceTension() const
+ { return surfaceTension_; }
+
+private:
+ Scalar mobility_;
+ Scalar surfaceTension_;
+ Scalar energyScale_;
+};
+} // end namespace Dumux
+```
+
+
+</details>
+
+## 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>
+
+
+```cpp
+namespace Dumux::Properties::TTag {
+struct CahnHilliardTest
+{
+ using InheritsFrom = std::tuple<CahnHilliardModel, BoxModel>;
+
+ using Grid = Dune::YaspGrid<2>;
+
+ 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>
+
+## 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>
+
+
+```cpp
+// functor for data communication with MPI
+struct MinScatter
+{
+ template<class A, class B>
+ static void apply(A& a, const B& b)
+ { a[0] = std::min(a[0], b[0]); }
+};
+
+// create the random initial solution
+template<class SolutionVector, class GridGeometry>
+SolutionVector createInitialSolution(const GridGeometry& gg)
+{
+ SolutionVector sol(gg.numDofs());
+
+ // Generate random number and add processor offset
+ // For sequential runs `rank` always returns `0`.
+ std::mt19937 gen(0); // seed is 0 for deterministic results
+ std::uniform_real_distribution<> dis(0.0, 1.0);
+ for (int n = 0; n < sol.size(); ++n)
+ {
+ sol[n][0] = 0.42 + 0.02*(0.5-dis(gen)) + gg.gridView().comm().rank();
+ sol[n][1] = 0.0;
+ }
+
+ // We take the value of the processor with the minimum rank
+ // and subtract the rank offset
+ if (gg.gridView().comm().size() > 1)
+ {
+ Dumux::VectorCommDataHandle<
+ typename GridGeometry::VertexMapper,
+ SolutionVector,
+ GridGeometry::GridView::dimension,
+ MinScatter
+ > minHandle(gg.vertexMapper(), sol);
+ gg.gridView().communicate(minHandle, Dune::All_All_Interface, Dune::ForwardCommunication);
+
+ // Remove processor offset
+ for (int n = 0; n < sol.size(); ++n)
+ sol[n][0] -= std::floor(sol[n][0]);
+ }
+ return sol;
+}
+```
+
+
+</details>
+
+## 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;
+```
+
+We initialize MPI and/or multithreading backend. When not running
+in parallel the MPI setup is skipped.
+
+```cpp
+ Dumux::initialize(argc, argv);
+```
+
+Then we initialize parameter tree.
+
+```cpp
+ Parameters::init(argc, argv);
+```
+
+We create an alias for the type tag for this problem
+and extract type information through properties.
+
+```cpp
+ 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>;
+ using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
+ using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
+```
+
+We initialize the grid, grid geometry, problem, solution vector, and grid variables.
+
+```cpp
+ 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);
+```
+
+We initialize the VTK output module and write out the initial concentration and chemical potential
+
+```cpp
+ 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);
+```
+
+We instantiate time loop using start and end time as well as
+the time step size from the parameter tree (`params.input`)
+
+```cpp
+ 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.
+
+```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
+ 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);
+```
+
+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
+ timeLoop->start(); do
+ {
+ // 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
+ oldSol = sol;
+ gridVariables->advanceTimeStep();
+
+ // advance to the time loop to the next step
+ timeLoop->advanceTimeStep();
+
+ // write VTK output (concentration field and chemical potential)
+ vtkWriter.write(timeLoop->time());
+
+ // report statistics of this time step
+ timeLoop->reportTimeStep();
+
+ // set new dt as suggested by the Newton solver
+ timeLoop->setTimeStepSize(solver.suggestTimeStepSize(timeLoop->timeStepSize()));
+
+ } while (!timeLoop->finished());
+
+ // print the final report
+ timeLoop->finalize(gridGeometry->gridView().comm());
+ return 0;
+}
+```
+
+
+</details>
+
+
+| [: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/model.md b/examples/cahn_hilliard/doc/model.md
new file mode 100644
index 0000000000000000000000000000000000000000..b9fc2a3d73c708ec5b98ca65d705200c4c4af265
--- /dev/null
+++ b/examples/cahn_hilliard/doc/model.md
@@ -0,0 +1,357 @@
+<!-- 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](main.md) |
+|---|---:|
+
+
+
+# Cahn-Hilliard equation model definition
+
+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
+
+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>
+
+
+```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
+{
+ using Scalar = typename Traits::PrimaryVariables::value_type;
+ static_assert(Traits::PrimaryVariables::dimension == Traits::ModelTraits::numEq());
+public:
+ //! export the type used for the primary variables
+ using PrimaryVariables = typename Traits::PrimaryVariables;
+ //! export the indices type
+ using Indices = typename Traits::ModelTraits::Indices;
+```
+
+
+**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
+ //! 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,
+ const Element& element,
+ const SubControlVolume& scv)
+ {
+ priVars_ = elemSol[scv.indexInElement()];
+ }
+```
+
+
+**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
+ { return priVars_[Indices::concentrationIdx]; }
+
+ Scalar chemicalPotential() const
+ { return priVars_[Indices::chemicalPotentialIdx]; }
+
+ Scalar priVar(const int pvIdx) const
+ { return priVars_[pvIdx]; }
+
+ const PrimaryVariables& priVars() const
+ { return priVars_; }
+
+ Scalar extrusionFactor() const
+ { return 1.0; }
+
+private:
+ PrimaryVariables priVars_;
+};
+} // end namespace Dumux
+```
+
+
+</details>
+
+## 3. The local residual
+
+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>
+
+
+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>
+{
+```
+
+
+**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
+ NumEqVector computeStorage(const Problem& problem,
+ const SubControlVolume& scv,
+ const VolumeVariables& volVars) const
+ {
+ NumEqVector storage;
+ storage[Indices::massBalanceEqIdx] = volVars.concentration();
+ storage[Indices::chemicalPotentialEqIdx] = 0.0;
+ return storage;
+ }
+```
+
+**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
+ NumEqVector computeFlux(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ 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];
+ // 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())
+ );
+ }
+
+ NumEqVector flux;
+ // 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 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.
+
+
+```cpp
+ NumEqVector computeSource(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ 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>
+
+## 4. The model properties/traits
+
+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>
+
+
+```cpp
+namespace Dumux::Properties {
+```
+
+The type of the local residual is the class defined above.
+
+```cpp
+template<class TypeTag>
+struct LocalResidual<TypeTag, TTag::CahnHilliardModel>
+{ using type = CahnHilliardModelLocalResidual<TypeTag>; };
+```
+
+The default scalar type is double.
+We compute with double precision floating point numbers.
+
+```cpp
+template<class TypeTag>
+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.
+
+```cpp
+template<class TypeTag>
+struct ModelTraits<TypeTag, TTag::CahnHilliardModel>
+{
+ struct type
+ {
+ struct Indices
+ {
+ static constexpr int concentrationIdx = 0;
+ static constexpr int chemicalPotentialIdx = 1;
+
+ static constexpr int massBalanceEqIdx = 0;
+ static constexpr int chemicalPotentialEqIdx = 1;
+ };
+
+ static constexpr int numEq() { return 2; }
+ };
+};
+```
+
+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 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.
+
+```cpp
+template<class TypeTag>
+struct VolumeVariables<TypeTag, TTag::CahnHilliardModel>
+{
+ struct Traits
+ {
+ using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
+ using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
+ };
+
+ using type = CahnHilliardModelVolumeVariables<Traits>;
+};
+
+} // end namespace Dumux::Properties
+```
+
+
+</details>
+
+
+| [:arrow_left: Back to the main documentation](../README.md) | [:arrow_right: Continue with part 2](main.md) |
+|---|---:|
+
diff --git a/examples/cahn_hilliard/img/animation.gif b/examples/cahn_hilliard/img/animation.gif
new file mode 100644
index 0000000000000000000000000000000000000000..6e266f5f8f146c7780dea6f6c1f4ddb950b49d46
Binary files /dev/null and b/examples/cahn_hilliard/img/animation.gif differ
diff --git a/examples/cahn_hilliard/main.cc b/examples/cahn_hilliard/main.cc
new file mode 100644
index 0000000000000000000000000000000000000000..b91a6d9acbcb8e692f655330baead3f2f60f540e
--- /dev/null
+++ b/examples/cahn_hilliard/main.cc
@@ -0,0 +1,336 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+
+// # 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
+//
+// __Table of contents__
+//
+// [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>
+#include <dumux/common/numeqvector.hh>
+#include <dumux/common/boundarytypes.hh>
+#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 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>
+{
+ // [[details]] boilerplate code
+ using ParentType = FVProblem<TypeTag>;
+ using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
+ using FVElementGeometry = typename GridGeometry::LocalView;
+ using SubControlVolume = typename GridGeometry::SubControlVolume;
+ using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
+ using NumEqVector = Dumux::NumEqVector<PrimaryVariables>;
+ using BoundaryTypes = Dumux::BoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>;
+ using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
+ // [[/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)
+ {
+ mobility_ = getParam<Scalar>("Problem.Mobility");
+ surfaceTension_ = getParam<Scalar>("Problem.SurfaceTension");
+ energyScale_ = getParam<Scalar>("Problem.EnergyScale");
+ }
+
+ // 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,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolume &scv) const
+ {
+ NumEqVector values(0.0);
+ const auto& c = elemVolVars[scv].concentration();
+ values[Indices::chemicalPotentialEqIdx] = -energyScale_*2.0*c*(2.0*c*c - 3*c + 1);
+ return values;
+ }
+
+ // 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
+ {
+ BoundaryTypes values;
+ values.setAllNeumann();
+ return values;
+ }
+
+ NumEqVector neumannAtPos(const GlobalPosition& globalPos) const
+ { return { 0.0, 0.0 }; }
+ // [[/codeblock]]
+
+ // 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_; }
+
+ Scalar surfaceTension() const
+ { return surfaceTension_; }
+
+private:
+ Scalar mobility_;
+ Scalar surfaceTension_;
+ Scalar energyScale_;
+};
+} // end namespace Dumux
+// [[/codeblock]]
+// [[/content]]
+
+//
+// ## 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)).
+//
+// [[content]]
+namespace Dumux::Properties::TTag {
+struct CahnHilliardTest
+{
+ using InheritsFrom = std::tuple<CahnHilliardModel, BoxModel>;
+
+ using Grid = Dune::YaspGrid<2>;
+
+ 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
+// [[/content]]
+//
+// ## 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).
+//
+// [[content]]
+// [[codeblock]]
+// functor for data communication with MPI
+struct MinScatter
+{
+ template<class A, class B>
+ static void apply(A& a, const B& b)
+ { a[0] = std::min(a[0], b[0]); }
+};
+
+// create the random initial solution
+template<class SolutionVector, class GridGeometry>
+SolutionVector createInitialSolution(const GridGeometry& gg)
+{
+ SolutionVector sol(gg.numDofs());
+
+ // Generate random number and add processor offset
+ // For sequential runs `rank` always returns `0`.
+ std::mt19937 gen(0); // seed is 0 for deterministic results
+ std::uniform_real_distribution<> dis(0.0, 1.0);
+ for (int n = 0; n < sol.size(); ++n)
+ {
+ sol[n][0] = 0.42 + 0.02*(0.5-dis(gen)) + gg.gridView().comm().rank();
+ sol[n][1] = 0.0;
+ }
+
+ // We take the value of the processor with the minimum rank
+ // and subtract the rank offset
+ if (gg.gridView().comm().size() > 1)
+ {
+ Dumux::VectorCommDataHandle<
+ typename GridGeometry::VertexMapper,
+ SolutionVector,
+ GridGeometry::GridView::dimension,
+ MinScatter
+ > minHandle(gg.vertexMapper(), sol);
+ gg.gridView().communicate(minHandle, Dune::All_All_Interface, Dune::ForwardCommunication);
+
+ // Remove processor offset
+ for (int n = 0; n < sol.size(); ++n)
+ sol[n][0] -= std::floor(sol[n][0]);
+ }
+ return sol;
+}
+// [[/codeblock]]
+// [[/content]]
+//
+// ## 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.
+//
+// [[content]]
+int main(int argc, char** argv)
+{
+ using namespace Dumux;
+
+ // We initialize MPI and/or multithreading backend. When not running
+ // in parallel the MPI setup is skipped.
+ Dumux::initialize(argc, argv);
+
+ // Then we initialize parameter tree.
+ Parameters::init(argc, argv);
+
+ // We create an alias for the type tag for this problem
+ // and extract type information through properties.
+ 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>;
+ using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
+ 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);
+
+ // 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);
+
+ // 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")
+ );
+
+ // We set the maximum time step size allowed in the adaptive time stepping scheme.
+ timeLoop->setMaxTimeStepSize(getParam<Scalar>("TimeLoop.MaxTimeStepSize"));
+
+ // 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);
+
+ // 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]]
+ timeLoop->start(); do
+ {
+ // 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
+ oldSol = sol;
+ gridVariables->advanceTimeStep();
+
+ // advance to the time loop to the next step
+ timeLoop->advanceTimeStep();
+
+ // write VTK output (concentration field and chemical potential)
+ vtkWriter.write(timeLoop->time());
+
+ // report statistics of this time step
+ timeLoop->reportTimeStep();
+
+ // set new dt as suggested by the Newton solver
+ timeLoop->setTimeStepSize(solver.suggestTimeStepSize(timeLoop->timeStepSize()));
+
+ } while (!timeLoop->finished());
+
+ // print the final report
+ timeLoop->finalize(gridGeometry->gridView().comm());
+ return 0;
+}
+// [[/codeblock]]
+// [[/content]]
diff --git a/examples/cahn_hilliard/model.hh b/examples/cahn_hilliard/model.hh
new file mode 100644
index 0000000000000000000000000000000000000000..f77c3d4f3f740e2c5cc6f146e7be61a3f3b1175f
--- /dev/null
+++ b/examples/cahn_hilliard/model.hh
@@ -0,0 +1,338 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+
+#ifndef DUMUX_EXAMPLES_CAHN_HILLIARD_MODEL_HH
+#define DUMUX_EXAMPLES_CAHN_HILLIARD_MODEL_HH
+
+// # Cahn-Hilliard equation model definition
+//
+// 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
+//
+// 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]]
+//
+// 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
+{
+ using Scalar = typename Traits::PrimaryVariables::value_type;
+ static_assert(Traits::PrimaryVariables::dimension == Traits::ModelTraits::numEq());
+public:
+ //! export the type used for the primary variables
+ using PrimaryVariables = typename Traits::PrimaryVariables;
+ //! export the indices type
+ using Indices = typename Traits::ModelTraits::Indices;
+// [[/codeblock]]
+ //
+ // **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]]
+ //! 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,
+ const Element& element,
+ const SubControlVolume& scv)
+ {
+ priVars_ = elemSol[scv.indexInElement()];
+ }
+ // [[/codeblock]]
+ //
+ // **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]; }
+
+ Scalar chemicalPotential() const
+ { return priVars_[Indices::chemicalPotentialIdx]; }
+
+ Scalar priVar(const int pvIdx) const
+ { return priVars_[pvIdx]; }
+
+ const PrimaryVariables& priVars() const
+ { return priVars_; }
+
+ Scalar extrusionFactor() const
+ { return 1.0; }
+
+private:
+ PrimaryVariables priVars_;
+};
+} // end namespace Dumux
+// [[/codeblock]]
+// [[/content]]
+//
+// ## 3. The local residual
+//
+// 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]]
+//
+// 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`.
+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>;
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using Problem = GetPropType<TypeTag, Properties::Problem>;
+ using NumEqVector = Dumux::NumEqVector<GetPropType<TypeTag, Properties::PrimaryVariables>>;
+
+ 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 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>;
+ using Indices = typename ModelTraits::Indices;
+ static constexpr int dimWorld = GridView::dimensionworld;
+public:
+ using ParentType::ParentType;
+ // [[/exclude]]
+ //
+ // **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]]
+ NumEqVector computeStorage(const Problem& problem,
+ const SubControlVolume& scv,
+ const VolumeVariables& volVars) const
+ {
+ NumEqVector storage;
+ storage[Indices::massBalanceEqIdx] = volVars.concentration();
+ storage[Indices::chemicalPotentialEqIdx] = 0.0;
+ return storage;
+ }
+ // [[/codeblock]]
+
+ // **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]]
+ NumEqVector computeFlux(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ 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];
+ // 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())
+ );
+ }
+
+ NumEqVector flux;
+ // 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 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.
+ //
+ // [[codeblock]]
+ NumEqVector computeSource(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ 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]]
+//
+// ## 4. The model properties/traits
+//
+// In the `Dumux::Properties` namespace, we specialize properties for
+// the created type tag `CahnHilliardModel`.
+//
+// [[content]]
+namespace Dumux::Properties {
+
+// The type of the local residual is the class defined above.
+template<class TypeTag>
+struct LocalResidual<TypeTag, TTag::CahnHilliardModel>
+{ using type = CahnHilliardModelLocalResidual<TypeTag>; };
+
+// The default scalar type is double.
+// We compute with double precision floating point numbers.
+template<class TypeTag>
+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 type
+ {
+ struct Indices
+ {
+ static constexpr int concentrationIdx = 0;
+ static constexpr int chemicalPotentialIdx = 1;
+
+ static constexpr int massBalanceEqIdx = 0;
+ static constexpr int chemicalPotentialEqIdx = 1;
+ };
+
+ static constexpr int numEq() { return 2; }
+ };
+};
+
+// 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 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.
+template<class TypeTag>
+struct VolumeVariables<TypeTag, TTag::CahnHilliardModel>
+{
+ struct Traits
+ {
+ using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
+ using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
+ };
+
+ using type = CahnHilliardModelVolumeVariables<Traits>;
+};
+
+} // end namespace Dumux::Properties
+// [[/content]]
+// [[exclude]]
+#endif
+// [[/exclude]]
diff --git a/examples/cahn_hilliard/params.input b/examples/cahn_hilliard/params.input
new file mode 100644
index 0000000000000000000000000000000000000000..f9d41a960d4dd43dd177449a3ffdba6ea82425c4
--- /dev/null
+++ b/examples/cahn_hilliard/params.input
@@ -0,0 +1,19 @@
+[TimeLoop]
+TEnd = 1.0
+InitialTimeStepSize = 0.0025
+MaxTimeStepSize = 0.01
+
+[Grid]
+LowerLeft = 0 0
+UpperRight = 1 1
+Cells = 100 100
+Overlap = 2
+
+[Problem]
+Name = cahn_hilliard
+Mobility = 0.0001
+SurfaceTension = 0.01
+EnergyScale = 100
+
+[Newton]
+EnablePartialReassembly = true
diff --git a/examples/diffusion/doc/main.md b/examples/diffusion/doc/main.md
index 28d3bd539ab5ae4faf80e587e086cb8d4ef1a25c..28791805adf788e9b36959401d4423087b806ac6 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 fb80e86507353b649afc9a898c5e105636503af5..a3585e7eda6afad123933362147c9aeecd0bbbda 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);
diff --git a/test/references/example_cahn_hilliard-reference.vtu b/test/references/example_cahn_hilliard-reference.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..09646223e7ad6c42a05c1f58c59541c0eb1bd1bd
--- /dev/null
+++ b/test/references/example_cahn_hilliard-reference.vtu
@@ -0,0 +1,143 @@
+<?xml version="1.0"?>
+<VTKFile type="UnstructuredGrid" version="0.1" byte_order="LittleEndian">
+ <UnstructuredGrid>
+ <Piece NumberOfCells="100" NumberOfPoints="121">
+ <PointData Scalars="c">
+ <DataArray type="Float32" Name="c" NumberOfComponents="1" format="ascii">
+ 0.0147483 0.0153922 0.0152169 0.0159215 0.0211903 0.0219842 0.0574734 0.059312 0.243869 0.249797 0.780021 0.785981
+ 0.933771 0.937786 0.919201 0.930348 0.727812 0.794766 0.217154 0.371188 0.089294 0.205316 0.0165708 0.0174294
+ 0.0241264 0.0638759 0.264274 0.799598 0.947169 0.952779 0.895596 0.777857 0.689978 0.0188589 0.0198609 0.0270257
+ 0.0679906 0.274159 0.808361 0.955311 0.970942 0.952698 0.923767 0.908282 0.0236373 0.0244423 0.0307829 0.0674349
+ 0.253451 0.784383 0.949036 0.976182 0.970998 0.960571 0.955706 0.0420799 0.0405039 0.0407857 0.063462 0.187209
+ 0.640986 0.887371 0.944205 0.948217 0.940673 0.936661 0.140059 0.120251 0.08502 0.0715591 0.116467 0.284734
+ 0.626543 0.764201 0.781002 0.767473 0.75901 0.57634 0.471128 0.260339 0.131959 0.0887285 0.112697 0.177008
+ 0.230102 0.238261 0.225351 0.218237 0.913277 0.880901 0.73408 0.319931 0.1164 0.0626875 0.0578299 0.059922
+ 0.0573087 0.0520993 0.0496347 0.995737 0.983019 0.91484 0.650198 0.178748 0.0581001 0.0342137 0.027712 0.0233939
+ 0.0199908 0.0186564 1.00917 1.00078 0.952561 0.748489 0.215301 0.06028 0.0307777 0.0230179 0.0186793 0.0156897
+ 0.0145771
+ </DataArray>
+ <DataArray type="Float32" Name="mu" NumberOfComponents="1" format="ascii">
+ 2.70616 2.62202 2.78879 2.71354 2.34932 2.46783 1.82308 1.98541 0.912757 1.12665 -0.645811 -0.389002
+ -2.94258 -2.62547 -5.9443 -5.52026 -9.56572 -8.97809 -12.1828 -12.2064 -12.9292 -12.8866 3.01891 2.96958
+ 2.80417 2.45808 1.7679 0.378813 -1.70733 -4.34003 -7.36216 -10.3029 -11.8652 3.34487 3.33489 3.29487
+ 3.17847 2.80737 1.63606 -0.278211 -2.63278 -5.13046 -7.2605 -8.15371 3.68841 3.7232 3.82718 3.98948
+ 4.08821 3.32235 1.53081 -0.675997 -2.84224 -4.49376 -5.12372 3.9604 4.03591 4.26374 4.65125 5.15866
+ 5.25912 3.48158 1.27611 -0.847849 -2.39324 -2.96585 4.08588 4.19506 4.50428 4.98384 5.58717 6.07524
+ 5.20498 2.9048 0.717505 -0.890181 -1.49583 4.07639 4.19553 4.54804 4.99335 5.43628 5.60167 5.04849
+ 3.62242 1.85036 0.491861 -0.00883526 4.19789 4.27215 4.48041 4.81194 5.0771 5.07193 4.60785 3.6696
+ 2.53143 1.62376 1.28207 4.26478 4.30765 4.42553 4.58578 4.78349 4.73602 4.35412 3.67582 2.88825
+ 2.26172 2.02395 4.28205 4.31461 4.3997 4.51565 4.67059 4.62439 4.27705 3.68104 3.00106 2.46303
+ 2.25887
+ </DataArray>
+ </PointData>
+ <CellData Scalars="process rank">
+ <DataArray type="Float32" Name="process rank" NumberOfComponents="1" format="ascii">
+ 0 0 0 0 0 0 0 0 0 0 0 0
+ 0 0 0 0 0 0 0 0 0 0 0 0
+ 0 0 0 0 0 0 0 0 0 0 0 0
+ 0 0 0 0 0 0 0 0 0 0 0 0
+ 0 0 0 0 0 0 0 0 0 0 0 0
+ 0 0 0 0 0 0 0 0 0 0 0 0
+ 0 0 0 0 0 0 0 0 0 0 0 0
+ 0 0 0 0 0 0 0 0 0 0 0 0
+ 0 0 0 0
+ </DataArray>
+ </CellData>
+ <Points>
+ <DataArray type="Float32" Name="Coordinates" NumberOfComponents="3" format="ascii">
+ 0 0 0 0.1 0 0 0 0.1 0 0.1 0.1 0
+ 0.2 0 0 0.2 0.1 0 0.3 0 0 0.3 0.1 0
+ 0.4 0 0 0.4 0.1 0 0.5 0 0 0.5 0.1 0
+ 0.6 0 0 0.6 0.1 0 0.7 0 0 0.7 0.1 0
+ 0.8 0 0 0.8 0.1 0 0.9 0 0 0.9 0.1 0
+ 1 0 0 1 0.1 0 0 0.2 0 0.1 0.2 0
+ 0.2 0.2 0 0.3 0.2 0 0.4 0.2 0 0.5 0.2 0
+ 0.6 0.2 0 0.7 0.2 0 0.8 0.2 0 0.9 0.2 0
+ 1 0.2 0 0 0.3 0 0.1 0.3 0 0.2 0.3 0
+ 0.3 0.3 0 0.4 0.3 0 0.5 0.3 0 0.6 0.3 0
+ 0.7 0.3 0 0.8 0.3 0 0.9 0.3 0 1 0.3 0
+ 0 0.4 0 0.1 0.4 0 0.2 0.4 0 0.3 0.4 0
+ 0.4 0.4 0 0.5 0.4 0 0.6 0.4 0 0.7 0.4 0
+ 0.8 0.4 0 0.9 0.4 0 1 0.4 0 0 0.5 0
+ 0.1 0.5 0 0.2 0.5 0 0.3 0.5 0 0.4 0.5 0
+ 0.5 0.5 0 0.6 0.5 0 0.7 0.5 0 0.8 0.5 0
+ 0.9 0.5 0 1 0.5 0 0 0.6 0 0.1 0.6 0
+ 0.2 0.6 0 0.3 0.6 0 0.4 0.6 0 0.5 0.6 0
+ 0.6 0.6 0 0.7 0.6 0 0.8 0.6 0 0.9 0.6 0
+ 1 0.6 0 0 0.7 0 0.1 0.7 0 0.2 0.7 0
+ 0.3 0.7 0 0.4 0.7 0 0.5 0.7 0 0.6 0.7 0
+ 0.7 0.7 0 0.8 0.7 0 0.9 0.7 0 1 0.7 0
+ 0 0.8 0 0.1 0.8 0 0.2 0.8 0 0.3 0.8 0
+ 0.4 0.8 0 0.5 0.8 0 0.6 0.8 0 0.7 0.8 0
+ 0.8 0.8 0 0.9 0.8 0 1 0.8 0 0 0.9 0
+ 0.1 0.9 0 0.2 0.9 0 0.3 0.9 0 0.4 0.9 0
+ 0.5 0.9 0 0.6 0.9 0 0.7 0.9 0 0.8 0.9 0
+ 0.9 0.9 0 1 0.9 0 0 1 0 0.1 1 0
+ 0.2 1 0 0.3 1 0 0.4 1 0 0.5 1 0
+ 0.6 1 0 0.7 1 0 0.8 1 0 0.9 1 0
+ 1 1 0
+ </DataArray>
+ </Points>
+ <Cells>
+ <DataArray type="Int32" Name="connectivity" NumberOfComponents="1" format="ascii">
+ 0 1 3 2 1 4 5 3 4 6 7 5
+ 6 8 9 7 8 10 11 9 10 12 13 11
+ 12 14 15 13 14 16 17 15 16 18 19 17
+ 18 20 21 19 2 3 23 22 3 5 24 23
+ 5 7 25 24 7 9 26 25 9 11 27 26
+ 11 13 28 27 13 15 29 28 15 17 30 29
+ 17 19 31 30 19 21 32 31 22 23 34 33
+ 23 24 35 34 24 25 36 35 25 26 37 36
+ 26 27 38 37 27 28 39 38 28 29 40 39
+ 29 30 41 40 30 31 42 41 31 32 43 42
+ 33 34 45 44 34 35 46 45 35 36 47 46
+ 36 37 48 47 37 38 49 48 38 39 50 49
+ 39 40 51 50 40 41 52 51 41 42 53 52
+ 42 43 54 53 44 45 56 55 45 46 57 56
+ 46 47 58 57 47 48 59 58 48 49 60 59
+ 49 50 61 60 50 51 62 61 51 52 63 62
+ 52 53 64 63 53 54 65 64 55 56 67 66
+ 56 57 68 67 57 58 69 68 58 59 70 69
+ 59 60 71 70 60 61 72 71 61 62 73 72
+ 62 63 74 73 63 64 75 74 64 65 76 75
+ 66 67 78 77 67 68 79 78 68 69 80 79
+ 69 70 81 80 70 71 82 81 71 72 83 82
+ 72 73 84 83 73 74 85 84 74 75 86 85
+ 75 76 87 86 77 78 89 88 78 79 90 89
+ 79 80 91 90 80 81 92 91 81 82 93 92
+ 82 83 94 93 83 84 95 94 84 85 96 95
+ 85 86 97 96 86 87 98 97 88 89 100 99
+ 89 90 101 100 90 91 102 101 91 92 103 102
+ 92 93 104 103 93 94 105 104 94 95 106 105
+ 95 96 107 106 96 97 108 107 97 98 109 108
+ 99 100 111 110 100 101 112 111 101 102 113 112
+ 102 103 114 113 103 104 115 114 104 105 116 115
+ 105 106 117 116 106 107 118 117 107 108 119 118
+ 108 109 120 119
+ </DataArray>
+ <DataArray type="Int32" Name="offsets" NumberOfComponents="1" format="ascii">
+ 4 8 12 16 20 24 28 32 36 40 44 48
+ 52 56 60 64 68 72 76 80 84 88 92 96
+ 100 104 108 112 116 120 124 128 132 136 140 144
+ 148 152 156 160 164 168 172 176 180 184 188 192
+ 196 200 204 208 212 216 220 224 228 232 236 240
+ 244 248 252 256 260 264 268 272 276 280 284 288
+ 292 296 300 304 308 312 316 320 324 328 332 336
+ 340 344 348 352 356 360 364 368 372 376 380 384
+ 388 392 396 400
+ </DataArray>
+ <DataArray type="UInt8" Name="types" NumberOfComponents="1" format="ascii">
+ 9 9 9 9 9 9 9 9 9 9 9 9
+ 9 9 9 9 9 9 9 9 9 9 9 9
+ 9 9 9 9 9 9 9 9 9 9 9 9
+ 9 9 9 9 9 9 9 9 9 9 9 9
+ 9 9 9 9 9 9 9 9 9 9 9 9
+ 9 9 9 9 9 9 9 9 9 9 9 9
+ 9 9 9 9 9 9 9 9 9 9 9 9
+ 9 9 9 9 9 9 9 9 9 9 9 9
+ 9 9 9 9
+ </DataArray>
+ </Cells>
+ </Piece>
+ </UnstructuredGrid>
+</VTKFile>
diff --git a/test/references/example_cahn_hilliard_parallel_p0-reference.vtu b/test/references/example_cahn_hilliard_parallel_p0-reference.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..33055a97eeb03d8d6140ce53d1c5286421968576
--- /dev/null
+++ b/test/references/example_cahn_hilliard_parallel_p0-reference.vtu
@@ -0,0 +1,90 @@
+<?xml version="1.0"?>
+<VTKFile type="UnstructuredGrid" version="0.1" byte_order="LittleEndian">
+ <UnstructuredGrid>
+ <Piece NumberOfCells="50" NumberOfPoints="66">
+ <PointData Scalars="c">
+ <DataArray type="Float32" Name="c" NumberOfComponents="1" format="ascii">
+ 0.987413 0.967792 0.987952 0.96968 0.835216 0.846695 0.315888 0.347015 0.07633 0.0842153 0.0310669 0.0318347
+ 0.98934 0.974433 0.876081 0.451628 0.108122 0.0343987 0.990994 0.979831 0.908563 0.599035 0.142578 0.0385368
+ 0.992267 0.983696 0.929883 0.688029 0.17195 0.0424578 0.992715 0.984891 0.935745 0.709774 0.181316 0.0438416
+ 0.992203 0.983164 0.926316 0.673731 0.166303 0.0418624 0.990913 0.978868 0.901396 0.565368 0.133737 0.0379597
+ 0.989316 0.97332 0.867197 0.415015 0.0999852 0.034596 0.988025 0.968704 0.838855 0.324748 0.0786051 0.0330494
+ 0.987533 0.966931 0.828442 0.299419 0.0719689 0.0327484
+ </DataArray>
+ <DataArray type="Float32" Name="mu" NumberOfComponents="1" format="ascii">
+ -0.582482 -0.579776 -0.5887 -0.593969 -0.515998 -0.568131 -0.114131 -0.238969 1.04176 0.96415 2.52288 2.46061
+ -0.59644 -0.613863 -0.646813 -0.49582 0.808299 2.31556 -0.588708 -0.60506 -0.635295 -0.505148 0.717069 2.17269
+ -0.566343 -0.570571 -0.549029 -0.307996 0.726686 2.0862 -0.545819 -0.543905 -0.497521 -0.22 0.746878 2.05582
+ -0.53554 -0.541183 -0.528833 -0.312476 0.723869 2.06566 -0.524652 -0.53964 -0.570403 -0.464335 0.727907 2.11388
+ -0.50005 -0.510371 -0.52179 -0.33466 0.833541 2.19285 -0.468666 -0.462752 -0.404704 -0.0578805 0.982785 2.26844
+ -0.453862 -0.438902 -0.343636 0.0588401 1.05113 2.29946
+ </DataArray>
+ </PointData>
+ <CellData Scalars="process rank">
+ <DataArray type="Float32" Name="process rank" NumberOfComponents="1" format="ascii">
+ 0 0 0 0 0 0 0 0 0 0 0 0
+ 0 0 0 0 0 0 0 0 0 0 0 0
+ 0 0 0 0 0 0 0 0 0 0 0 0
+ 0 0 0 0 0 0 0 0 0 0 0 0
+ 0 0
+ </DataArray>
+ </CellData>
+ <Points>
+ <DataArray type="Float32" Name="Coordinates" NumberOfComponents="3" format="ascii">
+ 0 0 0 0.1 0 0 0 0.1 0 0.1 0.1 0
+ 0.2 0 0 0.2 0.1 0 0.3 0 0 0.3 0.1 0
+ 0.4 0 0 0.4 0.1 0 0.5 0 0 0.5 0.1 0
+ 0 0.2 0 0.1 0.2 0 0.2 0.2 0 0.3 0.2 0
+ 0.4 0.2 0 0.5 0.2 0 0 0.3 0 0.1 0.3 0
+ 0.2 0.3 0 0.3 0.3 0 0.4 0.3 0 0.5 0.3 0
+ 0 0.4 0 0.1 0.4 0 0.2 0.4 0 0.3 0.4 0
+ 0.4 0.4 0 0.5 0.4 0 0 0.5 0 0.1 0.5 0
+ 0.2 0.5 0 0.3 0.5 0 0.4 0.5 0 0.5 0.5 0
+ 0 0.6 0 0.1 0.6 0 0.2 0.6 0 0.3 0.6 0
+ 0.4 0.6 0 0.5 0.6 0 0 0.7 0 0.1 0.7 0
+ 0.2 0.7 0 0.3 0.7 0 0.4 0.7 0 0.5 0.7 0
+ 0 0.8 0 0.1 0.8 0 0.2 0.8 0 0.3 0.8 0
+ 0.4 0.8 0 0.5 0.8 0 0 0.9 0 0.1 0.9 0
+ 0.2 0.9 0 0.3 0.9 0 0.4 0.9 0 0.5 0.9 0
+ 0 1 0 0.1 1 0 0.2 1 0 0.3 1 0
+ 0.4 1 0 0.5 1 0
+ </DataArray>
+ </Points>
+ <Cells>
+ <DataArray type="Int32" Name="connectivity" NumberOfComponents="1" format="ascii">
+ 0 1 3 2 1 4 5 3 4 6 7 5
+ 6 8 9 7 8 10 11 9 2 3 13 12
+ 3 5 14 13 5 7 15 14 7 9 16 15
+ 9 11 17 16 12 13 19 18 13 14 20 19
+ 14 15 21 20 15 16 22 21 16 17 23 22
+ 18 19 25 24 19 20 26 25 20 21 27 26
+ 21 22 28 27 22 23 29 28 24 25 31 30
+ 25 26 32 31 26 27 33 32 27 28 34 33
+ 28 29 35 34 30 31 37 36 31 32 38 37
+ 32 33 39 38 33 34 40 39 34 35 41 40
+ 36 37 43 42 37 38 44 43 38 39 45 44
+ 39 40 46 45 40 41 47 46 42 43 49 48
+ 43 44 50 49 44 45 51 50 45 46 52 51
+ 46 47 53 52 48 49 55 54 49 50 56 55
+ 50 51 57 56 51 52 58 57 52 53 59 58
+ 54 55 61 60 55 56 62 61 56 57 63 62
+ 57 58 64 63 58 59 65 64
+ </DataArray>
+ <DataArray type="Int32" Name="offsets" NumberOfComponents="1" format="ascii">
+ 4 8 12 16 20 24 28 32 36 40 44 48
+ 52 56 60 64 68 72 76 80 84 88 92 96
+ 100 104 108 112 116 120 124 128 132 136 140 144
+ 148 152 156 160 164 168 172 176 180 184 188 192
+ 196 200
+ </DataArray>
+ <DataArray type="UInt8" Name="types" NumberOfComponents="1" format="ascii">
+ 9 9 9 9 9 9 9 9 9 9 9 9
+ 9 9 9 9 9 9 9 9 9 9 9 9
+ 9 9 9 9 9 9 9 9 9 9 9 9
+ 9 9 9 9 9 9 9 9 9 9 9 9
+ 9 9
+ </DataArray>
+ </Cells>
+ </Piece>
+ </UnstructuredGrid>
+</VTKFile>
diff --git a/test/references/example_cahn_hilliard_parallel_p1-reference.vtu b/test/references/example_cahn_hilliard_parallel_p1-reference.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..654a428525ab4430f94642d8aacaa522a2cd99e2
--- /dev/null
+++ b/test/references/example_cahn_hilliard_parallel_p1-reference.vtu
@@ -0,0 +1,90 @@
+<?xml version="1.0"?>
+<VTKFile type="UnstructuredGrid" version="0.1" byte_order="LittleEndian">
+ <UnstructuredGrid>
+ <Piece NumberOfCells="50" NumberOfPoints="66">
+ <PointData Scalars="c">
+ <DataArray type="Float32" Name="c" NumberOfComponents="1" format="ascii">
+ 0.0310669 0.0454055 0.0318347 0.0417179 0.151223 0.122843 0.60494 0.44536 0.913399 0.842182 0.97646 0.9239
+ 0.0343987 0.0342536 0.0743882 0.197663 0.49711 0.669889 0.0385368 0.0284181 0.0439069 0.0838457 0.148188 0.187137
+ 0.0424578 0.0257314 0.0319396 0.0451773 0.059939 0.066897 0.0438416 0.0254117 0.0303727 0.040822 0.0512762 0.0558646
+ 0.0418624 0.0274558 0.0386987 0.0636586 0.0966672 0.114176 0.0379597 0.0335398 0.0679441 0.156054 0.307864 0.41219
+ 0.034596 0.0451329 0.133094 0.438475 0.796547 0.869974 0.0330494 0.0589986 0.221886 0.75329 0.949963 0.987869
+ 0.0327484 0.0656509 0.268929 0.820439 0.980281 1.00762
+ </DataArray>
+ <DataArray type="Float32" Name="mu" NumberOfComponents="1" format="ascii">
+ 2.52288 4.1823 2.46061 4.08921 6.01253 5.8368 7.70006 7.56053 8.05424 8.02932 8.11412 8.13276
+ 2.31556 3.86377 5.42192 6.84719 7.90439 8.06514 2.17269 3.62044 4.99448 6.19814 7.04558 7.35191
+ 2.0862 3.45163 4.71887 5.78889 6.52484 6.78999 2.05582 3.39119 4.64108 5.71237 6.46686 6.74404
+ 2.06566 3.42621 4.73554 5.92798 6.85212 7.22199 2.11388 3.51364 4.91738 6.26625 7.44102 8.01139
+ 2.19285 3.59492 5.03798 6.51027 7.23419 7.48124 2.26844 3.63308 5.01078 6.15979 6.78308 6.98309
+ 2.29946 3.63974 4.95753 5.99225 6.59989 6.79676
+ </DataArray>
+ </PointData>
+ <CellData Scalars="process rank">
+ <DataArray type="Float32" Name="process rank" NumberOfComponents="1" format="ascii">
+ 1 1 1 1 1 1 1 1 1 1 1 1
+ 1 1 1 1 1 1 1 1 1 1 1 1
+ 1 1 1 1 1 1 1 1 1 1 1 1
+ 1 1 1 1 1 1 1 1 1 1 1 1
+ 1 1
+ </DataArray>
+ </CellData>
+ <Points>
+ <DataArray type="Float32" Name="Coordinates" NumberOfComponents="3" format="ascii">
+ 0.5 0 0 0.6 0 0 0.5 0.1 0 0.6 0.1 0
+ 0.7 0 0 0.7 0.1 0 0.8 0 0 0.8 0.1 0
+ 0.9 0 0 0.9 0.1 0 1 0 0 1 0.1 0
+ 0.5 0.2 0 0.6 0.2 0 0.7 0.2 0 0.8 0.2 0
+ 0.9 0.2 0 1 0.2 0 0.5 0.3 0 0.6 0.3 0
+ 0.7 0.3 0 0.8 0.3 0 0.9 0.3 0 1 0.3 0
+ 0.5 0.4 0 0.6 0.4 0 0.7 0.4 0 0.8 0.4 0
+ 0.9 0.4 0 1 0.4 0 0.5 0.5 0 0.6 0.5 0
+ 0.7 0.5 0 0.8 0.5 0 0.9 0.5 0 1 0.5 0
+ 0.5 0.6 0 0.6 0.6 0 0.7 0.6 0 0.8 0.6 0
+ 0.9 0.6 0 1 0.6 0 0.5 0.7 0 0.6 0.7 0
+ 0.7 0.7 0 0.8 0.7 0 0.9 0.7 0 1 0.7 0
+ 0.5 0.8 0 0.6 0.8 0 0.7 0.8 0 0.8 0.8 0
+ 0.9 0.8 0 1 0.8 0 0.5 0.9 0 0.6 0.9 0
+ 0.7 0.9 0 0.8 0.9 0 0.9 0.9 0 1 0.9 0
+ 0.5 1 0 0.6 1 0 0.7 1 0 0.8 1 0
+ 0.9 1 0 1 1 0
+ </DataArray>
+ </Points>
+ <Cells>
+ <DataArray type="Int32" Name="connectivity" NumberOfComponents="1" format="ascii">
+ 0 1 3 2 1 4 5 3 4 6 7 5
+ 6 8 9 7 8 10 11 9 2 3 13 12
+ 3 5 14 13 5 7 15 14 7 9 16 15
+ 9 11 17 16 12 13 19 18 13 14 20 19
+ 14 15 21 20 15 16 22 21 16 17 23 22
+ 18 19 25 24 19 20 26 25 20 21 27 26
+ 21 22 28 27 22 23 29 28 24 25 31 30
+ 25 26 32 31 26 27 33 32 27 28 34 33
+ 28 29 35 34 30 31 37 36 31 32 38 37
+ 32 33 39 38 33 34 40 39 34 35 41 40
+ 36 37 43 42 37 38 44 43 38 39 45 44
+ 39 40 46 45 40 41 47 46 42 43 49 48
+ 43 44 50 49 44 45 51 50 45 46 52 51
+ 46 47 53 52 48 49 55 54 49 50 56 55
+ 50 51 57 56 51 52 58 57 52 53 59 58
+ 54 55 61 60 55 56 62 61 56 57 63 62
+ 57 58 64 63 58 59 65 64
+ </DataArray>
+ <DataArray type="Int32" Name="offsets" NumberOfComponents="1" format="ascii">
+ 4 8 12 16 20 24 28 32 36 40 44 48
+ 52 56 60 64 68 72 76 80 84 88 92 96
+ 100 104 108 112 116 120 124 128 132 136 140 144
+ 148 152 156 160 164 168 172 176 180 184 188 192
+ 196 200
+ </DataArray>
+ <DataArray type="UInt8" Name="types" NumberOfComponents="1" format="ascii">
+ 9 9 9 9 9 9 9 9 9 9 9 9
+ 9 9 9 9 9 9 9 9 9 9 9 9
+ 9 9 9 9 9 9 9 9 9 9 9 9
+ 9 9 9 9 9 9 9 9 9 9 9 9
+ 9 9
+ </DataArray>
+ </Cells>
+ </Piece>
+ </UnstructuredGrid>
+</VTKFile>