diff --git a/examples/2pinfiltration/.doc_config b/examples/2pinfiltration/.doc_config
index a8f7fb29c68864f042611612790f533869144d5e..2329831e93837e87c3bfb8b7c84dfbf47cdfa06d 100644
--- a/examples/2pinfiltration/.doc_config
+++ b/examples/2pinfiltration/.doc_config
@@ -1,9 +1,29 @@
{
"README.md" : [
- "doc/intro.md",
- "spatialparams.hh",
- "problem.hh",
+ "doc/_intro.md"
+ ],
+
+ "doc/2p.md" : [
+ "doc/2p_intro.md",
"properties.hh",
+ "problem.hh",
+ "spatialparams.hh"
+ ],
+
+ "doc/main.md" : [
+ "doc/main_intro.md",
"main.cc"
- ]
+ ],
+
+ "navigation" : {
+ "mainpage" : "README.md",
+ "subpages" : [
+ "doc/2p.md",
+ "doc/main.md"
+ ],
+ "subtitles" : [
+ "Two-phase infiltration set-up",
+ "Main program flow"
+ ]
+ }
}
diff --git a/examples/2pinfiltration/README.md b/examples/2pinfiltration/README.md
index 548cc6268d65cc8d98dba0bec26d911485b91d1a..371f6944d092b6fcf961ad81fedfa8751b0ab48d 100644
--- a/examples/2pinfiltration/README.md
+++ b/examples/2pinfiltration/README.md
@@ -6,10 +6,11 @@ __In this example, you will learn how to__
* solve a two-phase flow in porous media problem with two immiscible phases
* set boundary conditions and a simple injection well
-* specify a lens with different porous material parameters
-* use adaptive grid refinement around the saturation front
* specify a point source
* read the initial solution from a text file
+* specify a lens with different porous material parameters
+* use adaptive grid refinement around the saturation front
+
__Prerequisites__. You need [dune-alugrid](https://gitlab.dune-project.org/extensions/dune-alugrid) in order to compile and run this example.
@@ -71,1048 +72,15 @@ For more information about the discretization please have a look at the [handboo
├── spatialparams.hh -> spatial parameter fields
└── initialsolutioncc.txt -> text file with initial solution
```
+In the documentations behind the links provided in the following, you will find a detailed description how the the above mentioned set-up is realized in this example.
+## Part 1: Two-phase infiltration set-up
-## The file `spatialparams.hh`
-
-<details open>
-<summary><b>Click to hide/show the file documentation</b> (or inspect the [source code](spatialparams.hh))</summary>
-
-
-### Includes
-we include the basic spatial parameters for finite volumes file from which we will inherit
-
-```cpp
-#include <dumux/material/spatialparams/fv.hh>
-```
-
-we include all laws which are needed to define the interaction between the solid matrix and the fluids, e.g. laws for capillary pressure saturation relationships.
-
-```cpp
-#include <dumux/material/fluidmatrixinteractions/2p/regularizedvangenuchten.hh>
-#include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh>
-
-namespace Dumux {
-```
-
-In the TwoPTestSpatialParams class we define all functions needed to describe the porous matrix, e.g. porosity and permeability
-
-```cpp
-template<class GridGeometry, class Scalar>
-class TwoPTestSpatialParams
-: public FVSpatialParams<GridGeometry, Scalar, TwoPTestSpatialParams<GridGeometry, Scalar>>
-{
-```
-
-we introduce using declarations that are derived from the property system which we need in this class
-
-```cpp
- using GridView = typename GridGeometry::GridView;
- using Element = typename GridView::template Codim<0>::Entity;
- using FVElementGeometry = typename GridGeometry::LocalView;
- using SubControlVolume = typename FVElementGeometry::SubControlVolume;
- using ThisType = TwoPTestSpatialParams<GridGeometry, Scalar>;
- using ParentType = FVSpatialParams<GridGeometry, Scalar, ThisType>;
-
- static constexpr int dimWorld = GridView::dimensionworld;
- using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
-
- using EffectiveLaw = RegularizedVanGenuchten<Scalar>;
-
-public:
- using MaterialLaw = EffToAbsLaw<EffectiveLaw>;
- using MaterialLawParams = typename MaterialLaw::Params;
- using PermeabilityType = Scalar;
-
- TwoPTestSpatialParams(std::shared_ptr<const GridGeometry> gridGeometry)
- : ParentType(gridGeometry)
- {
-```
-
-we get the position of the lens from the params.input file. The lens is defined by the position of the lower left and the upper right corner
-
-```cpp
- lensLowerLeft_ = getParam<GlobalPosition>("SpatialParams.LensLowerLeft");
- lensUpperRight_ = getParam<GlobalPosition>("SpatialParams.LensUpperRight");
-```
-
-we set the parameters for the material law (here Van-Genuchten Law). First we set the residual saturations for the wetting phase and the non-wetting phase. lensMaterialParams_ define the material parameters for the lens while outerMaterialParams_ define material params for the rest of the domain.
-
-```cpp
- lensMaterialParams_.setSwr(0.18);
- lensMaterialParams_.setSnr(0.0);
- outerMaterialParams_.setSwr(0.05);
- outerMaterialParams_.setSnr(0.0);
-```
-
-we set the parameters for the Van Genuchten law alpha and n
-
-```cpp
- lensMaterialParams_.setVgAlpha(0.00045);
- lensMaterialParams_.setVgn(7.3);
- outerMaterialParams_.setVgAlpha(0.0037);
- outerMaterialParams_.setVgn(4.7);
-```
-
-here we get the permeabilities from the params.input file. In case that no parameter is set, the default parameters (9.05e-12 and 4.6e-10) are used
-
-```cpp
- lensK_ = getParam<Scalar>("SpatialParams.lensK", 9.05e-12);
- outerK_ = getParam<Scalar>("SpatialParams.outerK", 4.6e-10);
- }
-```
-
-We define the (intrinsic) permeability $`[m^2]`$. In this test, we use element-wise distributed permeabilities.
-
-```cpp
- template<class ElementSolution>
- PermeabilityType permeability(const Element& element,
- const SubControlVolume& scv,
- const ElementSolution& elemSol) const
- {
- if (isInLens_(element.geometry().center()))
- return lensK_;
- return outerK_;
- }
-```
-
-We set the porosity $`[-]`$ depending on the position
-
-```cpp
- Scalar porosityAtPos(const GlobalPosition& globalPos) const
- {
- if (isInLens_(globalPos))
- return 0.2;
- return 0.4;
- }
-```
-
-We set the parameter object for the Van Genuchten material law.
-
-```cpp
- template<class ElementSolution>
- const MaterialLawParams& materialLawParams(const Element& element,
- const SubControlVolume& scv,
- const ElementSolution& elemSol) const
- {
- if (isInLens_(element.geometry().center()))
- return lensMaterialParams_;
- return outerMaterialParams_;
- }
-```
-
-Here we can define which phase is to be considered as the wetting phase. Here the wetting phase is the the first phase of the fluidsystem. In this case that is water.
-
-```cpp
- template<class FluidSystem>
- int wettingPhaseAtPos(const GlobalPosition& globalPos) const
- { return FluidSystem::phase0Idx; }
-
-private:
-```
-
-we have a convenience definition of the position of the lens
-
-```cpp
- bool isInLens_(const GlobalPosition &globalPos) const
- {
- for (int i = 0; i < dimWorld; ++i) {
- if (globalPos[i] < lensLowerLeft_[i] + eps_ || globalPos[i] > lensUpperRight_[i] - eps_)
- return false;
- }
- return true;
- }
-
- GlobalPosition lensLowerLeft_;
- GlobalPosition lensUpperRight_;
-
- Scalar lensK_;
- Scalar outerK_;
- MaterialLawParams lensMaterialParams_;
- MaterialLawParams outerMaterialParams_;
-
- static constexpr Scalar eps_ = 1.5e-7;
-};
-
-} // end namespace Dumux
-```
-
-
-</details>
-
-
-
-## The file `problem.hh`
-
-<details open>
-<summary><b>Click to hide/show the file documentation</b> (or inspect the [source code](problem.hh))</summary>
-
-
-### Includes
-We start with includes for `PorousMediumFlowProblem` and `readFileToContainer` (used below).
-
-```cpp
-#include <dumux/porousmediumflow/problem.hh>
-#include <dumux/io/container.hh>
-```
-
-### Problem class
-The problem class `PointSourceProblem` implements boundary and initial conditions.
-It derives from the `PorousMediumFlowProblem` class.
-
-```cpp
-namespace Dumux {
-
-template <class TypeTag>
-class PointSourceProblem : public PorousMediumFlowProblem<TypeTag>
-{
-```
-
-The class implementation starts with some alias declarations and index definitions for convenience
-<details><summary>Click to show local alias declarations and indices</summary>
-
-```cpp
- using ParentType = PorousMediumFlowProblem<TypeTag>;
- using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
- using Element = typename GridView::template Codim<0>::Entity;
- using Vertex = typename GridView::template Codim<GridView::dimensionworld>::Entity;
- using Scalar = GetPropType<TypeTag, Properties::Scalar>;
- using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
- using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
- using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
- using PointSource = GetPropType<TypeTag, Properties::PointSource>;
- using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
- using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
- using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
- using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
-
- enum {
- pressureH2OIdx = Indices::pressureIdx,
- saturationDNAPLIdx = Indices::saturationIdx,
- contiDNAPLEqIdx = Indices::conti0EqIdx + FluidSystem::comp1Idx,
- waterPhaseIdx = FluidSystem::phase0Idx,
- dnaplPhaseIdx = FluidSystem::phase1Idx
- };
-```
-
-</details>
-
-In the constructor of the class, we call the parent type's constructor
-and read the intial values for the primary variables from a text file.
-The function `readFileToContainer` is implemented in the header `dumux/io/container.hh`.
-
-```cpp
-public:
- PointSourceProblem(std::shared_ptr<const GridGeometry> gridGeometry)
- : ParentType(gridGeometry)
- {
- initialValues_ = readFileToContainer<std::vector<PrimaryVariables>>("initialsolutioncc.txt");
- }
-```
-
-For isothermal problems, Dumux requires problem classes to implement a `temperature()`
-member function. Fluid properties that depend on temperature will be calculated with the specified temperature.
-
-```cpp
- Scalar temperature() const
- {
- return 293.15; // 10°C
- }
-```
-
-### 1. Boundary types
-We define the type of boundary conditions depending on location. Two types of boundary conditions
-can be specified: Dirichlet or Neumann boundary condition. On a Dirichlet boundary, the values of the
-primary variables need to be fixed. On a Neumann boundary condition, values for derivatives need to be fixed.
-Mixed boundary conditions (different types for different equations on the same boundary) are not accepted.
-
-```cpp
- BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
- {
- BoundaryTypes values;
- // Dirichlet boundaries on the left and right hand side of the domain
- if (onLeftBoundary_(globalPos) || onRightBoundary_(globalPos))
- values.setAllDirichlet();
- // and Neumann boundaries otherwise (top and bottom of the domain)
- else
- values.setAllNeumann();
- return values;
- }
-```
-
-### 2. Dirichlet boundaries
-We specify the values for the Dirichlet boundaries, depending on location.
-We need to fix values for the two primary variables: the water pressure
-and the DNAPL saturation.
-
-```cpp
- PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
- {
- // To determine the density of water for a given state, we build a fluid state with the given conditions:
- PrimaryVariables values;
- GetPropType<TypeTag, Properties::FluidState> fluidState;
- fluidState.setTemperature(temperature());
- fluidState.setPressure(waterPhaseIdx, /*pressure=*/1e5);
- fluidState.setPressure(dnaplPhaseIdx, /*pressure=*/1e5);
-
- // The density is then calculated by the fluid system:
- Scalar densityW = FluidSystem::density(fluidState, waterPhaseIdx);
-
- // The water phase pressure is the hydrostatic pressure, scaled with a factor:
- Scalar height = this->gridGeometry().bBoxMax()[1] - this->gridGeometry().bBoxMin()[1];
- Scalar depth = this->gridGeometry().bBoxMax()[1] - globalPos[1];
- Scalar alpha = 1 + 1.5/height;
- Scalar width = this->gridGeometry().bBoxMax()[0] - this->gridGeometry().bBoxMin()[0];
- Scalar factor = (width*alpha + (1.0 - alpha)*globalPos[0])/width;
-
- values[pressureH2OIdx] = 1e5 - factor*densityW*this->spatialParams().gravity(globalPos)[1]*depth;
- // The saturation of the DNAPL Trichlorethene is zero on our Dirichlet boundary:
- values[saturationDNAPLIdx] = 0.0;
-
- return values;
- }
-```
-
-### 3. Neumann boundaries
-In our case, we need to specify mass fluxes for our two liquid phases.
-Negative sign means influx and the unit of the boundary flux is $`kg/(m^2 s)`$.
-On the inlet area, we set a DNAPL influx of $`0.04 kg/(m^2 s)`$. On all other
-Neumann boundaries, the boundary flux is zero.
-
-```cpp
- NumEqVector neumannAtPos(const GlobalPosition &globalPos) const
- {
- NumEqVector values(0.0);
- if (onInlet_(globalPos))
- values[contiDNAPLEqIdx] = -0.04;
-
- return values;
- }
-```
-
-### 4. Initial conditions
-The initial condition needs to be set for all primary variables.
-Here, we take the data from the file that we read in previously.
-
-```cpp
- PrimaryVariables initial(const Element& element) const
- {
- // The input data is written for a uniform grid with discretization length delta.
- // Accordingly, we need to find the index of our cells, depending on the x and y coordinates,
- // that corresponds to the indices of the input data set.
- const auto delta = 0.0625;
- unsigned int cellsX = this->gridGeometry().bBoxMax()[0]/delta;
- const auto globalPos = element.geometry().center();
-
- unsigned int dataIdx = std::trunc(globalPos[1]/delta) * cellsX + std::trunc(globalPos[0]/delta);
- return initialValues_[dataIdx];
- }
-```
-
-### 5. Point source
-In this scenario, we set a point source (e.g. modeling a well). The point source value can be solution dependent.
-Point sources are added by pushing them into the vector `pointSources`.
-The `PointSource` constructor takes two arguments.
-The first argument is a coordinate array containing the position in space,
-the second argument is an array of source value for each equation (in units of $`kg/s`$).
-Recall that the first eqution is the water phase mass balance
-and the second equation is the DNAPL phase mass balance.
-
-```cpp
- void addPointSources(std::vector<PointSource>& pointSources) const
- {
- pointSources.push_back(PointSource({0.502, 3.02}, {0, 0.1}));
- }
-```
-
-In the private part of the class, we define some helper functions for
-the boundary conditions and local variables.
-<details><summary>Click to show private data members and functions</summary>
-
-```cpp
-private:
- bool onLeftBoundary_(const GlobalPosition &globalPos) const
- { return globalPos[0] < this->gridGeometry().bBoxMin()[0] + eps_; }
-
- bool onRightBoundary_(const GlobalPosition &globalPos) const
- { return globalPos[0] > this->gridGeometry().bBoxMax()[0] - eps_; }
-
- bool onUpperBoundary_(const GlobalPosition &globalPos) const
- { return globalPos[1] > this->gridGeometry().bBoxMax()[1] - eps_; }
-
- bool onInlet_(const GlobalPosition &globalPos) const
- {
- Scalar width = this->gridGeometry().bBoxMax()[0] - this->gridGeometry().bBoxMin()[0];
- Scalar lambda = (this->gridGeometry().bBoxMax()[0] - globalPos[0])/width;
- return onUpperBoundary_(globalPos) && 0.5 < lambda && lambda < 2.0/3.0;
- }
-
- static constexpr Scalar eps_ = 1e-6;
- std::vector<PrimaryVariables> initialValues_;
-};
-
-} // end namespace Dumux
-```
-
-</details>
-
-</details>
-
-
-
-## The file `properties.hh`
-
-<details open>
-<summary><b>Click to hide/show the file documentation</b> (or inspect the [source code](properties.hh))</summary>
-
-
-### Includes
-The header includes will be mentioned in the text below.
-<details><summary>Click to show the header includes</summary>
-
-```cpp
-#include <dune/alugrid/grid.hh>
-
-#include <dumux/common/properties.hh>
-#include <dumux/discretization/cctpfa.hh>
-#include <dumux/porousmediumflow/2p/model.hh>
-
-#include <dumux/material/components/trichloroethene.hh>
-#include <dumux/material/components/simpleh2o.hh>
-#include <dumux/material/fluidsystems/1pliquid.hh>
-#include <dumux/material/fluidsystems/2pimmiscible.hh>
-
-#include "spatialparams.hh"
-#include "problem.hh"
-```
-
-</details>
-
-All properties are defined in the (nested) namespace
-`Dumux::Properties`. To get and set properties, we need the definitions and implementations from the
-header `dumux/common/properties.hh` included above.
-
-```cpp
-namespace Dumux::Properties {
-```
-
-First, a so-called TypeTag is created. Properties are traits specialized for this TypeTag (a simple `struct`).
-
->>>
-:white_check_mark: The properties of other TypeTags are inherited if the alias `InheritsFrom` is present.
-These other TypeTags are listed in form of a `std::tuple` in order of importance.
->>>
-
-Here, properties from the two-phase flow model (`TTag::Twop`) and the
-cell-centered finite volume scheme with two-point-flux approximation (`TTag::CCTpfaModel`)
-are inherited. These other TypeTag definitions can be found in the included
-headers `dumux/porousmediumflow/2p/model.hh` and `dumux/discretization/cctpfa.hh`.
-
-```cpp
-namespace TTag {
-struct PointSourceExample { using InheritsFrom = std::tuple<TwoP, CCTpfaModel>; };
-}
-```
-
-Next, we specialize the properties `Problem` and `SpatialParams` for our new TypeTag and
-set the type to our problem and spatial parameter classes implemented
-in `problem.hh` and `spatialparams.hh`.
-
-```cpp
-template<class TypeTag>
-struct Problem<TypeTag, TTag::PointSourceExample>
-{ using type = PointSourceProblem<TypeTag>; };
-
-template<class TypeTag>
-struct SpatialParams<TypeTag, TTag::PointSourceExample>
-{
- // two local aliases for convenience and readability
- using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
- using Scalar = GetPropType<TypeTag, Properties::Scalar>;
-
- using type = TwoPTestSpatialParams<GridGeometry, Scalar>;
-};
-```
-
-The `Grid` property tells the
-simulator to use ALUGrid - an unstructured grid manager - here
-configured for grid and coordinate dimensions `2`,
-hexahedral element types (`Dune::cube`) and non-conforming refinement mode.
-`Dune::ALUGrid` is declared in the included header `dune/alugrid/grid.hh`
-from the Dune module `dune-alugrid`.
-
-```cpp
-template<class TypeTag>
-struct Grid<TypeTag, TTag::PointSourceExample>
-{ using type = Dune::ALUGrid<2, 2, Dune::cube, Dune::nonconforming>; };
-```
-
-The `FluidSystem` property specifies which fluids are used.
-This fluid system is composed of two immiscible liquid phases which are made up
-entirely of its respective main components `SimpleH2O` (a water component with constant properties)
-and `Trichloroethene` (a DNAPL). The components, phases, and the fluid system are implemented in
-the headers `dumux/material/components/simpleh2o.hh`,
-`dumux/material/components/trichloroethene.hh`,
-`dumux/material/fluidsystems/1pliquid.hh`,
-`dumux/material/fluidsystems/2pimmiscible.hh`
-included above.
-
-```cpp
-template<class TypeTag>
-struct FluidSystem<TypeTag, TTag::PointSourceExample>
-{
- using Scalar = GetPropType<TypeTag, Properties::Scalar>;
- using WettingPhase = FluidSystems::OnePLiquid<Scalar, Components::SimpleH2O<Scalar> >;
- using NonwettingPhase = FluidSystems::OnePLiquid<Scalar, Components::Trichloroethene<Scalar> >;
-
- using type = FluidSystems::TwoPImmiscible<Scalar, WettingPhase, NonwettingPhase>;
-};
-```
-
-The two-phase model implements different primary variable formulations.
-Here we choose the pressure of the first phase and the saturation of the second phase.
-The order of phases is specified by the fluid system.
-In this case that means that the primary variables are water pressure and DNAPL saturation.
-
-```cpp
-template<class TypeTag>
-struct Formulation<TypeTag, TTag::PointSourceExample>
-{ static constexpr auto value = TwoPFormulation::p0s1; };
-
-} // end namespace Dumux::Properties
-```
-
-
-</details>
-
-
-
-## The file `main.cc`
-
-<details open>
-<summary><b>Click to hide/show the file documentation</b> (or inspect the [source code](main.cc))</summary>
-
-
-This is the main file for the 2pinfiltration example. Here we can see the programme sequence and how the system is solved using Newton's method
-### Includes
-
-```cpp
-#include <config.h>
-```
-
-Standard header file for C++, to get time and date information.
-
-```cpp
-#include <ctime>
-```
-
-Standard header file for C++, for in- and output.
-
-```cpp
-#include <iostream>
-```
-
-Dumux is based on DUNE, the Distributed and Unified Numerics Environment, which provides several grid managers and linear solvers. So we need some includes from that.
-
-```cpp
-#include <dune/common/parallel/mpihelper.hh>
-#include <dune/common/timer.hh>
-#include <dune/grid/io/file/dgfparser/dgfexception.hh>
-#include <dune/grid/io/file/vtk.hh>
-#include <dune/istl/io.hh>
-```
-
-In Dumux, a property system is used to specify the model. For this, different properties are defined containing type definitions, values and methods. All properties are declared in the file properties.hh.
-
-```cpp
-#include <dumux/common/properties.hh>
-```
-
-The following file contains the parameter class, which manages the definition of input parameters by a default value, the inputfile or the command line.
-
-```cpp
-#include <dumux/common/parameters.hh>
-```
-
-The file dumuxmessage.hh contains the class defining the start and end message of the simulation.
-
-```cpp
-#include <dumux/common/dumuxmessage.hh>
-#include <dumux/common/defaultusagemessage.hh>
-```
-
-we include the linear solver to be used to solve the linear system
-
-```cpp
-#include <dumux/linear/amgbackend.hh>
-#include <dumux/linear/linearsolvertraits.hh>
-```
-
-we include the nonlinear Newton's method
-
-```cpp
-#include <dumux/nonlinear/newtonsolver.hh>
-```
-
-Further, we include assembler, which assembles the linear systems for finite volume schemes (box-scheme, tpfa-approximation, mpfa-approximation).
-
-```cpp
-#include <dumux/assembly/fvassembler.hh>
-```
-
-The containing class in the following file defines the different differentiation methods used to compute the derivatives of the residual.
-
-```cpp
-#include <dumux/assembly/diffmethod.hh>
-```
-
-we include the spatial discretization methods available in Dumux
-
-```cpp
-#include <dumux/discretization/method.hh>
-```
-
-We need the following class to simplify the writing of dumux simulation data to VTK format.
-
-```cpp
-#include <dumux/io/vtkoutputmodule.hh>
-```
-
-The gridmanager constructs a grid from the information in the input or grid file. There is a specification for the different supported grid managers.
-
-```cpp
-#include <dumux/io/grid/gridmanager.hh>
-```
-
-We include several files which are needed for the adaptive grid
-
-```cpp
-#include <dumux/adaptive/adapt.hh>
-#include <dumux/adaptive/markelements.hh>
-#include <dumux/adaptive/initializationindicator.hh>
-#include <dumux/porousmediumflow/2p/griddatatransfer.hh>
-#include <dumux/porousmediumflow/2p/gridadaptindicator.hh>
-```
-
-Finally, we include the properties which configure the simulation
-
-```cpp
-#include "properties.hh"
-```
-
-### Beginning of the main function
-
-```cpp
-int main(int argc, char** argv) try
-{
- using namespace Dumux;
-```
-
-we define the type tag for this problem
-
-```cpp
- using TypeTag = Properties::TTag::PointSourceExample;
-```
-
-We initialize MPI, finalize is done automatically on exit
-
-```cpp
- const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
-```
-
-We print dumux start message
-
-```cpp
- if (mpiHelper.rank() == 0)
- DumuxMessage::print(/*firstCall=*/true);
-```
-
-We parse command line arguments and input file
-
-```cpp
- Parameters::init(argc, argv);
-```
-
-### Create the grid
-A gridmanager tries to create the grid either from a grid file or the input file.
-
-```cpp
- GridManager<GetPropType<TypeTag, Properties::Grid>> gridManager;
- gridManager.init();
-```
-
-The instationary non-linear problem is run on this grid.
-
-we compute on the leaf grid view
-
-```cpp
- const auto& leafGridView = gridManager.grid().leafGridView();
-```
-
-### Set-up and solving of the problem
-
-#### Set-up
-We create and initialize the finite volume grid geometry, the problem, the linear system, including the jacobian matrix, the residual and the solution vector and the gridvariables.
-
-We need the finite volume geometry to build up the subcontrolvolumes (scv) and subcontrolvolume faces (scvf) for each element of the grid partition.
-
-```cpp
- using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
- auto gridGeometry = std::make_shared<GridGeometry>(leafGridView);
- gridGeometry->update();
-```
-
-In the problem, we define the boundary and initial conditions.
-
-```cpp
- using Problem = GetPropType<TypeTag, Properties::Problem>;
- auto problem = std::make_shared<Problem>(gridGeometry);
-```
-
-We call the `computePointSourceMap` method to compute the point sources. The `computePointSourceMap` method is inherited from the fvproblem and therefore specified in the `dumux/common/fvproblem.hh`. It calls the `addPointSources` method specified in the `problem.hh` file.
-
-```cpp
- problem->computePointSourceMap();
-```
-
-We initialize the solution vector,
-
-```cpp
- using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
- SolutionVector x(gridGeometry->numDofs());
- problem->applyInitialSolution(x);
- auto xOld = x;
-```
-
-and then use the solution vector to intialize the `gridVariables`.
-
-```cpp
- using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
- auto gridVariables = std::make_shared<GridVariables>(problem, gridGeometry);
- gridVariables->init(x);
-```
-
-We instantiate the indicator for grid adaption & the data transfer, we read some parameters for indicator from the input file.
-
-```cpp
- using Scalar = GetPropType<TypeTag, Properties::Scalar>;
- const Scalar refineTol = getParam<Scalar>("Adaptive.RefineTolerance");
- const Scalar coarsenTol = getParam<Scalar>("Adaptive.CoarsenTolerance");
-```
-
-We use an indicator for a two-phase flow problem that is saturation-dependent and defined in the file `dumux/porousmediumflow/2p/gridadaptindicator.hh.` It allows to set the minimum and maximum allowed refinement levels via the input parameters.
-
-```cpp
- TwoPGridAdaptIndicator<TypeTag> indicator(gridGeometry);
-```
-
-The data transfer performs the transfer of data on a grid from before to after adaptation and is defined in the file `dumux/porousmediumflow/2p/griddatatransfer.hh`. Its main functions are to store and reconstruct the primary variables.
-
-```cpp
- TwoPGridDataTransfer<TypeTag> dataTransfer(problem, gridGeometry, gridVariables, x);
-```
-
-We do an initial refinement around sources/BCs. We use the `GridAdaptInitializationIndicator` defined in `dumux/adaptive/initializationindicator.hh` for that.
-
-```cpp
- GridAdaptInitializationIndicator<TypeTag> initIndicator(problem, gridGeometry, gridVariables);
-```
-
-We refine up to the maximum level. For every level, the indicator used for the refinement/coarsening is calculated. If any grid cells have to be adapted, the gridvariables and the pointsourcemap are updated.
-
-```cpp
- const auto maxLevel = getParam<std::size_t>("Adaptive.MaxLevel", 0);
- for (std::size_t i = 0; i < maxLevel; ++i)
- {
-```
-
-we calculate the initial indicator for adaption for each grid cell using the initial solution x
-
-```cpp
- initIndicator.calculate(x);
-```
-
-and then we mark the elements that were adapted.
-
-```cpp
- bool wasAdapted = false;
- if (markElements(gridManager.grid(), initIndicator))
- wasAdapted = adapt(gridManager.grid(), dataTransfer);
-```
-
-In case of a grid adaptation, the gridvariables and the pointsourcemap are updated.
-
-```cpp
- if (wasAdapted)
- {
-```
-
-We overwrite the old solution with the new (resized & interpolated) one
-
-```cpp
- xOld = x;
-```
-
-We initialize the secondary variables to the new (and "new old") solution
-
-```cpp
- gridVariables->updateAfterGridAdaption(x);
-```
-
-we update the point source map after adaption
-
-```cpp
- problem->computePointSourceMap();
- }
- }
-```
-
-Depending on the initial conditions, another grid adaptation might be necessary. The gridadaptindicator uses the input parameters `Adaptive.RefineTolerance` and `Adaptive.CoarsenTolerance` for this step.
-
-```cpp
- indicator.calculate(x, refineTol, coarsenTol);
-```
-
-we mark the elements that were adapted
-
-```cpp
- bool wasAdapted = false;
- if (markElements(gridManager.grid(), indicator))
- wasAdapted = adapt(gridManager.grid(), dataTransfer);
-```
-
-In case of an additional grid adaptation, the gridvariables and the pointsourcemap are updated again.
-
-```cpp
- if (wasAdapted)
- {
-```
-
-Overwrite the old solution with the new (resized & interpolated) one
-
-```cpp
- xOld = x;
-```
-
-Initialize the secondary variables to the new (and "new old") solution
-
-```cpp
- gridVariables->updateAfterGridAdaption(x);
-```
-
-Update the point source map
-
-```cpp
- problem->computePointSourceMap();
- }
-```
-
-We get some time loop parameters from the input file params.input
-
-```cpp
- using Scalar = GetPropType<TypeTag, Properties::Scalar>;
- const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
- const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
- auto dt = getParam<Scalar>("TimeLoop.DtInitial");
-```
-
-and initialize the vtkoutput. Each model has a predefined model specific output with relevant parameters for that model.
-
-```cpp
- using IOFields = GetPropType<TypeTag, Properties::IOFields>;
- VtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables, x, problem->name());
- using VelocityOutput = GetPropType<TypeTag, Properties::VelocityOutput>;
- vtkWriter.addVelocityOutput(std::make_shared<VelocityOutput>(*gridVariables));
- IOFields::initOutputModule(vtkWriter); // Add model specific output fields
- vtkWriter.write(0.0);
-```
-
-we instantiate the time loop
-
-```cpp
- auto timeLoop = std::make_shared<TimeLoop<Scalar>>(0, dt, tEnd);
- timeLoop->setMaxTimeStepSize(maxDt);
-```
-
-we set the assembler with the time loop because we have an instationary problem
-
-```cpp
- using Assembler = FVAssembler<TypeTag, DiffMethod::numeric>;
- auto assembler = std::make_shared<Assembler>(problem, gridGeometry, gridVariables, timeLoop, xOld);
-```
-
-we set the linear solver
-
-```cpp
- using LinearSolver = AMGBiCGSTABBackend<LinearSolverTraits<GridGeometry>>;
- auto linearSolver = std::make_shared<LinearSolver>(leafGridView, gridGeometry->dofMapper());
-```
-
-additionally we set the non-linear solver.
-
-```cpp
- using NewtonSolver = Dumux::NewtonSolver<Assembler, LinearSolver>;
- NewtonSolver nonLinearSolver(assembler, linearSolver);
-```
-
-we start the time loop. In each time step before we start calculating a new solution we check if we have to refine the mesh again based on the solution of the previous time step.
-
-```cpp
- timeLoop->start(); do
- {
-```
-
-We only want to refine/coarsen after first time step is finished, not before. The initial refinement was already done before the start of the time loop. This means we only refine when the time is greater than 0.
-
-```cpp
- if (timeLoop->time() > 0)
- {
-```
-
-again we compute the refinement indicator with the `TwoPGridAdaptIndicator`
-
-```cpp
- indicator.calculate(x, refineTol, coarsenTol);
-```
-
-we mark elements and adapt grid if necessary
-
-```cpp
- wasAdapted = false;
- if (markElements(gridManager.grid(), indicator))
- wasAdapted = adapt(gridManager.grid(), dataTransfer);
-```
-
-In case of a grid adaptation, the gridvariables and the pointsourcemap are updated again.
-
-```cpp
- if (wasAdapted)
- {
-```
-
-We overwrite the old solution with the new (resized & interpolated) one
-
-```cpp
- xOld = x;
-```
-
-We tell the assembler to resize the matrix and set pattern
-
-```cpp
- assembler->setJacobianPattern();
-```
-
-We tell the assembler to resize the residual
-
-```cpp
- assembler->setResidualSize();
-```
-
-We initialize the secondary variables to the new (and "new old") solution
-
-```cpp
- gridVariables->updateAfterGridAdaption(x);
-```
-
-We update the point source map
-
-```cpp
- problem->computePointSourceMap();
- }
-```
-
-we leave the refinement step
-
-```cpp
- }
-```
-
-Now, we start to calculate the new solution of that time step. First, we define the old solution as the solution of the previous time step for storage evaluations.
-
-```cpp
- assembler->setPreviousSolution(xOld);
-```
-
-We solve the non-linear system with time step control.
-
-```cpp
- nonLinearSolver.solve(x, *timeLoop);
-```
-
-We make the new solution the old solution.
-
-```cpp
- xOld = x;
- gridVariables->advanceTimeStep();
-```
-
-We advance to the time loop to the next step.
-
-```cpp
- timeLoop->advanceTimeStep();
-```
-
-We write vtk output for each time step
-
-```cpp
- vtkWriter.write(timeLoop->time());
-```
-
-We report statistics of this time step
-
-```cpp
- timeLoop->reportTimeStep();
-```
-
-We set a new dt as suggested by the newton solver for the next time step
-
-```cpp
- timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()));
-
- } while (!timeLoop->finished());
-
- timeLoop->finalize(leafGridView.comm());
-```
-
-### Final Output
-
-print dumux end message
-
-```cpp
- if (mpiHelper.rank() == 0)
- {
- Parameters::print();
- DumuxMessage::print(/*firstCall=*/false);
- }
-
- return 0;
-} // end main
-catch (Dumux::ParameterException &e)
-{
- std::cerr << std::endl << e << " ---> Abort!" << std::endl;
- return 1;
-}
-catch (Dune::DGFException & e)
-{
- std::cerr << "DGF exception thrown (" << e <<
- "). Most likely, the DGF file name is wrong "
- "or the DGF file is corrupted, "
- "e.g. missing hash at end of file or wrong number (dimensions) of entries."
- << " ---> Abort!" << std::endl;
- return 2;
-}
-catch (Dune::Exception &e)
-{
- std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl;
- return 3;
-}
-catch (...)
-{
- std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl;
- return 4;
-}
-```
+| [:arrow_right: Click to continue with part 1 of the documentation](doc/2p.md) |
+|---:|
-</details>
+## Part 2: Main program flow
+| [:arrow_right: Click to continue with part 2 of the documentation](doc/main.md) |
+|---:|
\ No newline at end of file
diff --git a/examples/2pinfiltration/doc/2p.md b/examples/2pinfiltration/doc/2p.md
new file mode 100644
index 0000000000000000000000000000000000000000..a153ad9799ddc5bd9eb07930e19f352186d0dca8
--- /dev/null
+++ b/examples/2pinfiltration/doc/2p.md
@@ -0,0 +1,529 @@
+<!-- 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) |
+|---|---:|
+
+# Part 1: Implementation of the 2p set-up
+
+The two-phase flow set-up is implemented in the files `properties.hh`,
+`problem.hh` and `spatialparams.hh`. In the file `problem.hh` a description of the boundary and initial conditions can be found. Additionally, you can see how to implement an injection well as a point source and how to read the initial solution from a text file. In the file `spatialparams.hh` we create a lens within the porous medium, that has different spatial parameters than the surrounding material and set the parameters for the
+$`p_c - S_w`$ relationship.
+
+[[_TOC_]]
+
+
+## The file `properties.hh`
+
+<details open>
+<summary><b>Click to hide/show the file documentation</b> (or inspect the [source code](../properties.hh))</summary>
+
+
+### Includes
+The header includes will be mentioned in the text below.
+<details><summary>Click to show the header includes</summary>
+
+```cpp
+#include <dune/alugrid/grid.hh>
+
+#include <dumux/common/properties.hh>
+#include <dumux/discretization/cctpfa.hh>
+#include <dumux/porousmediumflow/2p/model.hh>
+
+#include <dumux/material/components/trichloroethene.hh>
+#include <dumux/material/components/simpleh2o.hh>
+#include <dumux/material/fluidsystems/1pliquid.hh>
+#include <dumux/material/fluidsystems/2pimmiscible.hh>
+
+#include "spatialparams.hh"
+#include "problem.hh"
+```
+
+</details>
+
+### Type tag definition
+All properties are defined in the (nested) namespace
+`Dumux::Properties`. To get and set properties, we need the definitions and implementations from the
+header `dumux/common/properties.hh` included above.
+
+```cpp
+namespace Dumux::Properties {
+```
+
+First, a so-called TypeTag is created. Properties are traits specialized for this TypeTag (a simple `struct`).
+
+>>>
+:white_check_mark: The properties of other TypeTags are inherited if the alias `InheritsFrom` is present.
+These other TypeTags are listed in form of a `std::tuple` in order of importance.
+>>>
+
+Here, properties from the two-phase flow model (`TTag::Twop`) and the
+cell-centered finite volume scheme with two-point-flux approximation (`TTag::CCTpfaModel`)
+are inherited. These other TypeTag definitions can be found in the included
+headers `dumux/porousmediumflow/2p/model.hh` and `dumux/discretization/cctpfa.hh`.
+
+```cpp
+namespace TTag {
+struct PointSourceExample { using InheritsFrom = std::tuple<TwoP, CCTpfaModel>; };
+}
+```
+
+### Property specializations
+Next, we specialize the properties `Problem` and `SpatialParams` for our new TypeTag and
+set the type to our problem and spatial parameter classes implemented
+in `problem.hh` and `spatialparams.hh`.
+
+```cpp
+template<class TypeTag>
+struct Problem<TypeTag, TTag::PointSourceExample>
+{ using type = PointSourceProblem<TypeTag>; };
+
+template<class TypeTag>
+struct SpatialParams<TypeTag, TTag::PointSourceExample>
+{
+ // two local aliases for convenience and readability
+ using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+
+ using type = TwoPTestSpatialParams<GridGeometry, Scalar>;
+};
+```
+
+The `Grid` property tells the
+simulator to use ALUGrid - an unstructured grid manager - here
+configured for grid and coordinate dimensions `2`,
+hexahedral element types (`Dune::cube`) and non-conforming refinement mode.
+`Dune::ALUGrid` is declared in the included header `dune/alugrid/grid.hh`
+from the Dune module `dune-alugrid`.
+
+```cpp
+template<class TypeTag>
+struct Grid<TypeTag, TTag::PointSourceExample>
+{ using type = Dune::ALUGrid<2, 2, Dune::cube, Dune::nonconforming>; };
+```
+
+The `FluidSystem` property specifies which fluids are used.
+This fluid system is composed of two immiscible liquid phases which are made up
+entirely of its respective main components `SimpleH2O` (a water component with constant properties)
+and `Trichloroethene` (a DNAPL). The components, phases, and the fluid system are implemented in
+the headers `dumux/material/components/simpleh2o.hh`,
+`dumux/material/components/trichloroethene.hh`,
+`dumux/material/fluidsystems/1pliquid.hh`,
+`dumux/material/fluidsystems/2pimmiscible.hh`
+included above.
+
+```cpp
+template<class TypeTag>
+struct FluidSystem<TypeTag, TTag::PointSourceExample>
+{
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using WettingPhase = FluidSystems::OnePLiquid<Scalar, Components::SimpleH2O<Scalar> >;
+ using NonwettingPhase = FluidSystems::OnePLiquid<Scalar, Components::Trichloroethene<Scalar> >;
+
+ using type = FluidSystems::TwoPImmiscible<Scalar, WettingPhase, NonwettingPhase>;
+};
+```
+
+The two-phase model implements different primary variable formulations.
+Here we choose the pressure of the first phase and the saturation of the second phase.
+The order of phases is specified by the fluid system.
+In this case that means that the primary variables are water pressure and DNAPL saturation.
+
+```cpp
+template<class TypeTag>
+struct Formulation<TypeTag, TTag::PointSourceExample>
+{ static constexpr auto value = TwoPFormulation::p0s1; };
+
+} // end namespace Dumux::Properties
+```
+
+
+</details>
+
+
+
+## The file `problem.hh`
+
+<details open>
+<summary><b>Click to hide/show the file documentation</b> (or inspect the [source code](../problem.hh))</summary>
+
+
+### Includes
+We start with includes for `PorousMediumFlowProblem` and `readFileToContainer` (used below).
+
+```cpp
+#include <dumux/porousmediumflow/problem.hh>
+#include <dumux/io/container.hh>
+```
+
+### Problem class
+The problem class `PointSourceProblem` implements boundary and initial conditions.
+It derives from the `PorousMediumFlowProblem` class.
+
+```cpp
+namespace Dumux {
+
+template <class TypeTag>
+class PointSourceProblem : public PorousMediumFlowProblem<TypeTag>
+{
+```
+
+The class implementation starts with some alias declarations and index definitions for convenience
+<details><summary>Click to show local alias declarations and indices</summary>
+
+```cpp
+ using ParentType = PorousMediumFlowProblem<TypeTag>;
+ using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using Vertex = typename GridView::template Codim<GridView::dimensionworld>::Entity;
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+ using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
+ using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
+ using PointSource = GetPropType<TypeTag, Properties::PointSource>;
+ using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+ using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
+ using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
+
+ enum {
+ pressureH2OIdx = Indices::pressureIdx,
+ saturationDNAPLIdx = Indices::saturationIdx,
+ contiDNAPLEqIdx = Indices::conti0EqIdx + FluidSystem::comp1Idx,
+ waterPhaseIdx = FluidSystem::phase0Idx,
+ dnaplPhaseIdx = FluidSystem::phase1Idx
+ };
+```
+
+</details>
+
+In the constructor of the class, we call the parent type's constructor
+and read the intial values for the primary variables from a text file.
+The function `readFileToContainer` is implemented in the header `dumux/io/container.hh`.
+
+```cpp
+public:
+ PointSourceProblem(std::shared_ptr<const GridGeometry> gridGeometry)
+ : ParentType(gridGeometry)
+ {
+ initialValues_ = readFileToContainer<std::vector<PrimaryVariables>>("initialsolutioncc.txt");
+ }
+```
+
+For isothermal problems, Dumux requires problem classes to implement a `temperature()`
+member function. Fluid properties that depend on temperature will be calculated with the specified temperature.
+
+```cpp
+ Scalar temperature() const
+ {
+ return 293.15; // 10°C
+ }
+```
+
+#### Boundary types
+We define the type of boundary conditions depending on location. Two types of boundary conditions
+can be specified: Dirichlet or Neumann boundary condition. On a Dirichlet boundary, the values of the
+primary variables need to be fixed. On a Neumann boundary condition, values for derivatives need to be fixed.
+Mixed boundary conditions (different types for different equations on the same boundary) are not accepted.
+
+```cpp
+ BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
+ {
+ BoundaryTypes values;
+ // Dirichlet boundaries on the left and right hand side of the domain
+ if (onLeftBoundary_(globalPos) || onRightBoundary_(globalPos))
+ values.setAllDirichlet();
+ // and Neumann boundaries otherwise (top and bottom of the domain)
+ else
+ values.setAllNeumann();
+ return values;
+ }
+```
+
+#### Dirichlet boundaries
+We specify the values for the Dirichlet boundaries, depending on location.
+We need to fix values for the two primary variables: the water pressure
+and the DNAPL saturation.
+
+```cpp
+ PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
+ {
+ // To determine the density of water for a given state, we build a fluid state with the given conditions:
+ PrimaryVariables values;
+ GetPropType<TypeTag, Properties::FluidState> fluidState;
+ fluidState.setTemperature(temperature());
+ fluidState.setPressure(waterPhaseIdx, /*pressure=*/1e5);
+ fluidState.setPressure(dnaplPhaseIdx, /*pressure=*/1e5);
+
+ // The density is then calculated by the fluid system:
+ Scalar densityW = FluidSystem::density(fluidState, waterPhaseIdx);
+
+ // The water phase pressure is the hydrostatic pressure, scaled with a factor:
+ Scalar height = this->gridGeometry().bBoxMax()[1] - this->gridGeometry().bBoxMin()[1];
+ Scalar depth = this->gridGeometry().bBoxMax()[1] - globalPos[1];
+ Scalar alpha = 1 + 1.5/height;
+ Scalar width = this->gridGeometry().bBoxMax()[0] - this->gridGeometry().bBoxMin()[0];
+ Scalar factor = (width*alpha + (1.0 - alpha)*globalPos[0])/width;
+
+ values[pressureH2OIdx] = 1e5 - factor*densityW*this->spatialParams().gravity(globalPos)[1]*depth;
+ // The saturation of the DNAPL Trichlorethene is zero on our Dirichlet boundary:
+ values[saturationDNAPLIdx] = 0.0;
+
+ return values;
+ }
+```
+
+#### Neumann boundaries
+In our case, we need to specify mass fluxes for our two liquid phases.
+Negative sign means influx and the unit of the boundary flux is $`kg/(m^2 s)`$.
+On the inlet area, we set a DNAPL influx of $`0.04 kg/(m^2 s)`$. On all other
+Neumann boundaries, the boundary flux is zero.
+
+```cpp
+ NumEqVector neumannAtPos(const GlobalPosition &globalPos) const
+ {
+ NumEqVector values(0.0);
+ if (onInlet_(globalPos))
+ values[contiDNAPLEqIdx] = -0.04;
+
+ return values;
+ }
+```
+
+#### Initial conditions
+The initial condition needs to be set for all primary variables.
+Here, we take the data from the file that we read in previously.
+
+```cpp
+ PrimaryVariables initial(const Element& element) const
+ {
+ // The input data is written for a uniform grid with discretization length delta.
+ // Accordingly, we need to find the index of our cells, depending on the x and y coordinates,
+ // that corresponds to the indices of the input data set.
+ const auto delta = 0.0625;
+ unsigned int cellsX = this->gridGeometry().bBoxMax()[0]/delta;
+ const auto globalPos = element.geometry().center();
+
+ unsigned int dataIdx = std::trunc(globalPos[1]/delta) * cellsX + std::trunc(globalPos[0]/delta);
+ return initialValues_[dataIdx];
+ }
+```
+
+#### Point source
+In this scenario, we set a point source (e.g. modeling a well). The point source value can be solution dependent.
+Point sources are added by pushing them into the vector `pointSources`.
+The `PointSource` constructor takes two arguments.
+The first argument is a coordinate array containing the position in space,
+the second argument is an array of source value for each equation (in units of $`kg/s`$).
+Recall that the first eqution is the water phase mass balance
+and the second equation is the DNAPL phase mass balance.
+
+```cpp
+ void addPointSources(std::vector<PointSource>& pointSources) const
+ {
+ pointSources.push_back(PointSource({0.502, 3.02}, {0, 0.1}));
+ }
+```
+
+In the private part of the class, we define some helper functions for
+the boundary conditions and local variables.
+<details><summary>Click to show private data members and functions</summary>
+
+```cpp
+private:
+ bool onLeftBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[0] < this->gridGeometry().bBoxMin()[0] + eps_; }
+
+ bool onRightBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[0] > this->gridGeometry().bBoxMax()[0] - eps_; }
+
+ bool onUpperBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[1] > this->gridGeometry().bBoxMax()[1] - eps_; }
+
+ bool onInlet_(const GlobalPosition &globalPos) const
+ {
+ Scalar width = this->gridGeometry().bBoxMax()[0] - this->gridGeometry().bBoxMin()[0];
+ Scalar lambda = (this->gridGeometry().bBoxMax()[0] - globalPos[0])/width;
+ return onUpperBoundary_(globalPos) && 0.5 < lambda && lambda < 2.0/3.0;
+ }
+
+ static constexpr Scalar eps_ = 1e-6;
+ std::vector<PrimaryVariables> initialValues_;
+};
+
+} // end namespace Dumux
+```
+
+</details>
+
+</details>
+
+
+
+## The file `spatialparams.hh`
+
+<details open>
+<summary><b>Click to hide/show the file documentation</b> (or inspect the [source code](../spatialparams.hh))</summary>
+
+
+### Includes
+We include the basic spatial parameters for finite volumes file from which we will inherit
+
+```cpp
+#include <dumux/material/spatialparams/fv.hh>
+```
+
+We include all laws which are needed to define the interaction between the solid matrix and the fluids, e.g. laws for capillary pressure saturation relationships.
+
+```cpp
+#include <dumux/material/fluidmatrixinteractions/2p/regularizedvangenuchten.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh>
+```
+
+### The spatial parameters class
+In the TwoPTestSpatialParams class we define all functions needed to describe the porous matrix, e.g. porosity and permeability. We inherit from the `FVSpatialParams` class, which is the base class for multiphase porous medium flow applications.
+
+```cpp
+namespace Dumux {
+
+template<class GridGeometry, class Scalar>
+class TwoPTestSpatialParams
+: public FVSpatialParams<GridGeometry, Scalar, TwoPTestSpatialParams<GridGeometry, Scalar>>
+{
+ //we introduce using declarations that are derived from the property system which we need in this class
+ using GridView = typename GridGeometry::GridView;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using FVElementGeometry = typename GridGeometry::LocalView;
+ using SubControlVolume = typename FVElementGeometry::SubControlVolume;
+ using ThisType = TwoPTestSpatialParams<GridGeometry, Scalar>;
+ using ParentType = FVSpatialParams<GridGeometry, Scalar, ThisType>;
+
+ static constexpr int dimWorld = GridView::dimensionworld;
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+
+ using EffectiveLaw = RegularizedVanGenuchten<Scalar>;
+
+public:
+ using MaterialLaw = EffToAbsLaw<EffectiveLaw>;
+ using MaterialLawParams = typename MaterialLaw::Params;
+ using PermeabilityType = Scalar;
+```
+
+Here, we get parameters for the position of the lens and porosity and permeability from the input file. Additionally, we set the parameters for the Van-Genuchten relationship.
+
+```cpp
+ TwoPTestSpatialParams(std::shared_ptr<const GridGeometry> gridGeometry)
+ : ParentType(gridGeometry)
+ {
+ // We get the position of the lens from the params.input file.
+ // The lens is defined by the position of the lower left and the upper right corner
+ lensLowerLeft_ = getParam<GlobalPosition>("SpatialParams.LensLowerLeft");
+ lensUpperRight_ = getParam<GlobalPosition>("SpatialParams.LensUpperRight");
+
+ // We set the parameters for the material law (here Van-Genuchten Law).
+ // First we set the residual saturations for the wetting phase and the non-wetting phase.
+ // lensMaterialParams_ define the material parameters for the lens while
+ // outerMaterialParams_ define material params for the rest of the domain.
+ lensMaterialParams_.setSwr(0.18);
+ lensMaterialParams_.setSnr(0.0);
+ outerMaterialParams_.setSwr(0.05);
+ outerMaterialParams_.setSnr(0.0);
+
+ //We set the parameters for the Van Genuchten law alpha and n
+ lensMaterialParams_.setVgAlpha(0.00045);
+ lensMaterialParams_.setVgn(7.3);
+ outerMaterialParams_.setVgAlpha(0.0037);
+ outerMaterialParams_.setVgn(4.7);
+
+ //Here, we get the permeabilities from the params.input file.
+ //In case that no parameter is set, the default parameters (9.05e-12 and 4.6e-10) are used
+ lensK_ = getParam<Scalar>("SpatialParams.lensK", 9.05e-12);
+ outerK_ = getParam<Scalar>("SpatialParams.outerK", 4.6e-10);
+ }
+```
+
+We define the (intrinsic) permeability $`[m^2]`$. In this test, we use element-wise distributed permeabilities.
+
+```cpp
+ template<class ElementSolution>
+ PermeabilityType permeability(const Element& element,
+ const SubControlVolume& scv,
+ const ElementSolution& elemSol) const
+ {
+ if (isInLens_(element.geometry().center()))
+ return lensK_;
+ return outerK_;
+ }
+```
+
+We set the porosity $`[-]`$ depending on the position
+
+```cpp
+ Scalar porosityAtPos(const GlobalPosition& globalPos) const
+ {
+ if (isInLens_(globalPos))
+ return 0.2;
+ return 0.4;
+ }
+```
+
+We set the parameter object for the Van Genuchten material law.
+
+```cpp
+ template<class ElementSolution>
+ const MaterialLawParams& materialLawParams(const Element& element,
+ const SubControlVolume& scv,
+ const ElementSolution& elemSol) const
+ {
+ if (isInLens_(element.geometry().center()))
+ return lensMaterialParams_;
+ return outerMaterialParams_;
+ }
+```
+
+Here we can define which phase is to be considered as the wetting phase. Here the wetting phase is the the first phase of the fluidsystem. In this case that is water.
+
+```cpp
+ template<class FluidSystem>
+ int wettingPhaseAtPos(const GlobalPosition& globalPos) const
+ { return FluidSystem::phase0Idx; }
+```
+
+The remainder of this class contains a convenient function to determine if
+a position is inside the lens and defines the data members.
+<details><summary> Click to show private data members and member functions</summary>
+
+```cpp
+private:
+ // we have a convenience definition of the position of the lens
+ bool isInLens_(const GlobalPosition &globalPos) const
+ {
+ for (int i = 0; i < dimWorld; ++i) {
+ if (globalPos[i] < lensLowerLeft_[i] + eps_ || globalPos[i] > lensUpperRight_[i] - eps_)
+ return false;
+ }
+ return true;
+ }
+
+ GlobalPosition lensLowerLeft_;
+ GlobalPosition lensUpperRight_;
+
+ Scalar lensK_;
+ Scalar outerK_;
+ MaterialLawParams lensMaterialParams_;
+ MaterialLawParams outerMaterialParams_;
+
+ static constexpr Scalar eps_ = 1.5e-7;
+};
+
+} // end namespace Dumux
+```
+
+</details>
+
+</details>
+
+
+| [:arrow_left: Back to the main documentation](../README.md) | [:arrow_right: Continue with part 2](main.md) |
+|---|---:|
+
diff --git a/examples/2pinfiltration/doc/2p_intro.md b/examples/2pinfiltration/doc/2p_intro.md
new file mode 100644
index 0000000000000000000000000000000000000000..235417ffedd8f177e4f5e6fe2905cea59525b023
--- /dev/null
+++ b/examples/2pinfiltration/doc/2p_intro.md
@@ -0,0 +1,7 @@
+# Part 1: Implementation of the 2p set-up
+
+The two-phase flow set-up is implemented in the files `properties.hh`,
+`problem.hh` and `spatialparams.hh`. In the file `problem.hh` a description of the boundary and initial conditions can be found. Additionally, you can see how to implement an injection well as a point source and how to read the initial solution from a text file. In the file `spatialparams.hh` we create a lens within the porous medium, that has different spatial parameters than the surrounding material and set the parameters for the
+$`p_c - S_w`$ relationship.
+
+[[_TOC_]]
diff --git a/examples/2pinfiltration/doc/intro.md b/examples/2pinfiltration/doc/_intro.md
similarity index 96%
rename from examples/2pinfiltration/doc/intro.md
rename to examples/2pinfiltration/doc/_intro.md
index d16dd16e1cf779942cc03d81d7376451cf28358b..6ed89c3a77bde3cd3a23e832102968c2cbe8e102 100644
--- a/examples/2pinfiltration/doc/intro.md
+++ b/examples/2pinfiltration/doc/_intro.md
@@ -4,10 +4,11 @@ __In this example, you will learn how to__
* solve a two-phase flow in porous media problem with two immiscible phases
* set boundary conditions and a simple injection well
-* specify a lens with different porous material parameters
-* use adaptive grid refinement around the saturation front
* specify a point source
* read the initial solution from a text file
+* specify a lens with different porous material parameters
+* use adaptive grid refinement around the saturation front
+
__Prerequisites__. You need [dune-alugrid](https://gitlab.dune-project.org/extensions/dune-alugrid) in order to compile and run this example.
@@ -69,3 +70,4 @@ For more information about the discretization please have a look at the [handboo
├── spatialparams.hh -> spatial parameter fields
└── initialsolutioncc.txt -> text file with initial solution
```
+In the documentations behind the links provided in the following, you will find a detailed description how the the above mentioned set-up is realized in this example.
diff --git a/examples/2pinfiltration/doc/main.md b/examples/2pinfiltration/doc/main.md
new file mode 100644
index 0000000000000000000000000000000000000000..04af4bc3d1300b158506e5759fb3dac993aa9cfd
--- /dev/null
+++ b/examples/2pinfiltration/doc/main.md
@@ -0,0 +1,419 @@
+<!-- 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](2p.md) |
+|---|---:|
+
+# Part 2: Main program flow
+
+This main program flow is implemented in the `main()` function
+of the program which is defined in the file `main.cc` described below. Additionally, we can see how to implement an adaptive grid with a saturation dependent indicator.
+
+The code documentation is structured as follows:
+
+[[_TOC_]]
+
+
+## The file `main.cc`
+
+<details open>
+<summary><b>Click to hide/show the file documentation</b> (or inspect the [source code](../main.cc))</summary>
+
+
+This is the main file for the 2pinfiltration example. Here we can see the programme sequence and how the system is solved using Newton's method
+### Included header files
+<details><summary> Click to show includes</summary>
+
+```cpp
+#include <config.h>
+
+// Standard header file for C++, to get time and date information.
+#include <ctime>
+
+// Standard header file for C++, for in- and output.
+#include <iostream>
+```
+
+Dumux is based on DUNE, the Distributed and Unified Numerics Environment, which provides several grid managers and linear solvers. So we need some includes from that.
+
+```cpp
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/timer.hh>
+#include <dune/grid/io/file/dgfparser/dgfexception.hh>
+#include <dune/grid/io/file/vtk.hh>
+#include <dune/istl/io.hh>
+```
+
+In Dumux, a property system is used to specify the model. For this, different properties are defined containing type definitions, values and methods. All properties are declared in the file properties.hh. Additionally, we include the parameter class, which manages the definition of input parameters by a default value, the inputfile or the command line.
+
+```cpp
+#include <dumux/common/properties.hh>
+#include <dumux/common/parameters.hh>
+```
+
+The file dumuxmessage.hh contains the class defining the start and end message of the simulation.
+
+```cpp
+#include <dumux/common/dumuxmessage.hh>
+#include <dumux/common/defaultusagemessage.hh>
+```
+
+We include the linear solver to be used to solve the linear system and the nonlinear Newton's method
+
+```cpp
+#include <dumux/linear/amgbackend.hh>
+#include <dumux/linear/linearsolvertraits.hh>
+#include <dumux/nonlinear/newtonsolver.hh>
+```
+
+Further, we include assembler, which assembles the linear systems for finite volume schemes (box-scheme, tpfa-approximation, mpfa-approximation) and a file that defines the different differentiation methods used to compute the derivatives of the residual
+
+```cpp
+#include <dumux/assembly/fvassembler.hh>
+#include <dumux/assembly/diffmethod.hh>
+```
+
+We include the spatial discretization methods available in Dumux
+
+```cpp
+#include <dumux/discretization/method.hh>
+```
+
+We need the following class to simplify the writing of dumux simulation data to VTK format.
+
+```cpp
+#include <dumux/io/vtkoutputmodule.hh>
+```
+
+The gridmanager constructs a grid from the information in the input or grid file. There is a specification for the different supported grid managers.
+
+```cpp
+#include <dumux/io/grid/gridmanager.hh>
+```
+
+We include several files which are needed for the adaptive grid
+
+```cpp
+#include <dumux/adaptive/adapt.hh>
+#include <dumux/adaptive/markelements.hh>
+#include <dumux/adaptive/initializationindicator.hh>
+#include <dumux/porousmediumflow/2p/griddatatransfer.hh>
+#include <dumux/porousmediumflow/2p/gridadaptindicator.hh>
+```
+
+Finally, we include the properties which configure the simulation
+
+```cpp
+#include "properties.hh"
+```
+
+</details>
+
+### The main function
+At the beginning of the simulation, we create a type tag with a suitable alias for our problem. // This will contain all the properties that are needed to define the problem set-up and the model we use. Additionally, we have to create an instance of `Dune::MPIHelper` and parse the run-time arguments:
+
+```cpp
+int main(int argc, char** argv) try
+{
+ using namespace Dumux;
+
+ // we define the type tag for this problem
+ using TypeTag = Properties::TTag::PointSourceExample;
+
+ //We initialize MPI, finalize is done automatically on exit
+ const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
+
+ //We print dumux start message
+ if (mpiHelper.rank() == 0)
+ DumuxMessage::print(/*firstCall=*/true);
+
+ //We parse command line arguments and input file
+ Parameters::init(argc, argv);
+```
+
+#### Create the grid
+A gridmanager tries to create the grid either from a grid file or the input file.
+
+```cpp
+ GridManager<GetPropType<TypeTag, Properties::Grid>> gridManager;
+ gridManager.init();
+
+ // The instationary non-linear problem is run on this grid.
+ //
+ // we compute on the leaf grid view
+ const auto& leafGridView = gridManager.grid().leafGridView();
+```
+
+#### Set-up of the problem
+
+We create and initialize the finite volume grid geometry, the problem, the linear system, including the jacobian matrix, the residual and the solution vector and the gridvariables.
+
+We need the finite volume geometry to build up the subcontrolvolumes (scv) and subcontrolvolume faces (scvf) for each element of the grid partition.
+
+```cpp
+ using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
+ auto gridGeometry = std::make_shared<GridGeometry>(leafGridView);
+ gridGeometry->update();
+```
+
+In the problem, we define the boundary and initial conditions and compute the point sources. The `computePointSourceMap` method is inherited from the fvproblem and therefore specified in the `dumux/common/fvproblem.hh`. It calls the `addPointSources` method specified in the `problem.hh` file.
+
+```cpp
+ using Problem = GetPropType<TypeTag, Properties::Problem>;
+ auto problem = std::make_shared<Problem>(gridGeometry);
+ // We call the `computePointSourceMap` method to compute the point sources.
+ problem->computePointSourceMap();
+```
+
+We initialize the solution vector and then use the solution vector to intialize the `gridVariables`.
+
+```cpp
+ using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
+ SolutionVector x(gridGeometry->numDofs());
+ problem->applyInitialSolution(x);
+ auto xOld = x;
+
+ using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
+ auto gridVariables = std::make_shared<GridVariables>(problem, gridGeometry);
+ gridVariables->init(x);
+```
+
+##### Grid adaption
+We instantiate the indicator for grid adaption & the data transfer, we read some parameters for indicator from the input file.
+
+```cpp
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ const Scalar refineTol = getParam<Scalar>("Adaptive.RefineTolerance");
+ const Scalar coarsenTol = getParam<Scalar>("Adaptive.CoarsenTolerance");
+ // We use an indicator for a two-phase flow problem that is saturation-dependent.
+ // It is defined in the file `dumux/porousmediumflow/2p/gridadaptindicator.hh.`
+ // and allows to set the minimum and maximum allowed refinement levels via the input parameters.
+ TwoPGridAdaptIndicator<TypeTag> indicator(gridGeometry);
+ // The data transfer performs the transfer of data on a grid from before to after adaptation
+ // and is defined in the file `dumux/porousmediumflow/2p/griddatatransfer.hh`.
+ // Its main functions are to store and reconstruct the primary variables.
+ TwoPGridDataTransfer<TypeTag> dataTransfer(problem, gridGeometry, gridVariables, x);
+```
+
+We do an initial refinement around sources/BCs. We use the `GridAdaptInitializationIndicator` defined in `dumux/adaptive/initializationindicator.hh` for that.
+
+```cpp
+ GridAdaptInitializationIndicator<TypeTag> initIndicator(problem, gridGeometry, gridVariables);
+```
+
+We refine up to the maximum level. For every level, the indicator used for the refinement/coarsening is calculated. If any grid cells have to be adapted, the gridvariables and the pointsourcemap are updated.
+
+```cpp
+ const auto maxLevel = getParam<std::size_t>("Adaptive.MaxLevel", 0);
+ for (std::size_t i = 0; i < maxLevel; ++i)
+ {
+ //we calculate the initial indicator for adaption for each grid cell using the initial solution x
+ initIndicator.calculate(x);
+
+ //and then we mark the elements that were adapted.
+ bool wasAdapted = false;
+ if (markElements(gridManager.grid(), initIndicator))
+ wasAdapted = adapt(gridManager.grid(), dataTransfer);
+
+ // In case of a grid adaptation, the gridvariables and the pointsourcemap are updated.
+ if (wasAdapted)
+ {
+ // We overwrite the old solution with the new (resized & interpolated) one
+ xOld = x;
+ //We initialize the secondary variables to the new (and "new old") solution
+ gridVariables->updateAfterGridAdaption(x);
+ // we update the point source map after adaption
+ problem->computePointSourceMap();
+ }
+ }
+```
+
+Depending on the initial conditions, another grid adaptation might be necessary.
+The gridadaptindicator uses the input parameters `Adaptive.RefineTolerance` and `Adaptive.CoarsenTolerance` for this step.
+Again, if elements were adapted, we mark them and update gridvariables and the pointsourcemap accordingly.
+
+```cpp
+ indicator.calculate(x, refineTol, coarsenTol);
+
+ //we mark the elements that were adapted
+ bool wasAdapted = false;
+ if (markElements(gridManager.grid(), indicator))
+ wasAdapted = adapt(gridManager.grid(), dataTransfer);
+
+ // In case of an additional grid adaptation, the gridvariables and the pointsourcemap are updated again.
+ if (wasAdapted)
+ {
+ // Overwrite the old solution with the new (resized & interpolated) one
+ xOld = x;
+ // Initialize the secondary variables to the new (and "new old") solution
+ gridVariables->updateAfterGridAdaption(x);
+ // Update the point source map
+ problem->computePointSourceMap();
+ }
+```
+
+#### Solving the problem
+We get some time loop parameters from the input file params.input
+
+```cpp
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
+ const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
+ auto dt = getParam<Scalar>("TimeLoop.DtInitial");
+```
+
+and initialize the vtkoutput. Each model has a predefined model specific output with relevant parameters for that model.
+
+```cpp
+ using IOFields = GetPropType<TypeTag, Properties::IOFields>;
+ VtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables, x, problem->name());
+ using VelocityOutput = GetPropType<TypeTag, Properties::VelocityOutput>;
+ vtkWriter.addVelocityOutput(std::make_shared<VelocityOutput>(*gridVariables));
+ IOFields::initOutputModule(vtkWriter); // Add model specific output fields
+ vtkWriter.write(0.0);
+```
+
+We instantiate the time loop
+
+```cpp
+ auto timeLoop = std::make_shared<TimeLoop<Scalar>>(0, dt, tEnd);
+ timeLoop->setMaxTimeStepSize(maxDt);
+```
+
+and set the assembler with the time loop because we have an instationary problem
+
+```cpp
+ using Assembler = FVAssembler<TypeTag, DiffMethod::numeric>;
+ auto assembler = std::make_shared<Assembler>(problem, gridGeometry, gridVariables, timeLoop, xOld);
+```
+
+We set the linear solver and the non-linear solver
+
+```cpp
+ using LinearSolver = AMGBiCGSTABBackend<LinearSolverTraits<GridGeometry>>;
+ auto linearSolver = std::make_shared<LinearSolver>(leafGridView, gridGeometry->dofMapper());
+
+ using NewtonSolver = Dumux::NewtonSolver<Assembler, LinearSolver>;
+ NewtonSolver nonLinearSolver(assembler, linearSolver);
+```
+
+##### The time loop
+We start the time loop. In each time step before we start calculating a new solution we check if we have to refine the mesh again based on the solution of the previous time step. If the grid is adapted we update the gridvariables and the pointsourcemap. Afterwards, the solution of that time step is calculated.
+
+```cpp
+ timeLoop->start(); do
+ {
+ // We only want to refine/coarsen after first time step is finished, not before.
+ //The initial refinement was already done before the start of the time loop.
+ //This means we only refine when the time is greater than 0.
+ if (timeLoop->time() > 0)
+ {
+ // again we compute the refinement indicator with the `TwoPGridAdaptIndicator`
+ indicator.calculate(x, refineTol, coarsenTol);
+
+ //we mark elements and adapt grid if necessary
+ wasAdapted = false;
+ if (markElements(gridManager.grid(), indicator))
+ wasAdapted = adapt(gridManager.grid(), dataTransfer);
+
+ //In case of a grid adaptation, the gridvariables and the pointsourcemap are updated again.
+ if (wasAdapted)
+ {
+ // We overwrite the old solution with the new (resized & interpolated) one
+ xOld = x;
+ // We tell the assembler to resize the matrix and set pattern
+ assembler->setJacobianPattern();
+ // We tell the assembler to resize the residual
+ assembler->setResidualSize();
+ // We initialize the secondary variables to the new (and "new old") solution
+ gridVariables->updateAfterGridAdaption(x);
+ // We update the point source map
+ problem->computePointSourceMap();
+ }
+ // we leave the refinement step
+ }
+
+ // Now, we start to calculate the new solution of that time step.
+ //First, we define the old solution as the solution of the previous time step for storage evaluations.
+ assembler->setPreviousSolution(xOld);
+
+ // We solve the non-linear system with time step control.
+ nonLinearSolver.solve(x, *timeLoop);
+
+ //We make the new solution the old solution.
+ xOld = x;
+ gridVariables->advanceTimeStep();
+
+ //We advance to the time loop to the next step.
+ timeLoop->advanceTimeStep();
+
+ //We write vtk output for each time step
+ vtkWriter.write(timeLoop->time());
+
+ //We report statistics of this time step
+ timeLoop->reportTimeStep();
+
+ //We set a new dt as suggested by the newton solver for the next time step
+ timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()));
+
+ } while (!timeLoop->finished());
+```
+
+The following piece of code prints a final status report of the time loop
+before the program is terminated and we print he dumux end message
+
+```cpp
+ timeLoop->finalize(leafGridView.comm());
+
+ // print dumux end message
+ if (mpiHelper.rank() == 0)
+ {
+ Parameters::print();
+ DumuxMessage::print(/*firstCall=*/false);
+ }
+
+ return 0;
+} // end main
+```
+
+### Exception handling
+In this part of the main file we catch and print possible exceptions that could
+occur during the simulation.
+<details><summary> Click to show exception handler</summary>
+
+```cpp
+// errors related to run-time parameters
+catch (Dumux::ParameterException &e)
+{
+ std::cerr << std::endl << e << " ---> Abort!" << std::endl;
+ return 1;
+}
+catch (Dune::DGFException & e)
+{
+ std::cerr << "DGF exception thrown (" << e <<
+ "). Most likely, the DGF file name is wrong "
+ "or the DGF file is corrupted, "
+ "e.g. missing hash at end of file or wrong number (dimensions) of entries."
+ << " ---> Abort!" << std::endl;
+ return 2;
+}
+catch (Dune::Exception &e)
+{
+ std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl;
+ return 3;
+}
+catch (...)
+{
+ std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl;
+ return 4;
+}
+```
+
+</details>
+
+</details>
+
+
+| [:arrow_left: Back to the main documentation](../README.md) | [:arrow_left: Go back to part 1](2p.md) |
+|---|---:|
+
diff --git a/examples/2pinfiltration/doc/main_intro.md b/examples/2pinfiltration/doc/main_intro.md
new file mode 100644
index 0000000000000000000000000000000000000000..cb7b9fef0fe04ee6f7155cb444a1fc79c1b4f5bf
--- /dev/null
+++ b/examples/2pinfiltration/doc/main_intro.md
@@ -0,0 +1,8 @@
+# Part 2: Main program flow
+
+This main program flow is implemented in the `main()` function
+of the program which is defined in the file `main.cc` described below. Additionally, we can see how to implement an adaptive grid with a saturation dependent indicator.
+
+The code documentation is structured as follows:
+
+[[_TOC_]]
diff --git a/examples/2pinfiltration/main.cc b/examples/2pinfiltration/main.cc
index 1c8d0d0560459c48d68058de624b677f74535018..b57af4cbacbfe9e8bfc62d673cfd98f369982301 100644
--- a/examples/2pinfiltration/main.cc
+++ b/examples/2pinfiltration/main.cc
@@ -20,15 +20,17 @@
// [[content]]
//
// This is the main file for the 2pinfiltration example. Here we can see the programme sequence and how the system is solved using Newton's method
-// ### Includes
+// ### Included header files
+// [[details]] includes
+// [[codeblock]]
#include <config.h>
-
// Standard header file for C++, to get time and date information.
#include <ctime>
// Standard header file for C++, for in- and output.
#include <iostream>
+// [[/codeblock]]
// Dumux is based on DUNE, the Distributed and Unified Numerics Environment, which provides several grid managers and linear solvers. So we need some includes from that.
#include <dune/common/parallel/mpihelper.hh>
@@ -37,24 +39,24 @@
#include <dune/grid/io/file/vtk.hh>
#include <dune/istl/io.hh>
-// In Dumux, a property system is used to specify the model. For this, different properties are defined containing type definitions, values and methods. All properties are declared in the file properties.hh.
+// In Dumux, a property system is used to specify the model. For this, different properties are defined containing type definitions, values and methods. All properties are declared in the file properties.hh. Additionally, we include the parameter class, which manages the definition of input parameters by a default value, the inputfile or the command line.
#include <dumux/common/properties.hh>
-// The following file contains the parameter class, which manages the definition of input parameters by a default value, the inputfile or the command line.
#include <dumux/common/parameters.hh>
+
// The file dumuxmessage.hh contains the class defining the start and end message of the simulation.
#include <dumux/common/dumuxmessage.hh>
#include <dumux/common/defaultusagemessage.hh>
-//we include the linear solver to be used to solve the linear system
+//We include the linear solver to be used to solve the linear system and the nonlinear Newton's method
#include <dumux/linear/amgbackend.hh>
#include <dumux/linear/linearsolvertraits.hh>
-//we include the nonlinear Newton's method
#include <dumux/nonlinear/newtonsolver.hh>
-// Further, we include assembler, which assembles the linear systems for finite volume schemes (box-scheme, tpfa-approximation, mpfa-approximation).
+
+// Further, we include assembler, which assembles the linear systems for finite volume schemes (box-scheme, tpfa-approximation, mpfa-approximation) and a file that defines the different differentiation methods used to compute the derivatives of the residual
#include <dumux/assembly/fvassembler.hh>
-// The containing class in the following file defines the different differentiation methods used to compute the derivatives of the residual.
#include <dumux/assembly/diffmethod.hh>
-// we include the spatial discretization methods available in Dumux
+
+// We include the spatial discretization methods available in Dumux
#include <dumux/discretization/method.hh>
// We need the following class to simplify the writing of dumux simulation data to VTK format.
#include <dumux/io/vtkoutputmodule.hh>
@@ -70,8 +72,12 @@
// Finally, we include the properties which configure the simulation
#include "properties.hh"
+// [[/details]]
+//
-// ### Beginning of the main function
+// ### The main function
+// At the beginning of the simulation, we create a type tag with a suitable alias for our problem. // This will contain all the properties that are needed to define the problem set-up and the model we use. Additionally, we have to create an instance of `Dune::MPIHelper` and parse the run-time arguments:
+// [[codeblock]]
int main(int argc, char** argv) try
{
using namespace Dumux;
@@ -88,9 +94,11 @@ int main(int argc, char** argv) try
//We parse command line arguments and input file
Parameters::init(argc, argv);
+ // [[/codeblock]]
- // ### Create the grid
+ // #### Create the grid
// A gridmanager tries to create the grid either from a grid file or the input file.
+ // [[codeblock]]
GridManager<GetPropType<TypeTag, Properties::Grid>> gridManager;
gridManager.init();
@@ -98,10 +106,10 @@ int main(int argc, char** argv) try
//
// we compute on the leaf grid view
const auto& leafGridView = gridManager.grid().leafGridView();
+ // [[/codeblock]]
- // ### Set-up and solving of the problem
+ // #### Set-up of the problem
//
- // #### Set-up
// We create and initialize the finite volume grid geometry, the problem, the linear system, including the jacobian matrix, the residual and the solution vector and the gridvariables.
//
// We need the finite volume geometry to build up the subcontrolvolumes (scv) and subcontrolvolume faces (scvf) for each element of the grid partition.
@@ -109,36 +117,45 @@ int main(int argc, char** argv) try
auto gridGeometry = std::make_shared<GridGeometry>(leafGridView);
gridGeometry->update();
- // In the problem, we define the boundary and initial conditions.
+ // In the problem, we define the boundary and initial conditions and compute the point sources. The `computePointSourceMap` method is inherited from the fvproblem and therefore specified in the `dumux/common/fvproblem.hh`. It calls the `addPointSources` method specified in the `problem.hh` file.
+ // [[codeblock]]
using Problem = GetPropType<TypeTag, Properties::Problem>;
auto problem = std::make_shared<Problem>(gridGeometry);
- // We call the `computePointSourceMap` method to compute the point sources. The `computePointSourceMap` method is inherited from the fvproblem and therefore specified in the `dumux/common/fvproblem.hh`. It calls the `addPointSources` method specified in the `problem.hh` file.
+ // We call the `computePointSourceMap` method to compute the point sources.
problem->computePointSourceMap();
+ // [[/codeblock]]
- // We initialize the solution vector,
+ // We initialize the solution vector and then use the solution vector to intialize the `gridVariables`.
using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
SolutionVector x(gridGeometry->numDofs());
problem->applyInitialSolution(x);
auto xOld = x;
- // and then use the solution vector to intialize the `gridVariables`.
using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
auto gridVariables = std::make_shared<GridVariables>(problem, gridGeometry);
gridVariables->init(x);
- //We instantiate the indicator for grid adaption & the data transfer, we read some parameters for indicator from the input file.
+ // ##### Grid adaption
+ // We instantiate the indicator for grid adaption & the data transfer, we read some parameters for indicator from the input file.
+ // [[codeblock]]
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
const Scalar refineTol = getParam<Scalar>("Adaptive.RefineTolerance");
const Scalar coarsenTol = getParam<Scalar>("Adaptive.CoarsenTolerance");
- // We use an indicator for a two-phase flow problem that is saturation-dependent and defined in the file `dumux/porousmediumflow/2p/gridadaptindicator.hh.` It allows to set the minimum and maximum allowed refinement levels via the input parameters.
+ // We use an indicator for a two-phase flow problem that is saturation-dependent.
+ // It is defined in the file `dumux/porousmediumflow/2p/gridadaptindicator.hh.`
+ // and allows to set the minimum and maximum allowed refinement levels via the input parameters.
TwoPGridAdaptIndicator<TypeTag> indicator(gridGeometry);
- // The data transfer performs the transfer of data on a grid from before to after adaptation and is defined in the file `dumux/porousmediumflow/2p/griddatatransfer.hh`. Its main functions are to store and reconstruct the primary variables.
+ // The data transfer performs the transfer of data on a grid from before to after adaptation
+ // and is defined in the file `dumux/porousmediumflow/2p/griddatatransfer.hh`.
+ // Its main functions are to store and reconstruct the primary variables.
TwoPGridDataTransfer<TypeTag> dataTransfer(problem, gridGeometry, gridVariables, x);
+ // [[/codeblock]]
// We do an initial refinement around sources/BCs. We use the `GridAdaptInitializationIndicator` defined in `dumux/adaptive/initializationindicator.hh` for that.
GridAdaptInitializationIndicator<TypeTag> initIndicator(problem, gridGeometry, gridVariables);
//We refine up to the maximum level. For every level, the indicator used for the refinement/coarsening is calculated. If any grid cells have to be adapted, the gridvariables and the pointsourcemap are updated.
+ // [[codeblock]]
const auto maxLevel = getParam<std::size_t>("Adaptive.MaxLevel", 0);
for (std::size_t i = 0; i < maxLevel; ++i)
{
@@ -161,8 +178,12 @@ int main(int argc, char** argv) try
problem->computePointSourceMap();
}
}
+ // [[/codeblock]]
- // Depending on the initial conditions, another grid adaptation might be necessary. The gridadaptindicator uses the input parameters `Adaptive.RefineTolerance` and `Adaptive.CoarsenTolerance` for this step.
+ // Depending on the initial conditions, another grid adaptation might be necessary.
+ //The gridadaptindicator uses the input parameters `Adaptive.RefineTolerance` and `Adaptive.CoarsenTolerance` for this step.
+ //Again, if elements were adapted, we mark them and update gridvariables and the pointsourcemap accordingly.
+ // [[codeblock]]
indicator.calculate(x, refineTol, coarsenTol);
//we mark the elements that were adapted
@@ -170,7 +191,7 @@ int main(int argc, char** argv) try
if (markElements(gridManager.grid(), indicator))
wasAdapted = adapt(gridManager.grid(), dataTransfer);
- //In case of an additional grid adaptation, the gridvariables and the pointsourcemap are updated again.
+ // In case of an additional grid adaptation, the gridvariables and the pointsourcemap are updated again.
if (wasAdapted)
{
// Overwrite the old solution with the new (resized & interpolated) one
@@ -180,14 +201,17 @@ int main(int argc, char** argv) try
// Update the point source map
problem->computePointSourceMap();
}
+ // [[/codeblock]]
+
+ // #### Solving the problem
- //We get some time loop parameters from the input file params.input
+ // We get some time loop parameters from the input file params.input
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
auto dt = getParam<Scalar>("TimeLoop.DtInitial");
- //and initialize the vtkoutput. Each model has a predefined model specific output with relevant parameters for that model.
+ // and initialize the vtkoutput. Each model has a predefined model specific output with relevant parameters for that model.
using IOFields = GetPropType<TypeTag, Properties::IOFields>;
VtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables, x, problem->name());
using VelocityOutput = GetPropType<TypeTag, Properties::VelocityOutput>;
@@ -195,26 +219,29 @@ int main(int argc, char** argv) try
IOFields::initOutputModule(vtkWriter); // Add model specific output fields
vtkWriter.write(0.0);
- //we instantiate the time loop
+ // We instantiate the time loop
auto timeLoop = std::make_shared<TimeLoop<Scalar>>(0, dt, tEnd);
timeLoop->setMaxTimeStepSize(maxDt);
- //we set the assembler with the time loop because we have an instationary problem
+ // and set the assembler with the time loop because we have an instationary problem
using Assembler = FVAssembler<TypeTag, DiffMethod::numeric>;
auto assembler = std::make_shared<Assembler>(problem, gridGeometry, gridVariables, timeLoop, xOld);
- // we set the linear solver
+ // We set the linear solver and the non-linear solver
using LinearSolver = AMGBiCGSTABBackend<LinearSolverTraits<GridGeometry>>;
auto linearSolver = std::make_shared<LinearSolver>(leafGridView, gridGeometry->dofMapper());
- // additionally we set the non-linear solver.
using NewtonSolver = Dumux::NewtonSolver<Assembler, LinearSolver>;
NewtonSolver nonLinearSolver(assembler, linearSolver);
- //we start the time loop. In each time step before we start calculating a new solution we check if we have to refine the mesh again based on the solution of the previous time step.
+ // ##### The time loop
+ // We start the time loop. In each time step before we start calculating a new solution we check if we have to refine the mesh again based on the solution of the previous time step. If the grid is adapted we update the gridvariables and the pointsourcemap. Afterwards, the solution of that time step is calculated.
+ // [[codeblock]]
timeLoop->start(); do
{
- // We only want to refine/coarsen after first time step is finished, not before. The initial refinement was already done before the start of the time loop. This means we only refine when the time is greater than 0.
+ // We only want to refine/coarsen after first time step is finished, not before.
+ //The initial refinement was already done before the start of the time loop.
+ //This means we only refine when the time is greater than 0.
if (timeLoop->time() > 0)
{
// again we compute the refinement indicator with the `TwoPGridAdaptIndicator`
@@ -242,7 +269,8 @@ int main(int argc, char** argv) try
// we leave the refinement step
}
- // Now, we start to calculate the new solution of that time step. First, we define the old solution as the solution of the previous time step for storage evaluations.
+ // Now, we start to calculate the new solution of that time step.
+ //First, we define the old solution as the solution of the previous time step for storage evaluations.
assembler->setPreviousSolution(xOld);
// We solve the non-linear system with time step control.
@@ -265,11 +293,13 @@ int main(int argc, char** argv) try
timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()));
} while (!timeLoop->finished());
+ // [[/codeblock]]
+ // The following piece of code prints a final status report of the time loop
+ // before the program is terminated and we print he dumux end message
+ // [[codeblock]]
timeLoop->finalize(leafGridView.comm());
- // ### Final Output
- //
// print dumux end message
if (mpiHelper.rank() == 0)
{
@@ -279,6 +309,13 @@ int main(int argc, char** argv) try
return 0;
} // end main
+//[[/codeblock]]
+// ### Exception handling
+// In this part of the main file we catch and print possible exceptions that could
+// occur during the simulation.
+// [[details]] exception handler
+// [[codeblock]]
+// errors related to run-time parameters
catch (Dumux::ParameterException &e)
{
std::cerr << std::endl << e << " ---> Abort!" << std::endl;
@@ -303,4 +340,6 @@ catch (...)
std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl;
return 4;
}
+// [[/codeblock]]
+// [[/details]]
// [[/content]]
diff --git a/examples/2pinfiltration/problem.hh b/examples/2pinfiltration/problem.hh
index 8fb6e92095e38d3d0be301fd896385977282dff1..feae64923bed1c85dadc8786ec5e854f73f18030 100644
--- a/examples/2pinfiltration/problem.hh
+++ b/examples/2pinfiltration/problem.hh
@@ -78,7 +78,7 @@ public:
return 293.15; // 10°C
}
- // ### 1. Boundary types
+ // #### Boundary types
// We define the type of boundary conditions depending on location. Two types of boundary conditions
// can be specified: Dirichlet or Neumann boundary condition. On a Dirichlet boundary, the values of the
// primary variables need to be fixed. On a Neumann boundary condition, values for derivatives need to be fixed.
@@ -97,7 +97,7 @@ public:
}
// [[/codeblock]]
- // ### 2. Dirichlet boundaries
+ // #### Dirichlet boundaries
// We specify the values for the Dirichlet boundaries, depending on location.
// We need to fix values for the two primary variables: the water pressure
// and the DNAPL saturation.
@@ -129,7 +129,7 @@ public:
}
// [[/codeblock]]
- // ### 3. Neumann boundaries
+ // #### Neumann boundaries
// In our case, we need to specify mass fluxes for our two liquid phases.
// Negative sign means influx and the unit of the boundary flux is $`kg/(m^2 s)`$.
// On the inlet area, we set a DNAPL influx of $`0.04 kg/(m^2 s)`$. On all other
@@ -145,7 +145,7 @@ public:
}
// [[/codeblock]]
- // ### 4. Initial conditions
+ // #### Initial conditions
// The initial condition needs to be set for all primary variables.
// Here, we take the data from the file that we read in previously.
// [[codeblock]]
@@ -163,7 +163,7 @@ public:
}
// [[/codeblock]]
- // ### 5. Point source
+ // #### Point source
// In this scenario, we set a point source (e.g. modeling a well). The point source value can be solution dependent.
// Point sources are added by pushing them into the vector `pointSources`.
// The `PointSource` constructor takes two arguments.
diff --git a/examples/2pinfiltration/properties.hh b/examples/2pinfiltration/properties.hh
index 31a977f21548d5bee363d91d5d454975f45cc838..9fe635dcd4f9d53e0d3ec8cf51a25781aa226861 100644
--- a/examples/2pinfiltration/properties.hh
+++ b/examples/2pinfiltration/properties.hh
@@ -41,7 +41,7 @@
#include "problem.hh"
// </details>
//
-
+// ### Type tag definition
// All properties are defined in the (nested) namespace
// `Dumux::Properties`. To get and set properties, we need the definitions and implementations from the
// header `dumux/common/properties.hh` included above.
@@ -61,7 +61,7 @@ namespace Dumux::Properties {
namespace TTag {
struct PointSourceExample { using InheritsFrom = std::tuple<TwoP, CCTpfaModel>; };
}
-
+// ### Property specializations
// Next, we specialize the properties `Problem` and `SpatialParams` for our new TypeTag and
// set the type to our problem and spatial parameter classes implemented
// in `problem.hh` and `spatialparams.hh`.
diff --git a/examples/2pinfiltration/spatialparams.hh b/examples/2pinfiltration/spatialparams.hh
index 14cd6a117e511db4e19081bbeafe02fa5f466349..ee1a68c240aa80b589c5f54623e8589b15d2e47d 100644
--- a/examples/2pinfiltration/spatialparams.hh
+++ b/examples/2pinfiltration/spatialparams.hh
@@ -24,16 +24,18 @@
// [[content]]
//
// ### Includes
-// we include the basic spatial parameters for finite volumes file from which we will inherit
+// We include the basic spatial parameters for finite volumes file from which we will inherit
#include <dumux/material/spatialparams/fv.hh>
-// we include all laws which are needed to define the interaction between the solid matrix and the fluids, e.g. laws for capillary pressure saturation relationships.
+// We include all laws which are needed to define the interaction between the solid matrix and the fluids, e.g. laws for capillary pressure saturation relationships.
#include <dumux/material/fluidmatrixinteractions/2p/regularizedvangenuchten.hh>
#include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh>
+// ### The spatial parameters class
+// In the TwoPTestSpatialParams class we define all functions needed to describe the porous matrix, e.g. porosity and permeability. We inherit from the `FVSpatialParams` class, which is the base class for multiphase porous medium flow applications.
+// [[codeblock]]
namespace Dumux {
-// In the TwoPTestSpatialParams class we define all functions needed to describe the porous matrix, e.g. porosity and permeability
template<class GridGeometry, class Scalar>
class TwoPTestSpatialParams
: public FVSpatialParams<GridGeometry, Scalar, TwoPTestSpatialParams<GridGeometry, Scalar>>
@@ -55,30 +57,39 @@ public:
using MaterialLaw = EffToAbsLaw<EffectiveLaw>;
using MaterialLawParams = typename MaterialLaw::Params;
using PermeabilityType = Scalar;
+ // [[/codeblock]]
+ //Here, we get parameters for the position of the lens and porosity and permeability from the input file. Additionally, we set the parameters for the Van-Genuchten relationship.
+ // [[codeblock]]
TwoPTestSpatialParams(std::shared_ptr<const GridGeometry> gridGeometry)
: ParentType(gridGeometry)
{
- //we get the position of the lens from the params.input file. The lens is defined by the position of the lower left and the upper right corner
+ // We get the position of the lens from the params.input file.
+ // The lens is defined by the position of the lower left and the upper right corner
lensLowerLeft_ = getParam<GlobalPosition>("SpatialParams.LensLowerLeft");
lensUpperRight_ = getParam<GlobalPosition>("SpatialParams.LensUpperRight");
- //we set the parameters for the material law (here Van-Genuchten Law). First we set the residual saturations for the wetting phase and the non-wetting phase. lensMaterialParams_ define the material parameters for the lens while outerMaterialParams_ define material params for the rest of the domain.
+ // We set the parameters for the material law (here Van-Genuchten Law).
+ // First we set the residual saturations for the wetting phase and the non-wetting phase.
+ // lensMaterialParams_ define the material parameters for the lens while
+ // outerMaterialParams_ define material params for the rest of the domain.
lensMaterialParams_.setSwr(0.18);
lensMaterialParams_.setSnr(0.0);
outerMaterialParams_.setSwr(0.05);
outerMaterialParams_.setSnr(0.0);
- //we set the parameters for the Van Genuchten law alpha and n
+ //We set the parameters for the Van Genuchten law alpha and n
lensMaterialParams_.setVgAlpha(0.00045);
lensMaterialParams_.setVgn(7.3);
outerMaterialParams_.setVgAlpha(0.0037);
outerMaterialParams_.setVgn(4.7);
- //here we get the permeabilities from the params.input file. In case that no parameter is set, the default parameters (9.05e-12 and 4.6e-10) are used
+ //Here, we get the permeabilities from the params.input file.
+ //In case that no parameter is set, the default parameters (9.05e-12 and 4.6e-10) are used
lensK_ = getParam<Scalar>("SpatialParams.lensK", 9.05e-12);
outerK_ = getParam<Scalar>("SpatialParams.outerK", 4.6e-10);
}
+ // [[/codeblock]]
// We define the (intrinsic) permeability $`[m^2]`$. In this test, we use element-wise distributed permeabilities.
template<class ElementSolution>
@@ -116,6 +127,10 @@ public:
int wettingPhaseAtPos(const GlobalPosition& globalPos) const
{ return FluidSystem::phase0Idx; }
+ // The remainder of this class contains a convenient function to determine if
+ // a position is inside the lens and defines the data members.
+ // [[details]] private data members and member functions
+ // [[codeblock]]
private:
// we have a convenience definition of the position of the lens
bool isInLens_(const GlobalPosition &globalPos) const
@@ -139,5 +154,7 @@ private:
};
} // end namespace Dumux
+// [[/codeblock]]
+// [[/details]]
// [[/content]]
#endif