diff --git a/examples/README.md b/examples/README.md
index beb935118d8db8a5d64f7bcedeb2a8c74f17f499..057987caa527b6fde04fa9dc369abf8d6814142e 100644
--- a/examples/README.md
+++ b/examples/README.md
@@ -224,3 +224,15 @@ __Discretization method:__ Cell-centered finite volumes with two-point flux appr
<figure><img src="embedded_network_1d3d/img/network.png" alt="blood vessel network"/></figure></td>
</a></td>
</tr></table>
+
+### [:open_file_folder: Example 11: The Cahn-Hilliard model, a nonlinear PDE](cahn_hilliard/README.md)
+
+<table><tr><td>
+
+In this example we simulate the separation of two immiscible phases using the Cahn-Hilliard model.
+
+</td>
+<td width="35%"><a href="cahn_hilliard/README.md">
+<figure><img src="cahn_hilliard/img/animation.gif" alt="phase distribution"/></figure></td>
+</a></td>
+</tr></table>
diff --git a/examples/cahn_hilliard/.doc_config b/examples/cahn_hilliard/.doc_config
new file mode 100644
index 0000000000000000000000000000000000000000..d1a87cc8a3bb170ad8f1db9821477b8b0f5747cd
--- /dev/null
+++ b/examples/cahn_hilliard/.doc_config
@@ -0,0 +1,27 @@
+{
+ "README.md" : [
+ "doc/_intro.md"
+ ],
+
+ "doc/mainfile.md" : [
+ "doc/mainfile_intro.md",
+ "main.cc"
+ ],
+
+ "doc/modelheader.md" : [
+ "doc/modelheader_intro.md",
+ "model.hh"
+ ],
+
+ "navigation" : {
+ "mainpage" : "README.md",
+ "subpages" : [
+ "doc/modelheader.md",
+ "doc/mainfile.md"
+ ],
+ "subtitles" : [
+ "Model implementation",
+ "Main program"
+ ]
+ }
+}
diff --git a/examples/cahn_hilliard/README.md b/examples/cahn_hilliard/README.md
new file mode 100644
index 0000000000000000000000000000000000000000..a192a3346b71a7b07ab85a37530f9ed5383889b1
--- /dev/null
+++ b/examples/cahn_hilliard/README.md
@@ -0,0 +1,70 @@
+<!-- Important: This file has been automatically generated by generate_example_docs.py. Do not edit this file directly! -->
+
+# Cahn-Hilliard model
+
+A random initial distribution of two phases separating according to the Cahn-Hilliard model.
+
+__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_]]
+
+__Result__. The result will look like below.
+You see the concentration variable $c$ capturing the phase distribution after 1 second.
+
+<figure>
+ <center>
+ <img src="img/results_phase_distribution.png" alt="Concentration field modeling phase distribution" width="50%"/>
+ </center>
+</figure>
+
+## Problem set-up
+
+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.
+
+## Model 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.
+
+```math
+\frac{\partial c}{\partial t} = M \nabla^2 \mu \\
+\mu = E (4 c^3 - 6 c^2 +2 c) - \gamma \nabla^2 c
+```
+
+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.
+
+## Implementation of a simple nonlinear PDE
+
+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.
+
+More details are given in [main.cc](doc/mainfile.md) and [model.hh](doc/modelheader.md).
+
+Additionally the folder contains `params.input` containing runtime parameters and the build
+system file `CMakeLists.txt` defining conditions for the accompanying test.
+
+## Part 1: Model implementation
+
+| [:arrow_right: Click to continue with part 1 of the documentation](doc/modelheader.md) |
+|---:|
+
+
+## Part 2: Main program
+
+| [:arrow_right: Click to continue with part 2 of the documentation](doc/mainfile.md) |
+|---:|
diff --git a/examples/cahn_hilliard/doc/_intro.md b/examples/cahn_hilliard/doc/_intro.md
new file mode 100644
index 0000000000000000000000000000000000000000..0c63dab29316bf6ccd1454d672d7b860d40daa61
--- /dev/null
+++ b/examples/cahn_hilliard/doc/_intro.md
@@ -0,0 +1,57 @@
+# Cahn-Hilliard model
+
+A random initial distribution of two phases separating according to the Cahn-Hilliard model.
+
+__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_]]
+
+__Result__. The result will look like below.
+You see the concentration variable $c$ capturing the phase distribution after 1 second.
+
+<figure>
+ <center>
+ <img src="img/results_phase_distribution.png" alt="Concentration field modeling phase distribution" width="50%"/>
+ </center>
+</figure>
+
+## Problem set-up
+
+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.
+
+## Model 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.
+
+```math
+\frac{\partial c}{\partial t} = M \nabla^2 \mu \\
+\mu = E (4 c^3 - 6 c^2 +2 c) - \gamma \nabla^2 c
+```
+
+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.
+
+## Implementation of a simple nonlinear PDE
+
+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.
+
+More details are given in [main.cc](doc/mainfile.md) and [model.hh](doc/modelheader.md).
+
+Additionally the folder contains `params.input` containing runtime parameters and the build
+system file `CMakeLists.txt` defining conditions for the accompanying test.
diff --git a/examples/cahn_hilliard/doc/mainfile.md b/examples/cahn_hilliard/doc/mainfile.md
new file mode 100644
index 0000000000000000000000000000000000000000..ef55ebc6065cc5459f824474642d647914ed1435
--- /dev/null
+++ b/examples/cahn_hilliard/doc/mainfile.md
@@ -0,0 +1,492 @@
+<!-- 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) |
+|---|---:|
+
+
+
+# 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.
+
+
+```cpp
+#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.
+
+
+<details open>
+<summary><b>Click to hide/show the file documentation</b> (or inspect the [source code](../main.cc))</summary>
+
+### Include headers
+
+
+```cpp
+// use the property system and runtime parameters
+#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>
+```
+
+
+### The problem class `CahnHilliardTestProblem`
+In this abstract problem we inherit from the generic `FVProblem`.
+
+
+```cpp
+namespace Dumux {
+
+template<class TypeTag>
+class CahnHilliardTestProblem : public FVProblem<TypeTag>
+{
+```
+
+<details><summary> Click to show alias definitions and local variables</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>
+
+```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");
+ }
+```
+
+
+### 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.
+
+
+```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;
+ }
+```
+
+
+### 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
+ {
+ BoundaryTypes values;
+ values.setAllNeumann();
+ return values;
+ }
+
+ NumEqVector neumannAtPos(const GlobalPosition& globalPos) const
+ { 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`.
+
+
+```cpp
+ Scalar mobility() const
+ { return mobility_; }
+
+ Scalar surfaceTension() const
+ { return surfaceTension_; }
+
+private:
+ Scalar mobility_;
+ 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.
+
+
+<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>; };
+
+// Enable caching
+template<class TypeTag>
+struct EnableGridVolumeVariablesCache<TypeTag, TTag::CahnHilliardTest>
+{ static constexpr bool value = true; };
+
+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; };
+
+} // 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.
+
+
+<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
+struct MinScatter
+{
+ template<class A, class B>
+ static void apply(A& a, const B& b)
+ { a[0] = std::min(a[0], b[0]); }
+};
+
+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;
+ }
+
+ // 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;
+}
+```
+
+
+### The main function
+
+The main function takes the command line arguments, optionally specifying an input file of
+parameters and/or individual key-value pairs of runtime parameters.
+
+
+```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
+
+
+```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
+
+
+```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);
+```
+
+
+### Problem setup
+
+We instantiate also the problem class according to the test properties
+
+
+```cpp
+ // the problem (initial and boundary conditions)
+ using Problem = GetPropType<TypeTag, Properties::Problem>;
+ auto problem = std::make_shared<Problem>(gridGeometry);
+```
+
+
+### 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.
+
+
+```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;
+
+ // the grid variables
+ using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
+ auto gridVariables = std::make_shared<GridVariables>(problem, gridGeometry);
+ gridVariables->init(sol);
+```
+
+
+### Initialize VTK output
+
+
+```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
+
+
+```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);
+```
+
+
+### Assembler, linear and nonlinear solver
+
+
+```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);
+
+ // the linear solver
+ using LinearSolver = SSORBiCGSTABIstlSolver<
+ LinearSolverTraits<GridGeometry>,
+ LinearAlgebraTraitsFromAssembler<Assembler>
+ >;
+ auto linearSolver = std::make_shared<LinearSolver>(gridGeometry->gridView(), gridGeometry->dofMapper());
+
+ // the solver
+ using Solver = Dumux::NewtonSolver<Assembler, LinearSolver>;
+ Solver solver(assembler, linearSolver);
+```
+
+
+### Time loop
+
+
+```cpp
+ // time loop
+ timeLoop->start(); do
+ {
+ // assemble & solve
+ 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
+ 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());
+```
+
+
+### Finalize
+
+
+```cpp
+ timeLoop->finalize(leafGridView.comm());
+
+ return 0;
+}
+```
+
+
+</details>
+
+
+| [:arrow_left: Back to the main documentation](../README.md) | [:arrow_left: Go back to part 1](modelheader.md) |
+|---|---:|
+
diff --git a/examples/cahn_hilliard/doc/mainfile_intro.md b/examples/cahn_hilliard/doc/mainfile_intro.md
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/examples/cahn_hilliard/doc/modelheader.md b/examples/cahn_hilliard/doc/modelheader.md
new file mode 100644
index 0000000000000000000000000000000000000000..125ea472c8b868a4e264198cddf2d1234aac2680
--- /dev/null
+++ b/examples/cahn_hilliard/doc/modelheader.md
@@ -0,0 +1,394 @@
+<!-- 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) |
+|---|---:|
+
+
+
+# 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.
+
+## VolumeVariables
+
+The volume variables store the local element volume variables, both primary and secondary.
+
+
+<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 {
+
+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 recomputes secondary variables.
+
+
+```cpp
+ /*!
+ * \brief 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
+
+
+```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_; }
+```
+
+
+### 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
+
+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.
+
+
+<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 model properties, depending on the discretization scheme.
+
+
+```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.
+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
+ {
+ NumEqVector storage;
+ storage[Indices::massBalanceEqIdx] = volVars.concentration();
+ storage[Indices::chemicalPotentialEqIdx] = 0.0;
+ return storage;
+ }
+```
+
+
+### The flux terms
+`computeFlux` computes the fluxes over a subcontrolvolumeface, including the integration over
+the area of the face.
+
+
+```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,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf,
+ const ElementFluxVariablesCache& elemFluxVarsCache) const
+ {
+ 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()));
+ }
+
+ 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();
+ return flux;
+ }
+```
+
+
+### Source terms
+
+`computeSource` defines the sources terms at 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.
+
+
+```cpp
+ /*!
+ * \brief Calculate the source term of the equation
+ */
+ NumEqVector computeSource(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolume &scv) const
+ {
+ NumEqVector source(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
+
+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.
+
+
+<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
+
+
+```cpp
+namespace TTag {
+struct CahnHilliardModel {};
+} // end namespace TTag
+```
+
+
+### 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.
+
+
+```cpp
+//! Set the default type of scalar values to double
+template<class TypeTag>
+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.
+
+
+```cpp
+//! Set the model traits property
+template<class TypeTag>
+struct ModelTraits<TypeTag, TTag::CahnHilliardModel>
+{
+ struct Traits
+ {
+ 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; }
+ };
+
+ using type = Traits;
+};
+```
+
+
+### Further model properties
+
+Define further properties of the model, selecting the local residual and volumevariables defined
+above.
+
+
+```cpp
+template<class TypeTag>
+struct LocalResidual<TypeTag, TTag::CahnHilliardModel>
+{ using type = CahnHilliardModelLocalResidual<TypeTag>; };
+
+//! Set the volume variables property
+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](mainfile.md) |
+|---|---:|
+
diff --git a/examples/cahn_hilliard/doc/modelheader_intro.md b/examples/cahn_hilliard/doc/modelheader_intro.md
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
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
index fa1bc08904d39562adfd5e244988c5ce95a3916c..4e18dcff990e0bf4b50dbc5fa9cdc45988868260 100644
--- a/examples/cahn_hilliard/main.cc
+++ b/examples/cahn_hilliard/main.cc
@@ -17,26 +17,43 @@
* 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.
+//
#include <config.h>
-
-//////////////////////////////////////////////////////////////////////////////
-//////////////////////////////////////////////////////////////////////////////
-// Problem //////// (often in a separate file problem.hh) ////////////////////
-//////////////////////////////////////////////////////////////////////////////
-//////////////////////////////////////////////////////////////////////////////
-
+// ## 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.
+//
+// [[content]]
+// ### Include headers
+//
+// [[codeblock]]
+// use the property system and runtime parameters
#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]]
+//
+// ### The problem class `CahnHilliardTestProblem`
+// In this abstract problem we inherit from the generic `FVProblem`.
+//
+// [[codeblock]]
namespace Dumux {
template<class TypeTag>
class CahnHilliardTestProblem : public FVProblem<TypeTag>
{
+// [[/codeblock]]
+// [[details]] alias definitions and local variables
+// [[codeblock]]
using ParentType = FVProblem<TypeTag>;
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
using FVElementGeometry = typename GridGeometry::LocalView;
@@ -50,6 +67,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]]
public:
CahnHilliardTestProblem(std::shared_ptr<const GridGeometry> gridGeometry)
: ParentType(gridGeometry)
@@ -58,8 +78,15 @@ public:
surfaceTension_ = getParam<Scalar>("Problem.SurfaceTension");
energyScale_ = getParam<Scalar>("Problem.EnergyScale");
}
-
- // here we implement the derivative of the free energy
+ // [[/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]]
template<class ElementVolumeVariables>
NumEqVector source(const Element &element,
const FVElementGeometry& fvGeometry,
@@ -71,7 +98,14 @@ 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
{
BoundaryTypes values;
@@ -81,7 +115,13 @@ 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`.
+ //
+ // [[codeblock]]
Scalar mobility() const
{ return mobility_; }
@@ -95,29 +135,52 @@ private:
};
} // end namespace Dumux
-
-//////////////////////////////////////////////////////////////////////////////
-
-//////////////////////////////////////////////////////////////////////////////
-//////////////////////////////////////////////////////////////////////////////
-// Test case properties/traits /// (often in a separate file properties.hh) //
-//////////////////////////////////////////////////////////////////////////////
-//////////////////////////////////////////////////////////////////////////////
-
+// [[/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
+//
+// We specify a grid to be used in the test, select our problem class and enable caching.
+//
+// [[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>; };
@@ -136,24 +199,52 @@ struct EnableGridGeometryCache<TypeTag, TTag::CahnHilliardTest>
{ static constexpr bool value = true; };
} // 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
+//
+// 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)
+//
+// [[codeblock]]
struct MinScatter
{
template<class A, class B>
@@ -161,19 +252,70 @@ struct MinScatter
{ a[0] = std::min(a[0], b[0]); }
};
+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;
+ }
+
+ // 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]]
+//
+// ### The main function
+//
+// The main function takes the command line arguments, optionally specifying an input file of
+// parameters and/or individual key-value pairs of runtime parameters.
+//
+// [[codeblock]]
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
Dumux::initialize(argc, argv);
// 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();
@@ -184,37 +326,27 @@ int main(int argc, char** argv)
// create the finite volume grid geometry
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>;
- SolutionVector sol(gridGeometry->numDofs());
-
- // create a random initial field
- 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)) + leafGridView.comm().rank();
- sol[n][1] = 0.0;
- }
-
- // take the value of the processor with the minimum rank
- VectorCommDataHandle<
- typename GridGeometry::VertexMapper,
- SolutionVector,
- GridGeometry::GridView::dimension,
- MinScatter
- > minHandle(gridGeometry->vertexMapper(), sol);
- leafGridView.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]);
-
+ auto sol = createInitialSolution<SolutionVector>(*gridGeometry);
// copy the vector to store state of previous time step
auto oldSol = sol;
@@ -222,13 +354,21 @@ int main(int argc, char** argv)
using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
auto gridVariables = std::make_shared<GridVariables>(problem, gridGeometry);
gridVariables->init(sol);
-
+ // [[/codeblock]]
+ //
+ // ### Initialize VTK output
+ //
+ // [[codeblock]]
// 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);
-
+ // [[/codeblock]]
+ //
+ // ### Set up time loop
+ //
+ // [[codeblock]]
// get some time loop parameters
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
@@ -238,7 +378,11 @@ int main(int argc, char** argv)
// 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);
@@ -253,7 +397,11 @@ int main(int argc, char** argv)
// the solver
using Solver = Dumux::NewtonSolver<Assembler, LinearSolver>;
Solver solver(assembler, linearSolver);
-
+ // [[/codeblock]]
+ //
+ // ### Time loop
+ //
+ // [[codeblock]]
// time loop
timeLoop->start(); do
{
@@ -277,8 +425,14 @@ int main(int argc, char** argv)
timeLoop->setTimeStepSize(solver.suggestTimeStepSize(timeLoop->timeStepSize()));
} while (!timeLoop->finished());
-
+ // [[/codeblock]]
+ //
+ // ### Finalize
+ //
+ // [[codeblock]]
timeLoop->finalize(leafGridView.comm());
return 0;
}
+// [[/codeblock]]
+// [[/content]]
diff --git a/examples/cahn_hilliard/model.hh b/examples/cahn_hilliard/model.hh
index 51b5eb62be61a8a886a9696ee5378ddad6ef2df8..9f69f9dadabc0a18945fcdfcddf432aaa6c90272 100644
--- a/examples/cahn_hilliard/model.hh
+++ b/examples/cahn_hilliard/model.hh
@@ -20,12 +20,19 @@
#ifndef DUMUX_EXAMPLES_CAHN_HILLIARD_MODEL_HH
#define DUMUX_EXAMPLES_CAHN_HILLIARD_MODEL_HH
-//////////////////////////////////////////////////////////////////////////////
-//////////////////////////////////////////////////////////////////////////////
-// VolumeVariables //////// (often in a separate file volumevariables.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.
+//
+// ## VolumeVariables
+//
+// The volume variables store the local element volume variables, both primary and secondary.
+//
+// [[content]]
+//
+// ### The `CahnHilliardModelVolumeVariables` class
+//
+// [[codeblock]]
namespace Dumux {
template <class Traits>
@@ -38,7 +45,14 @@ public:
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 recomputes secondary variables.
+ //
+ // [[codeblock]]
/*!
* \brief Update all quantities for a given control volume
*/
@@ -50,7 +64,13 @@ public:
{
priVars_ = elemSol[scv.indexInElement()];
}
-
+ // [[/codeblock]]
+ //
+ // ### Access functions
+ //
+ // Named and generic functions to access different primary variables
+ //
+ // [[codeblock]]
Scalar concentration() const
{ return priVars_[Indices::concentrationIdx]; }
@@ -62,29 +82,52 @@ 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
-
-//////////////////////////////////////////////////////////////////////////////
-
-//////////////////////////////////////////////////////////////////////////////
-//////////////////////////////////////////////////////////////////////////////
-// LocalResidual //////////// (often in a separate file localresidual.hh) ////
-//////////////////////////////////////////////////////////////////////////////
-//////////////////////////////////////////////////////////////////////////////
-
+// [[/codeblock]]
+// [[/content]]
+//
+// ## LocalResidual
+//
+// 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.
+//
+// [[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 model properties, depending on the discretization scheme.
+//
+// [[codeblock]]
namespace Dumux {
template<class TypeTag>
@@ -93,6 +136,10 @@ class CahnHilliardModelLocalResidual
{
// 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>;
@@ -113,7 +160,16 @@ class CahnHilliardModelLocalResidual
static constexpr int dimWorld = GridView::dimensionworld;
public:
using ParentType::ParentType;
-
+ // [[/codeblock]]
+ // [[/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.
+ // 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
*/
@@ -126,7 +182,13 @@ public:
storage[Indices::chemicalPotentialEqIdx] = 0.0;
return storage;
}
-
+ // [[/codeblock]]
+ //
+ // ### The flux terms
+ // `computeFlux` computes the fluxes over a subcontrolvolumeface, including the integration over
+ // the area of the face.
+ //
+ // [[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
@@ -159,7 +221,17 @@ public:
flux[Indices::chemicalPotentialEqIdx] = -1.0*vtmv(scvf.unitOuterNormal(), gamma, gradConcentration)*scvf.area();
return flux;
}
-
+ // [[/codeblock]]
+ //
+ // ### Source terms
+ //
+ // `computeSource` defines the sources terms at 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.
+ //
+ // [[codeblock]]
/*!
* \brief Calculate the source term of the equation
*/
@@ -181,22 +253,40 @@ public:
};
} // end namespace Dumux
-//////////////////////////////////////////////////////////////////////////////
-
-//////////////////////////////////////////////////////////////////////////////
-//////////////////////////////////////////////////////////////////////////////
-// Model properties/traits ///////////////////////////////////////////////////
-//////////////////////////////////////////////////////////////////////////////
-//////////////////////////////////////////////////////////////////////////////
-
+// [[/codeblock]]
+// [[/content]]
+//
+// ## 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.
+//
+// [[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
template<class TypeTag>
struct Scalar<TypeTag, TTag:: CahnHilliardModel >
@@ -211,7 +301,14 @@ struct PrimaryVariables<TypeTag, TTag:: CahnHilliardModel >
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
template<class TypeTag>
struct ModelTraits<TypeTag, TTag::CahnHilliardModel>
@@ -232,7 +329,14 @@ struct ModelTraits<TypeTag, TTag::CahnHilliardModel>
using type = Traits;
};
-
+// [[/codeblock]]
+//
+// ### Further model properties
+//
+// Define further properties of the model, selecting the local residual and volumevariables defined
+// above.
+//
+// [[codeblock]]
template<class TypeTag>
struct LocalResidual<TypeTag, TTag::CahnHilliardModel>
{ using type = CahnHilliardModelLocalResidual<TypeTag>; };
@@ -251,5 +355,7 @@ struct VolumeVariables<TypeTag, TTag::CahnHilliardModel>
};
} // end namespace Dumux::Properties
+// [[/codeblock]]
+// [[/content]]
#endif