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/cahn_hilliard/CMakeLists.txt b/examples/cahn_hilliard/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..580fe4111efe2719777fc6270fd383d3d37a22ee
--- /dev/null
+++ b/examples/cahn_hilliard/CMakeLists.txt
@@ -0,0 +1,2 @@
+dune_symlink_to_source_files(FILES "params.input")
+dumux_add_test(NAME example_cahn_hilliard SOURCES main.cc)
diff --git a/examples/cahn_hilliard/main.cc b/examples/cahn_hilliard/main.cc
new file mode 100644
index 0000000000000000000000000000000000000000..fa1bc08904d39562adfd5e244988c5ce95a3916c
--- /dev/null
+++ b/examples/cahn_hilliard/main.cc
@@ -0,0 +1,284 @@
+// -*- 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/>. *
+ *****************************************************************************/
+
+#include <config.h>
+
+//////////////////////////////////////////////////////////////////////////////
+//////////////////////////////////////////////////////////////////////////////
+// Problem //////// (often in a separate file problem.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>
+
+namespace Dumux {
+
+template<class TypeTag>
+class CahnHilliardTestProblem : public FVProblem<TypeTag>
+{
+ 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;
+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");
+ }
+
+ // here we implement the derivative of the free energy
+ 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;
+ }
+
+ BoundaryTypes boundaryTypesAtPos(const GlobalPosition& globalPos) const
+ {
+ BoundaryTypes values;
+ values.setAllNeumann();
+ return values;
+ }
+
+ NumEqVector neumannAtPos(const GlobalPosition& globalPos) const
+ { return { 0.0, 0.0 }; }
+
+ Scalar mobility() const
+ { return mobility_; }
+
+ Scalar surfaceTension() const
+ { return surfaceTension_; }
+
+private:
+ Scalar mobility_;
+ Scalar surfaceTension_;
+ Scalar energyScale_;
+};
+
+} // end namespace Dumux
+
+//////////////////////////////////////////////////////////////////////////////
+
+//////////////////////////////////////////////////////////////////////////////
+//////////////////////////////////////////////////////////////////////////////
+// Test case properties/traits /// (often in a separate file properties.hh) //
+//////////////////////////////////////////////////////////////////////////////
+//////////////////////////////////////////////////////////////////////////////
+
+#include <dune/grid/yaspgrid.hh>
+#include <dumux/discretization/box.hh>
+#include "model.hh"
+
+namespace Dumux::Properties {
+
+namespace TTag {
+struct CahnHilliardTest { using InheritsFrom = std::tuple<CahnHilliardModel, BoxModel>; };
+} // end namespace TTag
+
+template<class TypeTag>
+struct Grid<TypeTag, TTag::CahnHilliardTest>
+{ using type = Dune::YaspGrid<2>; };
+
+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
+
+//////////////////////////////////////////////////////////////////////////////
+
+#include <random>
+
+#include <dumux/common/initialize.hh>
+#include <dumux/common/properties.hh>
+#include <dumux/common/parameters.hh>
+
+#include <dumux/io/vtkoutputmodule.hh>
+#include <dumux/io/grid/gridmanager_yasp.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>
+
+struct MinScatter
+{
+ template<class A, class B>
+ static void apply(A& a, const B& b)
+ { a[0] = std::min(a[0], b[0]); }
+};
+
+int main(int argc, char** argv)
+{
+ using namespace Dumux;
+
+ // define the type tag for this problem
+ using TypeTag = Properties::TTag::CahnHilliardTest;
+
+ // maybe initialize MPI and/or multithreading backend
+ Dumux::initialize(argc, argv);
+
+ // initialize parameter tree
+ Parameters::init(argc, argv);
+
+ // 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);
+
+ // the problem (initial and boundary conditions)
+ using Problem = GetPropType<TypeTag, Properties::Problem>;
+ auto problem = std::make_shared<Problem>(gridGeometry);
+
+ // 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]);
+
+ // 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 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);
+
+ // 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);
+
+ // 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
+ 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());
+
+ timeLoop->finalize(leafGridView.comm());
+
+ return 0;
+}
diff --git a/examples/cahn_hilliard/model.hh b/examples/cahn_hilliard/model.hh
new file mode 100644
index 0000000000000000000000000000000000000000..51b5eb62be61a8a886a9696ee5378ddad6ef2df8
--- /dev/null
+++ b/examples/cahn_hilliard/model.hh
@@ -0,0 +1,255 @@
+// -*- 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
+
+//////////////////////////////////////////////////////////////////////////////
+//////////////////////////////////////////////////////////////////////////////
+// VolumeVariables //////// (often in a separate file volumevariables.hh) ////
+//////////////////////////////////////////////////////////////////////////////
+//////////////////////////////////////////////////////////////////////////////
+
+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;
+
+ /*!
+ * \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()];
+ }
+
+ 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_; }
+
+ // for compatibility with more general models
+ Scalar extrusionFactor() const
+ { return 1.0; }
+
+private:
+ PrimaryVariables priVars_;
+};
+
+} // end namespace Dumux
+
+//////////////////////////////////////////////////////////////////////////////
+
+//////////////////////////////////////////////////////////////////////////////
+//////////////////////////////////////////////////////////////////////////////
+// LocalResidual //////////// (often in a separate file localresidual.hh) ////
+//////////////////////////////////////////////////////////////////////////////
+//////////////////////////////////////////////////////////////////////////////
+
+#include <dumux/common/math.hh>
+#include <dumux/common/properties.hh>
+#include <dumux/common/numeqvector.hh>
+
+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>;
+ 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;
+
+ /*!
+ * \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;
+ }
+
+ /*!
+ * \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;
+ }
+
+ /*!
+ * \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
+//////////////////////////////////////////////////////////////////////////////
+
+//////////////////////////////////////////////////////////////////////////////
+//////////////////////////////////////////////////////////////////////////////
+// Model properties/traits ///////////////////////////////////////////////////
+//////////////////////////////////////////////////////////////////////////////
+//////////////////////////////////////////////////////////////////////////////
+
+#include <dumux/common/properties.hh>
+
+namespace Dumux::Properties {
+
+namespace TTag {
+struct CahnHilliardModel {};
+} // end namespace TTag
+
+//! 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()
+ >;
+};
+
+//! 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;
+};
+
+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
+
+#endif
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