diff --git a/examples/README.md b/examples/README.md
index 9d3e0ba9a2a58f442e56a7fbeb7c857665f62178..cff9e3ed7c7fc52bfff2683207c4fb6dc57a7245 100644
--- a/examples/README.md
+++ b/examples/README.md
@@ -29,7 +29,7 @@ You learn how to
### [:open_file_folder: Example 3: Shallow water model](shallowwaterfriction/README.md)
-The shallow water flow model is applied to simulate steady subcritical flow in a river including a bottom friction model.
+The shallow water flow model is applied to simulate steady subcritical flow in a channel including a bottom friction model.
You learn how to
* solve a shallow water flow problem including bottom friction
diff --git a/examples/shallowwaterfriction/.doc_config b/examples/shallowwaterfriction/.doc_config
index 7ac7201b34b8aae4b5b420aec30f17eed3688bfe..d5233fd849b2a624df0db5b86dcb89df13d00c22 100644
--- a/examples/shallowwaterfriction/.doc_config
+++ b/examples/shallowwaterfriction/.doc_config
@@ -1,9 +1,29 @@
{
"README.md" : [
- "doc/intro.md",
+ "doc/_intro.md"
+ ],
+
+ "doc/swe.md" : [
+ "doc/swe_intro.md",
"properties.hh",
- "spatialparams.hh",
"problem.hh",
+ "spatialparams.hh"
+ ],
+
+ "doc/main.md" : [
+ "doc/main_intro.md",
"main.cc"
- ]
+ ],
+
+ "navigation" : {
+ "mainpage" : "README.md",
+ "subpages" : [
+ "doc/swe.md",
+ "doc/main.md"
+ ],
+ "subtitles" : [
+ "Shallow water flow simulation setup",
+ "Main program flow"
+ ]
+ }
}
diff --git a/examples/shallowwaterfriction/README.md b/examples/shallowwaterfriction/README.md
index 753614f12311ba15194053ab77c3c9473d6d0733..29cb808d13e4dbfdb03cca74e3e1ccf26b8c443f 100644
--- a/examples/shallowwaterfriction/README.md
+++ b/examples/shallowwaterfriction/README.md
@@ -1,25 +1,60 @@
<!-- Important: This file has been automatically generated by generate_example_docs.py. Do not edit this file directly! -->
# Shallow water flow with bottom friction
-This example shows how the shallow water flow model can be
-applied to simulate steady subcritical flow including
-bottom friction (bed shear stress).
+In this example, the shallow water flow model is applied to simulate
+a steady subcritical flow including bottom friction (bed shear stress).
__You will learn how to__
* solve a shallow water flow problem including bottom friction
* compute and output (VTK) an analytical reference solution
-__Result__. The numerical and analytical solutions for the problem will look like this:
+__Result__. The numerical and analytical solutions for the free surface will look like this:
-
+<figure>
+ <center>
+ <img src="img/swe_result.png" alt="Shallow water result" width="60%"/>
+ <figcaption> <b> Fig.1 </b> - Setup and result for the shallow water problem with bottom friction.</figcaption>
+ </center>
+</figure>
__Table of contents__. This description is structured as follows:
[[_TOC_]]
-## Mathematical model
-The 2D shallow water equations (SWEs) are given by
+## Problem set-up
+### Model domain
+The model domain is given by a rough channel with a slope of 0.001.
+The domain is 500 meters long and 5 meters wide.
+The bottom altitude is 10 m at the inflow and hence 9.5 m at the outflow.
+Bottom friction is considered by applying
+Manning's law ($`n`$ = 0.025).
+
+### Boundary conditions
+At the lateral sides a no-flow boundary condition is applied. Also no friction is
+considered there and therefore a no slip boundary
+condition is applied. These are the default boundary condition for the shallow
+water model. At the left border a discharge boundary condition
+is applied as inflow boundary condition with $`q = -1.0 m^2 s^{-1}`$.
+At the right border a fixed water depth boundary condition
+is applied for the outflow. Normal flow is assumed, therefore the water
+depth at the right border is calculated using the equation
+of Gauckler, Manning and Strickler.
+
+### Initial conditons
+The initial water depth is set to 1 m, which is slightly higher than the normal flow
+water depth (0.87 m). Therefore, we expect a decreasing
+water level during the simulation until the normal flow condition is reached in
+the entire model domain. The inital velocity is set to zero.
+
+## Model description
+As mentioned above, this examples uses the shallow water equations (SWEs) to solve the problem.
+These are a depth averaged simplification of the Navier-Stokes equations. To calculate the
+bottom friction Manning's law is used. An alternative is Nikuradse's law, which is also implemented
+in DuMu<sup>x</sup>.
+
+### Shallow water model
+The shallow water equations are given as:
```math
\frac{\partial \mathbf{U}}{\partial t} +
@@ -35,54 +70,61 @@ where $`\mathbf{U}`$, $`\mathbf{F}`$ and $`\mathbf{G}`$ defined as
\mathbf{G} = \begin{bmatrix} hv \\ huv \\ hv^2 + \frac{1}{2} gh^2 \end{bmatrix}
```
-$`Z`$ is the bedSurface, $`h`$ the water depth, $`u`$ the velocity in
-x-direction and $`v`$ the velocity in y-direction, $`g`$ is the constant of gravity.
+$`h`$ the water depth, $`u`$ the velocity in x-direction and $`v`$ the velocity in y-direction,
+$`g`$ is the constant of gravity.
-The source terms for the bed friction $`\mathbf{S_b}`$ and bed slope
+The source terms for the bed slope $`\mathbf{S_b}`$ and friction
$`\mathbf{S_f}`$ are given as
```math
\mathbf{S_b} = \begin{bmatrix} 0 \\ -gh \frac{\partial z}{\partial x}
\\ -gh \frac{\partial z}{\partial y}\end{bmatrix},
-\mathbf{S_f} = \begin{bmatrix} 0 \\ -ghS_{fx} \\ -ghS_{fy}\end{bmatrix}.
+\mathbf{S_f} = \begin{bmatrix} 0 \\ghS_{fx} \\ghS_{fy}\end{bmatrix}.
```
-For this example, a cell-centered finite volume method (`cctpfa`) is applied to solve the SWEs
-in combination with a fully-implicit time discretization. For cases where no sharp fronts or
-traveling waves occur it is possible to apply time steps larger than CFL number = 1 to reduce
-the computation time. Even if a steady state solution is considered, an implicit time stepping method
-is applied.
+with the bedSurface $`z`$. $`S_{fx}`$ and $`S_{fy}`$ are the bed shear stess
+components in x- and y-direction, which are calculated by Manning's law.
-## Problem set-up
+### Mannings law
+The empirical Manning model specifies the bed shear stress by the following equations:
-The model domain is given by a rough channel with a slope of 0.001.
-The domain is 500 meters long and 10 meters wide.
-.
+```math
+S_{fx} = \frac{n^2u}{R_{hy}^{4/3}} \sqrt(u^2 + v^2),
-Bottom friction is considered by applying
-the friction law of Manning (Manning n = 0.025). At the lateral sides no friction is considered and a
-no-flow no slip boundary condition is applied. This is the default boundary condition for the shallow water model.
+S_{fy} = \frac{n^2v}{R_{hy}^{4/3}} \sqrt(u^2 + v^2)
+```
+
+$`n`$ is Manning's friction value and $`R_{hy}`$ is the hydraulic radius,
+which is assumed to be equal to the water depth $`h`$.
-At the left border a discharge boundary condition
-is applied as inflow boundary condition with q = -1.0 ($`m^2 s^{-1}`$). At the right border a water fixed depth boundary condition
-is applied for the outflow. Normal flow is assumed, therefore the water depth at the right border is calculated after
-the of Gaukler-Manning-Strickler equation:
+### Analytical solution
+Since normal flow conditions are assumed, the analytic solution is calculated using the equation
+of Gauckler, Manning and Strickler:
```math
v_m = n^{-1} R_{hy}^{2/3} I_s^{1/2}
```
-Where the mean velocity $`v_m`$ is given as $`v_m = \frac{q}{h}`$,
-$`n`$ is the friction value after Manning. $`R_{hy}`$ the hydraulic radius, which is assumed to be equal to
-the water depth. $`I_s`$ is the bed slope and $`q`$ the unity inflow discharge
+Where the mean velocity $`v_m`$ is given as
+
+```math
+v_m = \frac{q}{h}
+```
+
+$`I_s`$ is the bed slope and $`q`$ the unity inflow discharge.
+
+Hence, the water depth $`h`$ can be calculated by
-The water depth h can be calculated as
```math
h = \left(\frac{n q}{\sqrt{I_s}} \right)^{3/5}
```
-The formula of Gaukler Manning and Strickler is also used to calculate the analytic solution. All parameters
-for the simulation are given in the file `params.input`.
+### Discretisation
+For this example, a cell-centered finite volume method (cctpfa) is applied to solve the SWEs
+in combination with a fully-implicit time discretization. For cases where no sharp fronts or
+traveling waves occur it is possible to apply time steps larger than CFL number = 1 to reduce
+the computation time. Even if a steady state solution is considered, an implicit time stepping method
+is applied.
# Implementation
@@ -98,925 +140,13 @@ for the simulation are given in the file `params.input`.
└── spatialparams.hh -> spatial parameter fields
```
+## Part 1: Shallow water flow simulation setup
-## The file `properties.hh`
-
-The header includes will be mentioned in the text below.
-<details><summary>Click to show the header includes</summary>
-
-```cpp
-#include <dune/grid/yaspgrid.hh>
-
-#include <dumux/common/properties.hh>
-#include <dumux/discretization/cctpfa.hh>
-#include <dumux/freeflow/shallowwater/model.hh>
-
-#include "spatialparams.hh"
-#include "problem.hh"
-```
-
-</details>
-
-Let's define the properties for our simulation
-
-```cpp
-namespace Dumux::Properties {
-```
-
-First, a so-called TypeTag is created. Properties are traits specialized for this TypeTag (a simple `struct`).
-The properties of two other TypeTags are inherited by adding the alias `InheritsFrom`.
-Here, properties from the shallow water model (`TTag::ShallowWater`) 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/freeflow/shallowwater/model.hh` and `dumux/discretization/cctpfa.hh`.
-
-```cpp
-namespace TTag {
-struct RoughChannel { using InheritsFrom = std::tuple<ShallowWater, CCTpfaModel>; };
-}
-```
-
-We use a structured Cartesian grid with tensor product structure.
-`Dune::YaspGrid` (Yet Another Structure Parallel Grid) is defined in `dune/grid/yaspgrid.hh`
-in the Dune module `dune-grid`.
-
-```cpp
-template<class TypeTag>
-struct Grid<TypeTag, TTag::RoughChannel>
-{ using type = Dune::YaspGrid<2, Dune::TensorProductCoordinates<GetPropType<TypeTag, Properties::Scalar>, 2> >; };
-```
-
-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::RoughChannel>
-{ using type = Dumux::RoughChannelProblem<TypeTag>; };
-
-template<class TypeTag>
-struct SpatialParams<TypeTag, TTag::RoughChannel>
-{
- using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
- using Scalar = GetPropType<TypeTag, Properties::Scalar>;
- using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
- using VolumeVariables = typename ElementVolumeVariables::VolumeVariables;
-
- using type = RoughChannelSpatialParams<GridGeometry, Scalar, VolumeVariables>;
-};
-```
-
-Finally, we enable caching for the grid geometry. The cache
-stores values that were already calculated for later usage.
-This makes the simulation run faster but it uses more memory.
-
-```cpp
-template<class TypeTag>
-struct EnableGridGeometryCache<TypeTag, TTag::RoughChannel>
-{ static constexpr bool value = true; };
-
-} // end namespace Dumux::Properties
-```
-
-
-
-
-## The file `spatialparams.hh`
-
-
-We include the basic spatial parameters for finite volumes file from which we will inherit
-
-```cpp
-#include <dumux/material/spatialparams/fv.hh>
-```
-
-The parameters header is needed to retrieve run-time parameters.
-
-```cpp
-#include <dumux/common/parameters.hh>
-```
-
-We include all friction laws, between we can choose for the calculation of the bottom friction source.
-
-```cpp
-#include <dumux/material/fluidmatrixinteractions/frictionlaws/frictionlaw.hh>
-#include <dumux/material/fluidmatrixinteractions/frictionlaws/manning.hh>
-#include <dumux/material/fluidmatrixinteractions/frictionlaws/nikuradse.hh>
-```
-
-We enter the namespace Dumux. All Dumux functions and classes are in a namespace Dumux, to make sure they don`t clash with symbols from other libraries you may want to use in conjunction with Dumux.
-
-```cpp
-namespace Dumux {
-```
-
-In the RoughChannelSpatialParams class we define all functions needed to describe the spatial distributed parameters.
-
-```cpp
-template<class GridGeometry, class Scalar, class VolumeVariables>
-class RoughChannelSpatialParams
-: public FVSpatialParams<GridGeometry, Scalar,
- RoughChannelSpatialParams<GridGeometry, Scalar, VolumeVariables>>
-{
-```
-
-We introduce using declarations that are derived from the property system which we need in this class
-
-```cpp
- using ThisType = RoughChannelSpatialParams<GridGeometry, Scalar, VolumeVariables>;
- using ParentType = FVSpatialParams<GridGeometry, Scalar, ThisType>;
- using GridView = typename GridGeometry::GridView;
- using FVElementGeometry = typename GridGeometry::LocalView;
- using SubControlVolume = typename FVElementGeometry::SubControlVolume;
- using Element = typename GridView::template Codim<0>::Entity;
- using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
-
-public:
-```
-
-In the constructor be read some values from the `params.input` and initialize the friciton law.
-
-```cpp
- RoughChannelSpatialParams(std::shared_ptr<const GridGeometry> gridGeometry)
- : ParentType(gridGeometry)
- {
- gravity_ = getParam<Scalar>("Problem.Gravity");
- bedSlope_ = getParam<Scalar>("Problem.BedSlope");
- frictionLawType_ = getParam<std::string>("Problem.FrictionLaw");
- initFrictionLaw();
- }
-```
-
-We initialize the friction law based on the law specified in `params.input`.
-
-```cpp
- void initFrictionLaw()
- {
- if (frictionLawType_ == "Manning")
- {
- Scalar manningN = getParam<Scalar>("Problem.ManningN");
- frictionLaw_ = std::make_unique<FrictionLawManning<VolumeVariables>>(gravity_, manningN);
- }
- else if (frictionLawType_ == "Nikuradse")
- {
- Scalar ks = getParam<Scalar>("Problem.Ks");
- frictionLaw_ = std::make_unique<FrictionLawNikuradse<VolumeVariables>>(ks);
- }
- else
- {
- std::cout<<"The FrictionLaw in params.input is unknown. Valid entries are `Manning` and `Nikuradse`!"<<std::endl;
- }
- }
-```
-
-Use this function, if you want to vary the value for the gravity.
-
-```cpp
- Scalar gravity(const GlobalPosition& globalPos) const
- {
- return gravity_;
- }
-```
-
-Use this function for a constant gravity.
-
-```cpp
- Scalar gravity() const
- {
- return gravity_;
- }
-```
-
-This function returns an object of the friction law class, which is initialized with the appropriate friction values. If you want to use different friciton values or laws, you have to use a vector of unique_ptr for `frictionLaw_` and pick the right friction law instances via the `element` argument.
-
-```cpp
- const FrictionLaw<VolumeVariables>& frictionLaw(const Element& element,
- const SubControlVolume& scv) const
- {
- return *frictionLaw_;
- }
-```
-
-Define the bed surface based on the `bedSlope_`.
-
-```cpp
- Scalar bedSurface(const Element& element,
- const SubControlVolume& scv) const
- {
- return 10.0 - element.geometry().center()[0] * bedSlope_;
- }
-
-private:
- Scalar gravity_;
- Scalar bedSlope_;
- std::string frictionLawType_;
- std::unique_ptr<FrictionLaw<VolumeVariables>> frictionLaw_;
-};
-```
-
-end of namespace Dumux.
-
-```cpp
-}
-```
-
-
-
-
-## The file `problem.hh`
-We start with includes
-
-```cpp
-#include <dumux/common/parameters.hh>
-#include <dumux/common/properties.hh>
-#include <dumux/freeflow/shallowwater/problem.hh>
-#include <dumux/freeflow/shallowwater/boundaryfluxes.hh>
-```
-
-We enter the problem class where all necessary boundary conditions and initial conditions are set for our simulation.
-As this is a shallow water problem, we inherit from the basic ShallowWaterProblem.
-
-```cpp
-namespace Dumux {
-
-template <class TypeTag>
-class RoughChannelProblem : public ShallowWaterProblem<TypeTag>
-{
-```
-
-We use convenient declarations that we derive from the property system.
-
-```cpp
- using ParentType = ShallowWaterProblem<TypeTag>;
- using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
- using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
- using Scalar = GetPropType<TypeTag, Properties::Scalar>;
- using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
- using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
- using NeumannFluxes = GetPropType<TypeTag, Properties::NumEqVector>;
- using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
- using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
- using ElementFluxVariablesCache = typename GridVariables::GridFluxVariablesCache::LocalView;
- using VolumeVariables = typename ElementVolumeVariables::VolumeVariables;
- using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
- using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
- using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
- using Element = typename GridView::template Codim<0>::Entity;
- using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
- using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
- using SubControlVolume = typename FVElementGeometry::SubControlVolume;
-
-public:
-```
-
-This is the constructor of our problem class.
-
-```cpp
- RoughChannelProblem(std::shared_ptr<const GridGeometry> gridGeometry)
- : ParentType(gridGeometry)
- {
-```
-
-We read the parameters from the params.input file.
-
-```cpp
- name_ = getParam<std::string>("Problem.Name");
- constManningN_ = getParam<Scalar>("Problem.ManningN");
- bedSlope_ = getParam<Scalar>("Problem.BedSlope");
- discharge_ = getParam<Scalar>("Problem.Discharge");
-```
-
-We calculate the outflow boundary condition using the Gauckler-Manning-Strickler formula.
-
-```cpp
- hBoundary_ = this->gauklerManningStrickler(discharge_,constManningN_,bedSlope_);
-```
-
-We initialize the analytic solution to a verctor of the appropriate size filled with zeros.
-
-```cpp
- exactWaterDepth_.resize(gridGeometry->numDofs(), 0.0);
- exactVelocityX_.resize(gridGeometry->numDofs(), 0.0);
- }
-```
-
-Get the analytical water depth
-
-```cpp
- const std::vector<Scalar>& getExactWaterDepth()
- {
- return exactWaterDepth_;
- }
-```
-
-Get the analytical velocity
-
-```cpp
- const std::vector<Scalar>& getExactVelocityX()
- {
- return exactVelocityX_;
- }
-```
-
-Get the water depth with Gauckler-Manning-Strickler
-
-```cpp
- Scalar gauklerManningStrickler(Scalar discharge, Scalar manningN, Scalar bedSlope)
- {
- using std::pow;
- using std::abs;
- using std::sqrt;
-
- return pow(abs(discharge)*manningN/sqrt(bedSlope), 0.6);
- }
-```
-
-Get the analytical solution
-
-```cpp
- void analyticalSolution()
- {
- using std::abs;
-
- for (const auto& element : elements(this->gridGeometry().gridView()))
- {
- const Scalar h = this->gauklerManningStrickler(discharge_,constManningN_,bedSlope_);
- const Scalar u = abs(discharge_)/h;
-
- const auto eIdx = this->gridGeometry().elementMapper().index(element);
- exactWaterDepth_[eIdx] = h;
- exactVelocityX_[eIdx] = u;
- }
- }
-```
-
-Get the problem name. It is used as a prefix for files generated by the simulation.
-
-```cpp
- const std::string& name() const
- {
- return name_;
- }
-```
-
-Get the source term.
-
-```cpp
- NumEqVector source(const Element& element,
- const FVElementGeometry& fvGeometry,
- const ElementVolumeVariables& elemVolVars,
- const SubControlVolume &scv) const
- {
-
- NumEqVector source (0.0);
-```
-
-In this model the bottom friction is the only source.
-
-```cpp
- source += bottomFrictionSource(element, fvGeometry, elemVolVars, scv);
-
- return source;
- }
-```
-
-Get the source term due to bottom friction.
-
-```cpp
- NumEqVector bottomFrictionSource(const Element& element,
- const FVElementGeometry& fvGeometry,
- const ElementVolumeVariables& elemVolVars,
- const SubControlVolume &scv) const
- {
- NumEqVector bottomFrictionSource(0.0);
- const auto& volVars = elemVolVars[scv];
-```
-
-For the calculation of the source term due to bottom friction the two-dimensional bottom shear stess vector is needed. This is the force per area, which works between the flow and the bed. It is calculated within the `FrictionLaw`, which is a spatialParameter. In this model the `FrictionLawManning` is used (see `params.input`).
-
-```cpp
- Dune::FieldVector<Scalar, 2> bottomShearStress = this->spatialParams().frictionLaw(element, scv).shearStress(volVars);
-```
-
-The bottom shear stress causes a pure loss of momentum. Thus the first entry of the `bottomFrictionSource`, which is related to the mass balance equation is zero. The second entry of the `bottomFricitonSource` corresponds to the momentum equation in x-direction and is therefore equal to the first, the x-component, of the `bottomShearStress`. Accordingly the third entry of the `bottomFrictionSource` is equal to the second component of the `bottomShearStress`.
-
-```cpp
- bottomFrictionSource[0] = 0.0;
- bottomFrictionSource[1] = bottomShearStress[0];
- bottomFrictionSource[2] = bottomShearStress[1];
-
- return bottomFrictionSource;
- }
-```
-
-We specify the boundary condition type.
-
-```cpp
- BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
- {
- BoundaryTypes bcTypes;
-```
-
-Since we use a weak imposition all boundary conditions are of Neumann type.
-
-```cpp
- bcTypes.setAllNeumann();
- return bcTypes;
- }
-```
-
-We specify the neumann boundary. Due to the weak imposition we calculate the flux at the
-boundary, with a Rieman solver. For this the state of a virtual cell outside of the boundary
-is needed (`boundaryStateVariables`), wich is calculated with the Riemann invariants
-(see Yoon and Kang, Finite Volume Model for Two-Dimensional Shallow Water Flows on Unstructured Grids).
-The calculation of the Riemann invariants differ depending on the type of the boundary (h, q or no-flow boundary).
-
-```cpp
- NeumannFluxes neumann(const Element& element,
- const FVElementGeometry& fvGeometry,
- const ElementVolumeVariables& elemVolVars,
- const ElementFluxVariablesCache& elemFluxVarsCache,
- const SubControlVolumeFace& scvf) const
- {
- NeumannFluxes values(0.0);
-
- const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
- const auto& insideVolVars = elemVolVars[insideScv];
- const auto& nxy = scvf.unitOuterNormal();
- const auto gravity = this->spatialParams().gravity(scvf.center());
- std::array<Scalar, 3> boundaryStateVariables;
-```
-
-Calculate the rieman invariants for imposed discharge at the left side.
-
-```cpp
- if (scvf.center()[0] < 0.0 + eps_)
- {
- boundaryStateVariables = ShallowWater::fixedDischargeBoundary(discharge_,
- insideVolVars.waterDepth(),
- insideVolVars.velocity(0),
- insideVolVars.velocity(1),
- gravity,
- nxy);
- }
-```
-
-Calculate the rieman invariants for impose water depth at the right side.
-
-```cpp
- else if (scvf.center()[0] > 100.0 - eps_)
- {
- boundaryStateVariables = ShallowWater::fixedWaterDepthBoundary(hBoundary_,
- insideVolVars.waterDepth(),
- insideVolVars.velocity(0),
- insideVolVars.velocity(1),
- gravity,
- nxy);
- }
-```
-
-Calculate the rieman invarianty for the no-flow boundary.
-
-```cpp
- else
- {
- boundaryStateVariables[0] = insideVolVars.waterDepth();
- boundaryStateVariables[1] = -insideVolVars.velocity(0);
- boundaryStateVariables[2] = -insideVolVars.velocity(1);
- }
-```
-
-We calculate the boundary fluxes based on a Riemann problem.
-
-```cpp
- auto riemannFlux = ShallowWater::riemannProblem(insideVolVars.waterDepth(),
- boundaryStateVariables[0],
- insideVolVars.velocity(0),
- boundaryStateVariables[1],
- insideVolVars.velocity(1),
- boundaryStateVariables[2],
- insideVolVars.bedSurface(),
- insideVolVars.bedSurface(),
- gravity,
- nxy);
-
- values[Indices::massBalanceIdx] = riemannFlux[0];
- values[Indices::velocityXIdx] = riemannFlux[1];
- values[Indices::velocityYIdx] = riemannFlux[2];
-
- return values;
- }
-```
-
-We set the initial conditions. In this example constant initial conditions are used. Therefore the argument `globalPos` is not needed. If you want to impose spatial variable initial conditions, you have to use the `globalPos`.
-
-```cpp
- PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
- {
- PrimaryVariables values(0.0);
-```
-
-We set the initial water depth to one meter.
+| [:arrow_right: Click to continue with part 1 of the documentation](doc/swe.md) |
+|---:|
-```cpp
- values[0] = 1.0;
-```
-
-We set the x-component of the initial velocity to zero.
-
-```cpp
- values[1] = 0.0;
-```
-
-We set the y-component of the initial velocity to zero.
-
-```cpp
- values[2] = 0.0;
-
- return values;
- };
-```
-
-\}
-
-```cpp
-
-private:
-```
-
-We declare the private variables of the problem. They are initialized in the problems constructor.
-We declare the variable for the analytic solution.
-
-```cpp
- std::vector<Scalar> exactWaterDepth_;
- std::vector<Scalar> exactVelocityX_;
-```
-
-constant friction value. An analytic solution is only available for const friction. If you want to run the simulation with a non constant friciton value (specified in the spatialParams) you have to remove the analytic solution.
-
-```cpp
- Scalar constManningN_;
-```
-
-The constant bed slope.
-
-```cpp
- Scalar bedSlope_;
-```
-
-The water depth at the outflow boundary.
-
-```cpp
- Scalar hBoundary_;
-```
-
-The discharge at the inflow boundary.
-
-```cpp
- Scalar discharge_;
-```
-
-eps is used as a small value for the definition of the boundry conditions
-
-```cpp
- static constexpr Scalar eps_ = 1.0e-6;
- std::string name_;
-};
-
-} // end namespace Dumux
-```
-
-
-
-
-## The file `main.cc`
-
-
-This is the main file for the shallow water 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>
-```
-
-We need the following class to simplify the writing of dumux simulation data to VTK format.
-
-```cpp
-#include <dumux/io/vtkoutputmodule.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>
-```
-
-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 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 newtons 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>
-```
-
-We include the properties
-
-```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::RoughChannel;
-```
-
-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();
-```
-
-We compute on the leaf grid view
-
-```cpp
- const auto& leafGridView = gridManager.grid().leafGridView();
-```
-
-### Setup and solving of the problem
-#### Setup
-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 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 solutionvector to intialize the gridVariables.
-
-```cpp
- using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
- auto gridVariables = std::make_shared<GridVariables>(problem, gridGeometry);
- gridVariables->init(x);
-```
-
-We get some time loop parameters from the input file.
-
-```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");
-```
-
-We intialize the vtk output module. 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());
-```
-
-We add the analytical solution ("exactWaterDepth" and "exactVelocityX") to the predefined specific output.
-
-```cpp
- vtkWriter.addField(problem->getExactWaterDepth(), "exactWaterDepth");
- vtkWriter.addField(problem->getExactVelocityX(), "exactVelocityX");
-```
-
-We calculate the analytic solution.
-
-```cpp
- problem->analyticalSolution();
- IOFields::initOutputModule(vtkWriter);
- vtkWriter.write(0.0);
-```
-
-We instantiate time loop.
-
-```cpp
- auto timeLoop = std::make_shared<CheckPointTimeLoop<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());
-```
-
-Additionaly, we set the non-linear solver.
-
-```cpp
- using NewtonSolver = Dumux::NewtonSolver<Assembler, LinearSolver>;
- NewtonSolver nonLinearSolver(assembler, linearSolver);
-```
-
-We set some check point at the end of the time loop. The check point is used to trigger the vtk output.
-
-```cpp
- timeLoop->setCheckPoint(tEnd);
-```
-
-We start the time loop.
-
-```cpp
- timeLoop->start(); do
- {
-```
-
-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, if we reached the check point (end of time loop)
-
-```cpp
- if (timeLoop->isCheckPoint())
- vtkWriter.write(timeLoop->time());
-```
-
-We report statistics of this time step.
-
-```cpp
- timeLoop->reportTimeStep();
-```
-
-We set new dt as suggested by newton controller for the next time step.
-
-```cpp
- timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()));
-
-
- } while (!timeLoop->finished());
-
- timeLoop->finalize(leafGridView.comm());
-```
-
-### Final Output
-We print dumux end message.
-
-```cpp
- if (mpiHelper.rank() == 0)
- {
- Parameters::print();
- DumuxMessage::print(/*firstCall=*/false);
- }
-
- return 0;
-} // end main
-
-catch (const Dumux::ParameterException &e)
-{
- std::cerr << std::endl << e << " ---> Abort!" << std::endl;
- return 1;
-}
-catch (const 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 (const 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;
-}
-```
+## 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/shallowwaterfriction/doc/_intro.md b/examples/shallowwaterfriction/doc/_intro.md
new file mode 100644
index 0000000000000000000000000000000000000000..172ca8fb0631a7515098f68e5b07d34a59c121c2
--- /dev/null
+++ b/examples/shallowwaterfriction/doc/_intro.md
@@ -0,0 +1,139 @@
+# Shallow water flow with bottom friction
+In this example, the shallow water flow model is applied to simulate
+a steady subcritical flow including bottom friction (bed shear stress).
+
+__You will learn how to__
+
+* solve a shallow water flow problem including bottom friction
+* compute and output (VTK) an analytical reference solution
+
+__Result__. The numerical and analytical solutions for the free surface will look like this:
+
+<figure>
+ <center>
+ <img src="img/swe_result.png" alt="Shallow water result" width="60%"/>
+ <figcaption> <b> Fig.1 </b> - Setup and result for the shallow water problem with bottom friction.</figcaption>
+ </center>
+</figure>
+
+__Table of contents__. This description is structured as follows:
+
+[[_TOC_]]
+
+## Problem set-up
+### Model domain
+The model domain is given by a rough channel with a slope of 0.001.
+The domain is 500 meters long and 5 meters wide.
+The bottom altitude is 10 m at the inflow and hence 9.5 m at the outflow.
+Bottom friction is considered by applying
+Manning's law ($`n`$ = 0.025).
+
+### Boundary conditions
+At the lateral sides a no-flow boundary condition is applied. Also no friction is
+considered there and therefore a no slip boundary
+condition is applied. These are the default boundary condition for the shallow
+water model. At the left border a discharge boundary condition
+is applied as inflow boundary condition with $`q = -1.0 m^2 s^{-1}`$.
+At the right border a fixed water depth boundary condition
+is applied for the outflow. Normal flow is assumed, therefore the water
+depth at the right border is calculated using the equation
+of Gauckler, Manning and Strickler.
+
+### Initial conditons
+The initial water depth is set to 1 m, which is slightly higher than the normal flow
+water depth (0.87 m). Therefore, we expect a decreasing
+water level during the simulation until the normal flow condition is reached in
+the entire model domain. The inital velocity is set to zero.
+
+## Model description
+As mentioned above, this examples uses the shallow water equations (SWEs) to solve the problem.
+These are a depth averaged simplification of the Navier-Stokes equations. To calculate the
+bottom friction Manning's law is used. An alternative is Nikuradse's law, which is also implemented
+in DuMu<sup>x</sup>.
+
+### Shallow water model
+The shallow water equations are given as:
+
+```math
+\frac{\partial \mathbf{U}}{\partial t} +
+\frac{\partial \mathbf{F}}{\partial x} +
+\frac{\partial \mathbf{G}}{\partial y} - \mathbf{S_b} - \mathbf{S_f} = 0
+```
+
+where $`\mathbf{U}`$, $`\mathbf{F}`$ and $`\mathbf{G}`$ defined as
+
+```math
+\mathbf{U} = \begin{bmatrix} h \\ uh \\ vh \end{bmatrix},
+\mathbf{F} = \begin{bmatrix} hu \\ hu^2 + \frac{1}{2} gh^2 \\ huv \end{bmatrix},
+\mathbf{G} = \begin{bmatrix} hv \\ huv \\ hv^2 + \frac{1}{2} gh^2 \end{bmatrix}
+```
+
+$`h`$ the water depth, $`u`$ the velocity in x-direction and $`v`$ the velocity in y-direction,
+$`g`$ is the constant of gravity.
+
+The source terms for the bed slope $`\mathbf{S_b}`$ and friction
+$`\mathbf{S_f}`$ are given as
+
+```math
+\mathbf{S_b} = \begin{bmatrix} 0 \\ -gh \frac{\partial z}{\partial x}
+ \\ -gh \frac{\partial z}{\partial y}\end{bmatrix},
+\mathbf{S_f} = \begin{bmatrix} 0 \\ghS_{fx} \\ghS_{fy}\end{bmatrix}.
+```
+
+with the bedSurface $`z`$. $`S_{fx}`$ and $`S_{fy}`$ are the bed shear stess
+components in x- and y-direction, which are calculated by Manning's law.
+
+### Mannings law
+The empirical Manning model specifies the bed shear stress by the following equations:
+
+```math
+S_{fx} = \frac{n^2u}{R_{hy}^{4/3}} \sqrt(u^2 + v^2),
+
+S_{fy} = \frac{n^2v}{R_{hy}^{4/3}} \sqrt(u^2 + v^2)
+```
+
+$`n`$ is Manning's friction value and $`R_{hy}`$ is the hydraulic radius,
+which is assumed to be equal to the water depth $`h`$.
+
+### Analytical solution
+Since normal flow conditions are assumed, the analytic solution is calculated using the equation
+of Gauckler, Manning and Strickler:
+
+```math
+v_m = n^{-1} R_{hy}^{2/3} I_s^{1/2}
+```
+
+Where the mean velocity $`v_m`$ is given as
+
+```math
+v_m = \frac{q}{h}
+```
+
+$`I_s`$ is the bed slope and $`q`$ the unity inflow discharge.
+
+Hence, the water depth $`h`$ can be calculated by
+
+```math
+h = \left(\frac{n q}{\sqrt{I_s}} \right)^{3/5}
+```
+
+### Discretisation
+For this example, a cell-centered finite volume method (cctpfa) is applied to solve the SWEs
+in combination with a fully-implicit time discretization. For cases where no sharp fronts or
+traveling waves occur it is possible to apply time steps larger than CFL number = 1 to reduce
+the computation time. Even if a steady state solution is considered, an implicit time stepping method
+is applied.
+
+# Implementation
+
+## Folder layout and files
+
+```
+└── shallowwaterfriction/
+ ├── CMakeLists.txt -> build system file
+ ├── main.cc -> main program flow
+ ├── params.input -> runtime parameters
+ ├── properties.hh -> compile time configuration
+ ├── problem.hh -> boundary & initial conditions
+ └── spatialparams.hh -> spatial parameter fields
+```
diff --git a/examples/shallowwaterfriction/doc/intro.md b/examples/shallowwaterfriction/doc/intro.md
deleted file mode 100644
index f0f5e0a155dce73ff47df3b4d029de20963e7165..0000000000000000000000000000000000000000
--- a/examples/shallowwaterfriction/doc/intro.md
+++ /dev/null
@@ -1,97 +0,0 @@
-# Shallow water flow with bottom friction
-This example shows how the shallow water flow model can be
-applied to simulate steady subcritical flow including
-bottom friction (bed shear stress).
-
-__You will learn how to__
-
-* solve a shallow water flow problem including bottom friction
-* compute and output (VTK) an analytical reference solution
-
-__Result__. The numerical and analytical solutions for the problem will look like this:
-
-
-
-__Table of contents__. This description is structured as follows:
-
-[[_TOC_]]
-
-## Mathematical model
-The 2D shallow water equations (SWEs) are given by
-
-```math
-\frac{\partial \mathbf{U}}{\partial t} +
-\frac{\partial \mathbf{F}}{\partial x} +
-\frac{\partial \mathbf{G}}{\partial y} - \mathbf{S_b} - \mathbf{S_f} = 0
-```
-
-where $`\mathbf{U}`$, $`\mathbf{F}`$ and $`\mathbf{G}`$ defined as
-
-```math
-\mathbf{U} = \begin{bmatrix} h \\ uh \\ vh \end{bmatrix},
-\mathbf{F} = \begin{bmatrix} hu \\ hu^2 + \frac{1}{2} gh^2 \\ huv \end{bmatrix},
-\mathbf{G} = \begin{bmatrix} hv \\ huv \\ hv^2 + \frac{1}{2} gh^2 \end{bmatrix}
-```
-
-$`Z`$ is the bedSurface, $`h`$ the water depth, $`u`$ the velocity in
-x-direction and $`v`$ the velocity in y-direction, $`g`$ is the constant of gravity.
-
-The source terms for the bed friction $`\mathbf{S_b}`$ and bed slope
-$`\mathbf{S_f}`$ are given as
-
-```math
-\mathbf{S_b} = \begin{bmatrix} 0 \\ -gh \frac{\partial z}{\partial x}
- \\ -gh \frac{\partial z}{\partial y}\end{bmatrix},
-\mathbf{S_f} = \begin{bmatrix} 0 \\ -ghS_{fx} \\ -ghS_{fy}\end{bmatrix}.
-```
-
-For this example, a cell-centered finite volume method (`cctpfa`) is applied to solve the SWEs
-in combination with a fully-implicit time discretization. For cases where no sharp fronts or
-traveling waves occur it is possible to apply time steps larger than CFL number = 1 to reduce
-the computation time. Even if a steady state solution is considered, an implicit time stepping method
-is applied.
-
-## Problem set-up
-
-The model domain is given by a rough channel with a slope of 0.001.
-The domain is 500 meters long and 10 meters wide.
-.
-
-Bottom friction is considered by applying
-the friction law of Manning (Manning n = 0.025). At the lateral sides no friction is considered and a
-no-flow no slip boundary condition is applied. This is the default boundary condition for the shallow water model.
-
-At the left border a discharge boundary condition
-is applied as inflow boundary condition with q = -1.0 ($`m^2 s^{-1}`$). At the right border a water fixed depth boundary condition
-is applied for the outflow. Normal flow is assumed, therefore the water depth at the right border is calculated after
-the of Gaukler-Manning-Strickler equation:
-
-```math
-v_m = n^{-1} R_{hy}^{2/3} I_s^{1/2}
-```
-
-Where the mean velocity $`v_m`$ is given as $`v_m = \frac{q}{h}`$,
-$`n`$ is the friction value after Manning. $`R_{hy}`$ the hydraulic radius, which is assumed to be equal to
-the water depth. $`I_s`$ is the bed slope and $`q`$ the unity inflow discharge
-
-The water depth h can be calculated as
-```math
-h = \left(\frac{n q}{\sqrt{I_s}} \right)^{3/5}
-```
-
-The formula of Gaukler Manning and Strickler is also used to calculate the analytic solution. All parameters
-for the simulation are given in the file `params.input`.
-
-# Implementation
-
-## Folder layout and files
-
-```
-└── shallowwaterfriction/
- ├── CMakeLists.txt -> build system file
- ├── main.cc -> main program flow
- ├── params.input -> runtime parameters
- ├── properties.hh -> compile time configuration
- ├── problem.hh -> boundary & initial conditions
- └── spatialparams.hh -> spatial parameter fields
-```
diff --git a/examples/shallowwaterfriction/doc/main.md b/examples/shallowwaterfriction/doc/main.md
new file mode 100644
index 0000000000000000000000000000000000000000..52ad48457c0364d705e3415ad60f0c679285db99
--- /dev/null
+++ b/examples/shallowwaterfriction/doc/main.md
@@ -0,0 +1,301 @@
+<!-- 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](swe.md) |
+|---|---:|
+
+# Part 2: Main program flow
+
+We want to solve a shallow water flow problem in a rough
+channel to obtain the water table and
+compare it to the analytic solution. This is done with the
+`main()` function
+of the program which is defined in the file `main.cc` described below.
+
+The code documentation is structured as follows:
+
+[[_TOC_]]
+
+
+
+## The main file `main.cc`
+
+<details open>
+<summary><b>Click to hide/show the file documentation</b> (or inspect the [source code](../main.cc))</summary>
+
+
+### Included header files
+<details><summary> Click to show includes</summary>
+
+These are DUNE helper classes related to parallel computations, time measurements and file I/O
+
+```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>
+```
+
+The following headers include functionality related to property definition or retrieval, as well as
+the retrieval of input parameters specified in the input file or via the command line.
+
+```cpp
+#include <dumux/common/properties.hh>
+#include <dumux/common/parameters.hh>
+```
+
+The following files contains the available linear solver backends, the non linear Newton Solver
+and the assembler for the linear systems arising from finite volume discretizations
+(box-scheme, tpfa-approximation, mpfa-approximation).
+
+```cpp
+#include <dumux/linear/amgbackend.hh>
+#include <dumux/nonlinear/newtonsolver.hh>
+#include <dumux/assembly/fvassembler.hh>
+```
+
+The following class provides a convenient way of writing of dumux simulation results to VTK format.
+
+```cpp
+#include <dumux/io/vtkoutputmodule.hh>
+```
+
+The gridmanager constructs a grid from the information in the input or grid file.
+Many different Dune grid implementations are supported, of which a list can be found
+in `gridmanager.hh`.
+
+```cpp
+#include <dumux/io/grid/gridmanager.hh>
+```
+
+We include the header file specifing the properties of this example
+
+```cpp
+#include "properties.hh"
+```
+
+</details>
+
+### The main function
+We will now discuss the main program flow implemented within the `main` function.
+At the beginning of each program using Dune, an instance of `Dune::MPIHelper` has to
+be created. Moreover, we parse the run-time arguments from the command line and the
+input file:
+
+```cpp
+int main(int argc, char** argv) try
+{
+ using namespace Dumux;
+
+ // The Dune MPIHelper must be instantiated for each program using Dune
+ Dune::MPIHelper::instance(argc, argv);
+
+ // We parse command line arguments and input file
+ Parameters::init(argc, argv);
+```
+
+We define a convenience alias for the type tag of the problem. The type
+tag contains all the properties that are needed to define the model and the problem
+setup. Throughout the main file, we will obtain types defined for these type tag
+using the property system, i.e. with `GetPropType`.
+
+```cpp
+ using TypeTag = Properties::TTag::RoughChannel;
+```
+
+#### Step 1: Create the grid
+The `GridManager` class creates the grid from information given in the input file.
+This can either be a grid file or, in the case of structured grids, one can specify the coordinates
+of the corners of the grid and the number of cells to be used to discretize each spatial direction.
+
+```cpp
+ GridManager<GetPropType<TypeTag, Properties::Grid>> gridManager;
+ gridManager.init();
+
+ // We compute on the leaf grid view
+ const auto& leafGridView = gridManager.grid().leafGridView();
+```
+
+#### Step 2: Solving the shallow water problem
+First, a finite volume grid geometry is constructed from the grid that was created above.
+This builds the sub-control volumes (scv) and sub-control volume 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();
+```
+
+We now instantiate the problem, in which we define the boundary and initial conditions.
+
+```cpp
+ using Problem = GetPropType<TypeTag, Properties::Problem>;
+ auto problem = std::make_shared<Problem>(gridGeometry);
+```
+
+We initialize the solution vector. The shallow water problem is transient,
+therefore the initial solution defined in the problem is applied to the
+solution vector. On the basis of this solution, we initialize then the grid variables.
+
+```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);
+```
+
+Let us now instantiate the time loop. Therefore, we read in some time loop parameters from the input file.
+The parameter `tEnd` defines the duration of the simulation, `dt` the initial time step size and `maxDt` the maximal time step size.
+Moreover, we define the end of the simulation `tEnd` as check point in the time loop at which we will write the solution to vtk files.
+
+```cpp
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>; // type for scalar values
+ const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
+ const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
+ auto dt = getParam<Scalar>("TimeLoop.DtInitial");
+
+ // We instantiate time loop.
+ auto timeLoop = std::make_shared<CheckPointTimeLoop<Scalar>>(0.0, dt, tEnd);
+ timeLoop->setMaxTimeStepSize(maxDt);
+ timeLoop->setCheckPoint(tEnd);
+```
+
+We initialize the assembler with a time loop for the transient problem.
+Within the time loop, we will use this assembler in each time step to assemble the linear system.
+
+```cpp
+ using Assembler = FVAssembler<TypeTag, DiffMethod::numeric>;
+ auto assembler = std::make_shared<Assembler>(problem, gridGeometry, gridVariables, timeLoop, xOld);
+```
+
+We initialize the linear solver.
+
+```cpp
+ using LinearSolver = AMGBiCGSTABBackend<LinearSolverTraits<GridGeometry>>;
+ auto linearSolver = std::make_shared<LinearSolver>(leafGridView, gridGeometry->dofMapper());
+```
+
+We initialize the non-linear solver.
+
+```cpp
+ using NewtonSolver = Dumux::NewtonSolver<Assembler, LinearSolver>;
+ NewtonSolver nonLinearSolver(assembler, linearSolver);
+```
+
+The following lines of code initialize the vtk output module, add the velocity output facility
+and write out the initial solution. At each checkpoint, we will use the output module to write
+the solution of a time step into a corresponding vtk file.
+
+```cpp
+ VtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables,x, problem->name());
+
+ // add model-specific output fields to the writer
+ using IOFields = GetPropType<TypeTag, Properties::IOFields>;
+ IOFields::initOutputModule(vtkWriter);
+
+ // We add the analytical solution ("exactWaterDepth" and "exactVelocityX") to the predefined specific output.
+ vtkWriter.addField(problem->getExactWaterDepth(), "exactWaterDepth");
+ vtkWriter.addField(problem->getExactVelocityX(), "exactVelocityX");
+
+ // We calculate the analytic solution.
+ problem->analyticalSolution();
+
+ // write initial solution (including the above calculated analytical solution.
+ vtkWriter.write(0.0);
+```
+
+##### The time loop
+We start the time loop and solve a new time step as long as `tEnd` is not reached. In every time step,
+the problem is assembled and solved, the solution is updated, and when a checkpoint is reached the solution
+is written to a new vtk file. In addition, statistics related to CPU time, the current simulation time
+and the time step sizes used is printed to the terminal.
+
+```cpp
+ timeLoop->start(); do
+ {
+ // 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, using Newthon's method.
+ 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, if we reached the check point (end of time loop)
+ if (timeLoop->isCheckPoint())
+ vtkWriter.write(timeLoop->time());
+
+ // We report statistics of this time step.
+ timeLoop->reportTimeStep();
+
+ // We set new dt as suggested by newton controller 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.
+
+```cpp
+ timeLoop->finalize(leafGridView.comm());
+
+ return 0;
+}
+```
+
+#### 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 (const Dumux::ParameterException &e)
+{
+ std::cerr << std::endl << e << " ---> Abort!" << std::endl;
+ return 1;
+}
+// errors related to the parsing of Dune grid files
+catch (const 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;
+}
+// generic error handling with Dune::Exception
+catch (const Dune::Exception &e)
+{
+ std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl;
+ return 3;
+}
+// other exceptions
+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](swe.md) |
+|---|---:|
+
diff --git a/examples/shallowwaterfriction/doc/main_intro.md b/examples/shallowwaterfriction/doc/main_intro.md
new file mode 100644
index 0000000000000000000000000000000000000000..15c6b8515d544cf54a2ef9663e7732683fd49b9f
--- /dev/null
+++ b/examples/shallowwaterfriction/doc/main_intro.md
@@ -0,0 +1,12 @@
+# Part 2: Main program flow
+
+We want to solve a shallow water flow problem in a rough
+channel to obtain the water table and
+compare it to the analytic solution. This is done with the
+`main()` function
+of the program which is defined in the file `main.cc` described below.
+
+The code documentation is structured as follows:
+
+[[_TOC_]]
+
diff --git a/examples/shallowwaterfriction/doc/swe.md b/examples/shallowwaterfriction/doc/swe.md
new file mode 100644
index 0000000000000000000000000000000000000000..8a94e8440be20369989531bd712c9c333f49b649
--- /dev/null
+++ b/examples/shallowwaterfriction/doc/swe.md
@@ -0,0 +1,570 @@
+<!-- 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 shallow water flow simulation setup
+
+The shallow water flow setup, including the bottom friction,
+is implemented in the files `properties.hh`,
+`problem.hh` and `spatialparams.hh`. In the first of these files, a new
+type tag is declared for this problem. This allows the specialization
+of DuMu<sup>x</sup> `properties` for this type tag, which can be used to customize compile-time
+settings for the simulation. Two exemplary `properties`, that are mandatory to be
+specialized, are `Problem` and `SpatialParams`. With the first, one sets the
+`Problem` class to be used, in which users can define initial and boundary conditions.
+Similarly, in the `SpatialParams` class one implements the parameter distributions
+(e.g. friction value) that should be used by the model.
+
+The documentation provided in the sequel is structured as follows:
+
+[[_TOC_]]
+
+
+
+## Compile-time settings (`properties.hh`)
+
+In this file, the type tag used for the shallow water flow simulation is defined,
+for which we then specialize `properties` to the needs of the desired setup.
+
+
+<details open>
+<summary><b>Click to hide/show the file documentation</b> (or inspect the [source code](../properties.hh))</summary>
+
+
+### Includes
+<details><summary> Click to show include files</summary>
+
+The `ShallowWater` type tag specializes most of the `properties` required for a
+shallow water flow simulation in DuMu<sup>x</sup>. We will use this in the following to inherit the
+respective properties and subsequently specialize those `properties` for our
+type tag, which we want to modify or for which no meaningful default can be set.
+
+```cpp
+#include <dumux/freeflow/shallowwater/model.hh>
+```
+
+We want to use `YaspGrid`, an implementation of the dune grid interface for structured grids:
+
+```cpp
+#include <dune/grid/yaspgrid.hh>
+```
+
+In this example, we want to discretize the equations with the cell centered finite volume
+scheme using two-point-flux approximation:
+
+```cpp
+#include <dumux/discretization/cctpfa.hh>
+```
+
+We include the problem and spatial parameters headers used for this simulation.
+
+```cpp
+#include "problem.hh"
+#include "spatialparams.hh"
+```
+
+</details>
+
+### Type tag definition
+
+First, a so-called type tag is created. Properties are traits specialized for this type tag (a simple `struct`).
+The properties of two other type tags are inherited by adding the alias `InheritsFrom`.
+Here, properties from the shallow water model (`TTag::ShallowWater`) and the
+cell-centered finite volume scheme with two-point-flux approximation (`TTag::CCTpfaModel`)
+are inherited. These other type tag definitions can be found in the included
+headers `dumux/freeflow/shallowwater/model.hh` and `dumux/discretization/cctpfa.hh`.
+
+```cpp
+// We enter the namespace Dumux::Properties in order to import the entire Dumux namespace for general use:
+namespace Dumux::Properties {
+
+namespace TTag {
+struct RoughChannel { using InheritsFrom = std::tuple<ShallowWater, CCTpfaModel>; };
+}
+```
+
+### Property specializations
+
+We use a structured Cartesian grid with tensor product structure.
+`Dune::YaspGrid` (Yet Another Structure Parallel Grid) is defined in `dune/grid/yaspgrid.hh`
+in the Dune module `dune-grid`.
+
+```cpp
+template<class TypeTag>
+struct Grid<TypeTag, TTag::RoughChannel>
+{ using type = Dune::YaspGrid<2, Dune::TensorProductCoordinates<GetPropType<TypeTag, Properties::Scalar>, 2> >; };
+```
+
+Next, we specialize the properties `Problem` and `SpatialParams` for our new type tag 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::RoughChannel>
+{ using type = Dumux::RoughChannelProblem<TypeTag>; };
+
+template<class TypeTag>
+struct SpatialParams<TypeTag, TTag::RoughChannel>
+{
+ using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
+ using VolumeVariables = typename ElementVolumeVariables::VolumeVariables;
+
+ using type = RoughChannelSpatialParams<GridGeometry, Scalar, VolumeVariables>;
+};
+```
+
+Finally, we enable caching for the grid geometry. The cache
+stores values that were already calculated for later usage.
+This makes the simulation run faster but it uses more memory.
+
+```cpp
+template<class TypeTag>
+struct EnableGridGeometryCache<TypeTag, TTag::RoughChannel>
+{ static constexpr bool value = true; };
+
+} // end namespace Dumux::Properties
+```
+
+
+</details>
+
+
+
+## The problem file (`problem.hh`)
+
+This file contains the __problem class__ which defines the initial and boundary
+conditions for the shallow water flow simulation with bottom friction.
+In addition, the analytical solution is defined here.
+
+
+<details open>
+<summary><b>Click to hide/show the file documentation</b> (or inspect the [source code](../problem.hh))</summary>
+
+
+### Include files
+
+The first include we need here is the `ShallowWaterProblem` class, the base
+class from which we will derive.
+
+```cpp
+#include <dumux/freeflow/shallowwater/problem.hh>
+```
+
+In addition, we need the boundaryflux header, which handels the flux over
+the model boundaries.
+
+```cpp
+#include <dumux/freeflow/shallowwater/boundaryfluxes.hh>
+```
+
+### The problem class
+We enter the problem class where all necessary boundary conditions and initial conditions are set for our simulation.
+In addition the analytical solution of the problem is calculated.
+As this is a shallow water problem, we inherit from the basic ShallowWaterProblem.
+
+```cpp
+namespace Dumux {
+
+template <class TypeTag>
+class RoughChannelProblem : public ShallowWaterProblem<TypeTag>
+{
+ // A few convenience aliases used throughout this class.
+ using ParentType = ShallowWaterProblem<TypeTag>;
+ using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
+ using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
+ using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
+ using NeumannFluxes = GetPropType<TypeTag, Properties::NumEqVector>;
+ using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
+ using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
+ using ElementFluxVariablesCache = typename GridVariables::GridFluxVariablesCache::LocalView;
+ using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
+ using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
+ using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+ using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
+ using SubControlVolume = typename FVElementGeometry::SubControlVolume;
+
+public:
+ // This is the constructor of our problem class.
+ RoughChannelProblem(std::shared_ptr<const GridGeometry> gridGeometry)
+ : ParentType(gridGeometry)
+ {
+ // We read the parameters from the params.input file.
+ constManningN_ = getParam<Scalar>("Problem.ManningN");
+ bedSlope_ = getParam<Scalar>("Problem.BedSlope");
+ discharge_ = getParam<Scalar>("Problem.Discharge");
+ // We calculate the outflow boundary condition using the Gauckler-Manning-Strickler formula.
+ hBoundary_ = this->gaucklerManningStrickler(discharge_,constManningN_,bedSlope_);
+ // We initialize the analytic solution to a verctor of the appropriate size filled with zeros.
+ exactWaterDepth_.resize(gridGeometry->numDofs(), 0.0);
+ exactVelocityX_.resize(gridGeometry->numDofs(), 0.0);
+ }
+```
+
+#### Analytical Solution
+
+The analytical solution is calculated using the equation of Gauckler, Manning and Strickler.
+
+```cpp
+
+ // Equation of Gauckler, Manning and Strickler
+ Scalar gaucklerManningStrickler(Scalar discharge, Scalar manningN, Scalar bedSlope)
+ {
+ using std::pow;
+ using std::abs;
+ using std::sqrt;
+
+ return pow(abs(discharge)*manningN/sqrt(bedSlope), 0.6);
+ }
+
+ // Calculate the analytical solution
+ void analyticalSolution()
+ {
+ using std::abs;
+
+ for (const auto& element : elements(this->gridGeometry().gridView()))
+ {
+ const Scalar h = this->gaucklerManningStrickler(discharge_,constManningN_,bedSlope_);
+ const Scalar u = abs(discharge_)/h;
+
+ const auto eIdx = this->gridGeometry().elementMapper().index(element);
+ exactWaterDepth_[eIdx] = h;
+ exactVelocityX_[eIdx] = u;
+ }
+ }
+
+ // Getter function for the analytical solution of the water depth
+ const std::vector<Scalar>& getExactWaterDepth()
+ {
+ return exactWaterDepth_;
+ }
+
+ // Getter function for the analytical solution of the velocity in x-direction
+ const std::vector<Scalar>& getExactVelocityX()
+ {
+ return exactVelocityX_;
+ }
+```
+
+#### Bottom friction
+
+The bottom friction is a source term and therefore handled by the `source` function.
+
+```cpp
+ NumEqVector source(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolume &scv) const
+ {
+
+ NumEqVector source (0.0);
+
+ // Since the bed slope source term is handels within the flux computation,
+ // in this model the bottom friction is the only source term.
+ source += bottomFrictionSource(element, fvGeometry, elemVolVars, scv);
+
+ return source;
+ }
+```
+
+The calculation of the source term due to bottom friction needs the bottom shear stess.
+This is the force per area, which works between the flow and the channel bed
+(1D vector with two entries) and is calculated within the `FrictionLaw` class.
+The bottom friction causes a loss of momentum. Thus the first entry of the `bottomFrictionSource`,
+which is related to the mass balance equation is zero.
+The second entry of the `bottomFricitonSource` corresponds to the momentum equation in x-direction
+and is therefore equal to the first, the x-component, of the `bottomShearStress`.
+Accordingly, the third entry of the `bottomFrictionSource` is equal to the second component of the `bottomShearStress`.
+
+```cpp
+ NumEqVector bottomFrictionSource(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolume &scv) const
+ {
+ NumEqVector bottomFrictionSource(0.0);
+ const auto& volVars = elemVolVars[scv];
+
+ // bottom shear stress vector
+ Dune::FieldVector<Scalar, 2> bottomShearStress = this->spatialParams().frictionLaw(element, scv).shearStress(volVars);
+
+ // source term due to bottom friction
+ bottomFrictionSource[0] = 0.0;
+ bottomFrictionSource[1] = bottomShearStress[0];
+ bottomFrictionSource[2] = bottomShearStress[1];
+
+ return bottomFrictionSource;
+ }
+```
+
+#### Boundary conditions
+
+We define the __type of all boundary conditions__ as neumann-type,
+because we use a weak imposition.
+
+```cpp
+ BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
+ {
+ BoundaryTypes bcTypes;
+ bcTypes.setAllNeumann();
+ return bcTypes;
+ }
+```
+
+In the following function we implement the __Neumann boundary conditions__.
+Due to the weak imposition we calculate the flux at the boundary with a Riemann solver.
+This needs the state of a virtual cell outside of the boundary (`boundaryStateVariables`),
+wich is calculated with the Riemann invariants
+(see: Yoon and Kang, "Finite Volume Model for Two-Dimensional Shallow Water Flows on Unstructured Grids").
+
+```cpp
+ NeumannFluxes neumann(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const ElementFluxVariablesCache& elemFluxVarsCache,
+ const SubControlVolumeFace& scvf) const
+ {
+ NeumannFluxes values(0.0);
+
+ const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
+ const auto& insideVolVars = elemVolVars[insideScv];
+ const auto& nxy = scvf.unitOuterNormal();
+ const auto gravity = this->spatialParams().gravity(scvf.center());
+ std::array<Scalar, 3> boundaryStateVariables;
+
+ // Calculate the Riemann invariants for imposed discharge at the left side.
+ if (scvf.center()[0] < 0.0 + eps_)
+ {
+ boundaryStateVariables = ShallowWater::fixedDischargeBoundary(discharge_,
+ insideVolVars.waterDepth(),
+ insideVolVars.velocity(0),
+ insideVolVars.velocity(1),
+ gravity,
+ nxy);
+ }
+ // Calculate the Riemann invariants for imposed water depth at the right side.
+ else if (scvf.center()[0] > 500.0 - eps_)
+ {
+ boundaryStateVariables = ShallowWater::fixedWaterDepthBoundary(hBoundary_,
+ insideVolVars.waterDepth(),
+ insideVolVars.velocity(0),
+ insideVolVars.velocity(1),
+ gravity,
+ nxy);
+ }
+ // Calculate the Riemann invariants for the no-flow boundary.
+ else
+ {
+ boundaryStateVariables[0] = insideVolVars.waterDepth();
+ boundaryStateVariables[1] = -insideVolVars.velocity(0);
+ boundaryStateVariables[2] = -insideVolVars.velocity(1);
+ }
+ // Calculate the boundary fluxes based on a Riemann problem.
+ auto riemannFlux = ShallowWater::riemannProblem(insideVolVars.waterDepth(),
+ boundaryStateVariables[0],
+ insideVolVars.velocity(0),
+ boundaryStateVariables[1],
+ insideVolVars.velocity(1),
+ boundaryStateVariables[2],
+ insideVolVars.bedSurface(),
+ insideVolVars.bedSurface(),
+ gravity,
+ nxy);
+
+ values[Indices::massBalanceIdx] = riemannFlux[0];
+ values[Indices::velocityXIdx] = riemannFlux[1];
+ values[Indices::velocityYIdx] = riemannFlux[2];
+
+ return values;
+ }
+```
+
+#### Initial conditions
+
+We specify the initial conditions for the primary variables (water depth, velocity in y-direction
+and velocity in x-direction). In this example constant initial conditions are used. Therefore the
+argument `globalPos` is not needed. If you want to impose spatial variable initial conditions,
+you have to use the `globalPos` argument.
+
+```cpp
+ PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
+ {
+ PrimaryVariables initialValues(0.0);
+ // We set the initial water depth to one meter.
+ initialValues[0] = 1.0;
+ // We set the x-component of the initial velocity to zero.
+ initialValues[1] = 0.0;
+ // We set the y-component of the initial velocity to zero.
+ initialValues[2] = 0.0;
+
+ return initialValues;
+ }
+```
+
+We declare the private variables of the problem.
+
+```cpp
+private:
+ // variables for the analytic solution.
+ std::vector<Scalar> exactWaterDepth_;
+ std::vector<Scalar> exactVelocityX_;
+ // constant friction value (an analytic solution is only available for const friction).
+ Scalar constManningN_;
+ // The constant channel bed slope.
+ Scalar bedSlope_;
+ // The water depth at the outflow boundary.
+ Scalar hBoundary_;
+ // The discharge at the inflow boundary.
+ Scalar discharge_;
+ // We assign a private global variable for the epsilon:
+ static constexpr Scalar eps_ = 1.0e-6;
+
+}; // end class definition RoughChannelProblem
+} // end namespace Dumux
+```
+
+
+</details>
+
+
+
+## Parameter distributions (`spatialparams.hh`)
+
+This file contains the __spatial parameters class__ which defines the
+the friction law, including it's friction parameter, the acceleration
+due to gravity and the altitude of the channel bed surface. In this example only the bed
+surface has a non constant distribution.
+
+
+<details open>
+<summary><b>Click to hide/show the file documentation</b> (or inspect the [source code](../spatialparams.hh))</summary>
+
+
+### Include files
+We include the basic spatial parameters file for finite volumes, from which we will inherit.
+
+```cpp
+#include <dumux/material/spatialparams/fv.hh>
+```
+
+We include all friction laws.
+
+```cpp
+#include <dumux/material/fluidmatrixinteractions/frictionlaws/frictionlaw.hh>
+#include <dumux/material/fluidmatrixinteractions/frictionlaws/manning.hh>
+#include <dumux/material/fluidmatrixinteractions/frictionlaws/nikuradse.hh>
+#include <dumux/material/fluidmatrixinteractions/frictionlaws/nofriction.hh>
+```
+
+### The spatial parameters class
+
+In the `RoughChannelSpatialParams` class, we define all functions needed to describe
+the rough channel for the shallow water problem.
+We inherit from the `FVSpatialParams` class, which is the base class
+for spatial parameters in the context of
+applications using finite volume discretization schemes.
+
+```cpp
+namespace Dumux {
+
+template<class GridGeometry, class Scalar, class VolumeVariables>
+class RoughChannelSpatialParams
+: public FVSpatialParams<GridGeometry, Scalar,
+ RoughChannelSpatialParams<GridGeometry, Scalar, VolumeVariables>>
+{
+ // This convenience aliases will be used throughout this class
+ using ThisType = RoughChannelSpatialParams<GridGeometry, Scalar, VolumeVariables>;
+ using ParentType = FVSpatialParams<GridGeometry, Scalar, ThisType>;
+ using GridView = typename GridGeometry::GridView;
+ using FVElementGeometry = typename GridGeometry::LocalView;
+ using SubControlVolume = typename FVElementGeometry::SubControlVolume;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+```
+
+In the following, the properties of the the rough channel are set. Namely, these are
+the friction law, including it's friction parameter, the acceleration
+due to gravity and the altitude of the channel bed surface.
+
+```cpp
+public:
+ // In the constructor we read some values from the `params.input` and initialize the friciton law.
+ RoughChannelSpatialParams(std::shared_ptr<const GridGeometry> gridGeometry)
+ : ParentType(gridGeometry)
+ {
+ gravity_ = getParam<Scalar>("Problem.Gravity");
+ bedSlope_ = getParam<Scalar>("Problem.BedSlope");
+ frictionLawType_ = getParam<std::string>("Problem.FrictionLaw");
+ initFrictionLaw();
+ }
+
+ // This function handles the initialization of the friction law based on the settings
+ // specified in `params.input`.
+ void initFrictionLaw()
+ {
+ if (frictionLawType_ == "Manning")
+ {
+ Scalar manningN = getParam<Scalar>("Problem.ManningN");
+ frictionLaw_ = std::make_unique<FrictionLawManning<VolumeVariables>>(gravity_, manningN);
+ }
+ else if (frictionLawType_ == "Nikuradse")
+ {
+ Scalar ks = getParam<Scalar>("Problem.Ks");
+ frictionLaw_ = std::make_unique<FrictionLawNikuradse<VolumeVariables>>(ks);
+ }
+ else if (frictionLawType_ == "None")
+ {
+ frictionLaw_ = std::make_unique<FrictionLawNoFriction<VolumeVariables>>();
+ }
+ else
+ {
+ std::cout<<"The FrictionLaw in params.input is unknown. Valid entries are `Manning`,"
+ " `Nikuradse` and `None`!"<<std::endl;
+ }
+ }
+
+ // This function returns an object of the friction law class, already initialized with a friction value.
+ const FrictionLaw<VolumeVariables>& frictionLaw(const Element& element,
+ const SubControlVolume& scv) const
+ {
+ return *frictionLaw_;
+ }
+
+ // This function returns the acceleration due to gravity.
+ Scalar gravity(const GlobalPosition& globalPos) const
+ {
+ return gravity_;
+ }
+
+ // Define the bed surface based on the bed slope and the bed level at the inflow (10 m).
+ Scalar bedSurface(const Element& element,
+ const SubControlVolume& scv) const
+ {
+ return 10.0 - element.geometry().center()[0] * bedSlope_;
+ }
+
+// We declare the private variables of the problem.
+private:
+ Scalar gravity_;
+ Scalar bedSlope_;
+ std::string frictionLawType_;
+ std::unique_ptr<FrictionLaw<VolumeVariables>> frictionLaw_;
+}; // end class definition of RoughChannelSpatialParams
+} // end of namespace Dumux.
+```
+
+
+</details>
+
+
+| [:arrow_left: Back to the main documentation](../README.md) | [:arrow_right: Continue with part 2](main.md) |
+|---|---:|
+
diff --git a/examples/shallowwaterfriction/doc/swe_intro.md b/examples/shallowwaterfriction/doc/swe_intro.md
new file mode 100644
index 0000000000000000000000000000000000000000..c0e5511287903793cfc4f4954cb9186e948e04a2
--- /dev/null
+++ b/examples/shallowwaterfriction/doc/swe_intro.md
@@ -0,0 +1,17 @@
+# Part 1: Implementation of the shallow water flow simulation setup
+
+The shallow water flow setup, including the bottom friction,
+is implemented in the files `properties.hh`,
+`problem.hh` and `spatialparams.hh`. In the first of these files, a new
+type tag is declared for this problem. This allows the specialization
+of DuMu<sup>x</sup> `properties` for this type tag, which can be used to customize compile-time
+settings for the simulation. Two exemplary `properties`, that are mandatory to be
+specialized, are `Problem` and `SpatialParams`. With the first, one sets the
+`Problem` class to be used, in which users can define initial and boundary conditions.
+Similarly, in the `SpatialParams` class one implements the parameter distributions
+(e.g. friction value) that should be used by the model.
+
+The documentation provided in the sequel is structured as follows:
+
+[[_TOC_]]
+
diff --git a/examples/shallowwaterfriction/img/domain.png b/examples/shallowwaterfriction/img/domain.png
deleted file mode 100644
index 27029ab28e1de01c7736fb54be698db508f5777f..0000000000000000000000000000000000000000
Binary files a/examples/shallowwaterfriction/img/domain.png and /dev/null differ
diff --git a/examples/shallowwaterfriction/img/result.png b/examples/shallowwaterfriction/img/result.png
deleted file mode 100644
index 1ef7c2c9798930cfd4d6b1b8a98b78e0cb718406..0000000000000000000000000000000000000000
Binary files a/examples/shallowwaterfriction/img/result.png and /dev/null differ
diff --git a/examples/shallowwaterfriction/img/swe_result.png b/examples/shallowwaterfriction/img/swe_result.png
new file mode 100644
index 0000000000000000000000000000000000000000..0c5c8025ffbf8d476d146920508c563e1cb4a951
Binary files /dev/null and b/examples/shallowwaterfriction/img/swe_result.png differ
diff --git a/examples/shallowwaterfriction/main.cc b/examples/shallowwaterfriction/main.cc
index e3092b0eeccecb03cf1d6446e7f67db32063e61d..eba1b871dafe52c578f9f6acda27549ba589c4d4 100644
--- a/examples/shallowwaterfriction/main.cc
+++ b/examples/shallowwaterfriction/main.cc
@@ -17,138 +17,168 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
-// ## The file `main.cc`
+// ## The main file `main.cc`
+// [[content]]
//
-//
-// This is the main file for the shallow water example. Here we can see the programme sequence and how the system is solved using newton's method.
-// ### Includes
+// ### Included header files
+// [[details]] includes
+// [[exclude]]
+// Some generic includes.
#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.
+// [[/exclude]]
+//
+// These are DUNE helper classes related to parallel computations, time measurements and file I/O
#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>
-// We need the following class to simplify the writing of dumux simulation data to VTK format.
-#include <dumux/io/vtkoutputmodule.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.
+// The following headers include functionality related to property definition or retrieval, as well as
+// the retrieval of input parameters specified in the input file or via 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>
-// The gridmanager constructs a grid from the information in the input or grid file. There is a specification for the different supported grid managers.
-#include <dumux/io/grid/gridmanager.hh>
-// We include the linear solver to be used to solve the linear system
+
+// The following files contains the available linear solver backends, the non linear Newton Solver
+// and the assembler for the linear systems arising from finite volume discretizations
+// (box-scheme, tpfa-approximation, mpfa-approximation).
#include <dumux/linear/amgbackend.hh>
-#include <dumux/linear/linearsolvertraits.hh>
-// We include the nonlinear newtons 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)
#include <dumux/assembly/fvassembler.hh>
-// We include the properties
+
+// The following class provides a convenient way of writing of dumux simulation results to VTK format.
+#include <dumux/io/vtkoutputmodule.hh>
+// The gridmanager constructs a grid from the information in the input or grid file.
+// Many different Dune grid implementations are supported, of which a list can be found
+// in `gridmanager.hh`.
+#include <dumux/io/grid/gridmanager.hh>
+
+// We include the header file specifing the properties of this example
#include "properties.hh"
+// [[/details]]
+//
-// ### Beginning of the main function
+// ### The main function
+// We will now discuss the main program flow implemented within the `main` function.
+// At the beginning of each program using Dune, an instance of `Dune::MPIHelper` has to
+// be created. Moreover, we parse the run-time arguments from the command line and the
+// input file:
+// [[codeblock]]
int main(int argc, char** argv) try
{
using namespace Dumux;
- // We define the type tag for this problem
- using TypeTag = Properties::TTag::RoughChannel;
-
- // 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);
+ // The Dune MPIHelper must be instantiated for each program using Dune
+ Dune::MPIHelper::instance(argc, argv);
// We parse command line arguments and input file
Parameters::init(argc, argv);
+ // [[/codeblock]]
+
+ // We define a convenience alias for the type tag of the problem. The type
+ // tag contains all the properties that are needed to define the model and the problem
+ // setup. Throughout the main file, we will obtain types defined for these type tag
+ // using the property system, i.e. with `GetPropType`.
+ using TypeTag = Properties::TTag::RoughChannel;
- // ### Create the grid
- // A gridmanager tries to create the grid either from a grid file or the input file.
+ // #### Step 1: Create the grid
+ // The `GridManager` class creates the grid from information given in the input file.
+ // This can either be a grid file or, in the case of structured grids, one can specify the coordinates
+ // of the corners of the grid and the number of cells to be used to discretize each spatial direction.
+ //[[codeblock]]
GridManager<GetPropType<TypeTag, Properties::Grid>> gridManager;
gridManager.init();
// We compute on the leaf grid view
const auto& leafGridView = gridManager.grid().leafGridView();
+ //[[/codeblock]]
- // ### Setup and solving of the problem
- // #### Setup
- // 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.
+ // #### Step 2: Solving the shallow water problem
+ // First, a finite volume grid geometry is constructed from the grid that was created above.
+ // This builds the sub-control volumes (scv) and sub-control volume faces (scvf) for each element
+ // of the grid partition.
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.
+ // We now instantiate the problem, in which we define the boundary and initial conditions.
using Problem = GetPropType<TypeTag, Properties::Problem>;
auto problem = std::make_shared<Problem>(gridGeometry);
- // We initialize the solution vector
+ // We initialize the solution vector. The shallow water problem is transient,
+ // therefore the initial solution defined in the problem is applied to the
+ // solution vector. On the basis of this solution, we initialize then the grid variables.
using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
SolutionVector x(gridGeometry->numDofs());
problem->applyInitialSolution(x);
auto xOld = x;
- // And then use the solutionvector to intialize the gridVariables.
using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
auto gridVariables = std::make_shared<GridVariables>(problem, gridGeometry);
gridVariables->init(x);
- // We get some time loop parameters from the input file.
- using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ // Let us now instantiate the time loop. Therefore, we read in some time loop parameters from the input file.
+ // The parameter `tEnd` defines the duration of the simulation, `dt` the initial time step size and `maxDt` the maximal time step size.
+ // Moreover, we define the end of the simulation `tEnd` as check point in the time loop at which we will write the solution to vtk files.
+ // [[codeblock]]
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>; // type for scalar values
const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
auto dt = getParam<Scalar>("TimeLoop.DtInitial");
- // We intialize the vtk output module. 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());
- // We add the analytical solution ("exactWaterDepth" and "exactVelocityX") to the predefined specific output.
- vtkWriter.addField(problem->getExactWaterDepth(), "exactWaterDepth");
- vtkWriter.addField(problem->getExactVelocityX(), "exactVelocityX");
- // We calculate the analytic solution.
- problem->analyticalSolution();
- IOFields::initOutputModule(vtkWriter);
- vtkWriter.write(0.0);
-
// We instantiate time loop.
- auto timeLoop = std::make_shared<CheckPointTimeLoop<Scalar>>(0, dt, tEnd);
+ auto timeLoop = std::make_shared<CheckPointTimeLoop<Scalar>>(0.0, dt, tEnd);
timeLoop->setMaxTimeStepSize(maxDt);
+ timeLoop->setCheckPoint(tEnd);
+ // [[/codeblock]]
- //we set the assembler with the time loop because we have an instationary problem.
+ // We initialize the assembler with a time loop for the transient problem.
+ // Within the time loop, we will use this assembler in each time step to assemble the linear system.
using Assembler = FVAssembler<TypeTag, DiffMethod::numeric>;
auto assembler = std::make_shared<Assembler>(problem, gridGeometry, gridVariables, timeLoop, xOld);
- // We set the linear solver.
+ // We initialize the linear solver.
using LinearSolver = AMGBiCGSTABBackend<LinearSolverTraits<GridGeometry>>;
auto linearSolver = std::make_shared<LinearSolver>(leafGridView, gridGeometry->dofMapper());
- // Additionaly, we set the non-linear solver.
+ // We initialize the non-linear solver.
using NewtonSolver = Dumux::NewtonSolver<Assembler, LinearSolver>;
NewtonSolver nonLinearSolver(assembler, linearSolver);
- // We set some check point at the end of the time loop. The check point is used to trigger the vtk output.
- timeLoop->setCheckPoint(tEnd);
+ // The following lines of code initialize the vtk output module, add the velocity output facility
+ // and write out the initial solution. At each checkpoint, we will use the output module to write
+ // the solution of a time step into a corresponding vtk file.
+ // [[codeblock]]
+ VtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables,x, problem->name());
+
+ // add model-specific output fields to the writer
+ using IOFields = GetPropType<TypeTag, Properties::IOFields>;
+ IOFields::initOutputModule(vtkWriter);
- // We start the time loop.
+ // We add the analytical solution ("exactWaterDepth" and "exactVelocityX") to the predefined specific output.
+ vtkWriter.addField(problem->getExactWaterDepth(), "exactWaterDepth");
+ vtkWriter.addField(problem->getExactVelocityX(), "exactVelocityX");
+
+ // We calculate the analytic solution.
+ problem->analyticalSolution();
+
+ // write initial solution (including the above calculated analytical solution.
+ vtkWriter.write(0.0);
+ // [[/codeblock]]
+
+ // ##### The time loop
+ // We start the time loop and solve a new time step as long as `tEnd` is not reached. In every time step,
+ // the problem is assembled and solved, the solution is updated, and when a checkpoint is reached the solution
+ // is written to a new vtk file. In addition, statistics related to CPU time, the current simulation time
+ // and the time step sizes used is printed to the terminal.
+ // [[codeblock]]
timeLoop->start(); do
{
- // 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.
+ // 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.
+ // We solve the non-linear system with time step control, using Newthon's method.
nonLinearSolver.solve(x,*timeLoop);
// We make the new solution the old solution.
@@ -170,25 +200,26 @@ int main(int argc, char** argv) try
} while (!timeLoop->finished());
+ // [[/codeblock]]
+ // The following piece of code prints a final status report of the time loop
+ // before the program is terminated.
timeLoop->finalize(leafGridView.comm());
- // ### Final Output
- // We 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]] exception handler
+// [[codeblock]]
+// errors related to run-time parameters
catch (const Dumux::ParameterException &e)
{
std::cerr << std::endl << e << " ---> Abort!" << std::endl;
return 1;
}
+// errors related to the parsing of Dune grid files
catch (const Dune::DGFException & e)
{
std::cerr << "DGF exception thrown (" << e <<
@@ -198,13 +229,18 @@ catch (const Dune::DGFException & e)
<< " ---> Abort!" << std::endl;
return 2;
}
+// generic error handling with Dune::Exception
catch (const Dune::Exception &e)
{
std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl;
return 3;
}
+// other exceptions
catch (...)
{
std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl;
return 4;
}
+// [[/codeblock]]
+// [[/details]]
+// [[/content]]
diff --git a/examples/shallowwaterfriction/problem.hh b/examples/shallowwaterfriction/problem.hh
index 12bfbe3261223c699bf04878d28ca3c5ead26a86..caa370dee7cd5088684602d9fcadd0c1665f1736 100644
--- a/examples/shallowwaterfriction/problem.hh
+++ b/examples/shallowwaterfriction/problem.hh
@@ -20,21 +20,34 @@
#ifndef DUMUX_EXAMPLE_SHALLOWWATER_FRICTION_PROBLEM_HH
#define DUMUX_EXAMPLE_SHALLOWWATER_FRICTION_PROBLEM_HH
-// ## The file `problem.hh`
-// We start with includes
-#include <dumux/common/parameters.hh>
-#include <dumux/common/properties.hh>
+// ## The problem file (`problem.hh`)
+//
+// This file contains the __problem class__ which defines the initial and boundary
+// conditions for the shallow water flow simulation with bottom friction.
+// In addition, the analytical solution is defined here.
+//
+// [[content]]
+//
+// ### Include files
+//
+// The first include we need here is the `ShallowWaterProblem` class, the base
+// class from which we will derive.
#include <dumux/freeflow/shallowwater/problem.hh>
+// In addition, we need the boundaryflux header, which handels the flux over
+// the model boundaries.
#include <dumux/freeflow/shallowwater/boundaryfluxes.hh>
+// ### The problem class
// We enter the problem class where all necessary boundary conditions and initial conditions are set for our simulation.
+// In addition the analytical solution of the problem is calculated.
// As this is a shallow water problem, we inherit from the basic ShallowWaterProblem.
+// [[codeblock]]
namespace Dumux {
template <class TypeTag>
class RoughChannelProblem : public ShallowWaterProblem<TypeTag>
{
- // We use convenient declarations that we derive from the property system.
+ // A few convenience aliases used throughout this class.
using ParentType = ShallowWaterProblem<TypeTag>;
using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
@@ -45,7 +58,6 @@ class RoughChannelProblem : public ShallowWaterProblem<TypeTag>
using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
using ElementFluxVariablesCache = typename GridVariables::GridFluxVariablesCache::LocalView;
- using VolumeVariables = typename ElementVolumeVariables::VolumeVariables;
using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
@@ -60,31 +72,24 @@ public:
: ParentType(gridGeometry)
{
// We read the parameters from the params.input file.
- name_ = getParam<std::string>("Problem.Name");
constManningN_ = getParam<Scalar>("Problem.ManningN");
bedSlope_ = getParam<Scalar>("Problem.BedSlope");
discharge_ = getParam<Scalar>("Problem.Discharge");
// We calculate the outflow boundary condition using the Gauckler-Manning-Strickler formula.
- hBoundary_ = this->gauklerManningStrickler(discharge_,constManningN_,bedSlope_);
+ hBoundary_ = this->gaucklerManningStrickler(discharge_,constManningN_,bedSlope_);
// We initialize the analytic solution to a verctor of the appropriate size filled with zeros.
exactWaterDepth_.resize(gridGeometry->numDofs(), 0.0);
exactVelocityX_.resize(gridGeometry->numDofs(), 0.0);
}
+ // [[/codeblock]]
- // Get the analytical water depth
- const std::vector<Scalar>& getExactWaterDepth()
- {
- return exactWaterDepth_;
- }
-
- // Get the analytical velocity
- const std::vector<Scalar>& getExactVelocityX()
- {
- return exactVelocityX_;
- }
+ // #### Analytical Solution
+ //
+ // The analytical solution is calculated using the equation of Gauckler, Manning and Strickler.
+ // [[codeblock]]
- // Get the water depth with Gauckler-Manning-Strickler
- Scalar gauklerManningStrickler(Scalar discharge, Scalar manningN, Scalar bedSlope)
+ // Equation of Gauckler, Manning and Strickler
+ Scalar gaucklerManningStrickler(Scalar discharge, Scalar manningN, Scalar bedSlope)
{
using std::pow;
using std::abs;
@@ -93,14 +98,14 @@ public:
return pow(abs(discharge)*manningN/sqrt(bedSlope), 0.6);
}
- // Get the analytical solution
+ // Calculate the analytical solution
void analyticalSolution()
{
using std::abs;
for (const auto& element : elements(this->gridGeometry().gridView()))
{
- const Scalar h = this->gauklerManningStrickler(discharge_,constManningN_,bedSlope_);
+ const Scalar h = this->gaucklerManningStrickler(discharge_,constManningN_,bedSlope_);
const Scalar u = abs(discharge_)/h;
const auto eIdx = this->gridGeometry().elementMapper().index(element);
@@ -109,61 +114,87 @@ public:
}
}
- // Get the problem name. It is used as a prefix for files generated by the simulation.
- const std::string& name() const
+ // Getter function for the analytical solution of the water depth
+ const std::vector<Scalar>& getExactWaterDepth()
{
- return name_;
+ return exactWaterDepth_;
}
- // Get the source term.
- NumEqVector source(const Element& element,
- const FVElementGeometry& fvGeometry,
- const ElementVolumeVariables& elemVolVars,
- const SubControlVolume &scv) const
+ // Getter function for the analytical solution of the velocity in x-direction
+ const std::vector<Scalar>& getExactVelocityX()
+ {
+ return exactVelocityX_;
+ }
+ // [[/codeblock]]
+
+ // #### Bottom friction
+ //
+ // The bottom friction is a source term and therefore handled by the `source` function.
+ // [[codeblock]]
+ NumEqVector source(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolume &scv) const
{
NumEqVector source (0.0);
- // In this model the bottom friction is the only source.
+ // Since the bed slope source term is handels within the flux computation,
+ // in this model the bottom friction is the only source term.
source += bottomFrictionSource(element, fvGeometry, elemVolVars, scv);
return source;
}
-
- // Get the source term due to bottom friction.
- NumEqVector bottomFrictionSource(const Element& element,
+ // [[/codeblock]]
+
+ // The calculation of the source term due to bottom friction needs the bottom shear stess.
+ // This is the force per area, which works between the flow and the channel bed
+ // (1D vector with two entries) and is calculated within the `FrictionLaw` class.
+ // The bottom friction causes a loss of momentum. Thus the first entry of the `bottomFrictionSource`,
+ // which is related to the mass balance equation is zero.
+ // The second entry of the `bottomFricitonSource` corresponds to the momentum equation in x-direction
+ // and is therefore equal to the first, the x-component, of the `bottomShearStress`.
+ // Accordingly, the third entry of the `bottomFrictionSource` is equal to the second component of the `bottomShearStress`.
+ // [[codeblock]]
+ NumEqVector bottomFrictionSource(const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const SubControlVolume &scv) const
- {
+ {
NumEqVector bottomFrictionSource(0.0);
const auto& volVars = elemVolVars[scv];
- // For the calculation of the source term due to bottom friction the two-dimensional bottom shear stess vector is needed. This is the force per area, which works between the flow and the bed. It is calculated within the `FrictionLaw`, which is a spatialParameter. In this model the `FrictionLawManning` is used (see `params.input`).
+ // bottom shear stress vector
Dune::FieldVector<Scalar, 2> bottomShearStress = this->spatialParams().frictionLaw(element, scv).shearStress(volVars);
- // The bottom shear stress causes a pure loss of momentum. Thus the first entry of the `bottomFrictionSource`, which is related to the mass balance equation is zero. The second entry of the `bottomFricitonSource` corresponds to the momentum equation in x-direction and is therefore equal to the first, the x-component, of the `bottomShearStress`. Accordingly the third entry of the `bottomFrictionSource` is equal to the second component of the `bottomShearStress`.
+ // source term due to bottom friction
bottomFrictionSource[0] = 0.0;
bottomFrictionSource[1] = bottomShearStress[0];
bottomFrictionSource[2] = bottomShearStress[1];
return bottomFrictionSource;
- }
+ }
+ // [[/codeblock]]
- // We specify the boundary condition type.
+ // #### Boundary conditions
+ //
+ // We define the __type of all boundary conditions__ as neumann-type,
+ // because we use a weak imposition.
+ // [[codeblock]]
BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
{
BoundaryTypes bcTypes;
- // Since we use a weak imposition all boundary conditions are of Neumann type.
bcTypes.setAllNeumann();
return bcTypes;
}
-
- // We specify the neumann boundary. Due to the weak imposition we calculate the flux at the
- // boundary, with a Rieman solver. For this the state of a virtual cell outside of the boundary
- // is needed (`boundaryStateVariables`), wich is calculated with the Riemann invariants
- // (see Yoon and Kang, Finite Volume Model for Two-Dimensional Shallow Water Flows on Unstructured Grids).
- // The calculation of the Riemann invariants differ depending on the type of the boundary (h, q or no-flow boundary).
+ // [[/codeblock]]
+
+ // In the following function we implement the __Neumann boundary conditions__.
+ // Due to the weak imposition we calculate the flux at the boundary with a Riemann solver.
+ // This needs the state of a virtual cell outside of the boundary (`boundaryStateVariables`),
+ // wich is calculated with the Riemann invariants
+ // (see: Yoon and Kang, "Finite Volume Model for Two-Dimensional Shallow Water Flows on Unstructured Grids").
+ // [[codeblock]]
NeumannFluxes neumann(const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
@@ -178,7 +209,7 @@ public:
const auto gravity = this->spatialParams().gravity(scvf.center());
std::array<Scalar, 3> boundaryStateVariables;
- // Calculate the rieman invariants for imposed discharge at the left side.
+ // Calculate the Riemann invariants for imposed discharge at the left side.
if (scvf.center()[0] < 0.0 + eps_)
{
boundaryStateVariables = ShallowWater::fixedDischargeBoundary(discharge_,
@@ -188,8 +219,8 @@ public:
gravity,
nxy);
}
- // Calculate the rieman invariants for impose water depth at the right side.
- else if (scvf.center()[0] > 100.0 - eps_)
+ // Calculate the Riemann invariants for imposed water depth at the right side.
+ else if (scvf.center()[0] > 500.0 - eps_)
{
boundaryStateVariables = ShallowWater::fixedWaterDepthBoundary(hBoundary_,
insideVolVars.waterDepth(),
@@ -198,14 +229,14 @@ public:
gravity,
nxy);
}
- // Calculate the rieman invarianty for the no-flow boundary.
+ // Calculate the Riemann invariants for the no-flow boundary.
else
{
boundaryStateVariables[0] = insideVolVars.waterDepth();
boundaryStateVariables[1] = -insideVolVars.velocity(0);
boundaryStateVariables[2] = -insideVolVars.velocity(1);
}
- // We calculate the boundary fluxes based on a Riemann problem.
+ // Calculate the boundary fluxes based on a Riemann problem.
auto riemannFlux = ShallowWater::riemannProblem(insideVolVars.waterDepth(),
boundaryStateVariables[0],
insideVolVars.velocity(0),
@@ -223,41 +254,48 @@ public:
return values;
}
-
- // We set the initial conditions. In this example constant initial conditions are used. Therefore the argument `globalPos` is not needed. If you want to impose spatial variable initial conditions, you have to use the `globalPos`.
+ // [[/codeblock]]
+
+ // #### Initial conditions
+ //
+ // We specify the initial conditions for the primary variables (water depth, velocity in y-direction
+ // and velocity in x-direction). In this example constant initial conditions are used. Therefore the
+ //argument `globalPos` is not needed. If you want to impose spatial variable initial conditions,
+ // you have to use the `globalPos` argument.
+ // [[codeblock]]
PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
{
- PrimaryVariables values(0.0);
+ PrimaryVariables initialValues(0.0);
// We set the initial water depth to one meter.
- values[0] = 1.0;
+ initialValues[0] = 1.0;
// We set the x-component of the initial velocity to zero.
- values[1] = 0.0;
+ initialValues[1] = 0.0;
// We set the y-component of the initial velocity to zero.
- values[2] = 0.0;
-
- return values;
- };
+ initialValues[2] = 0.0;
- // \}
+ return initialValues;
+ }
+ // [[/codeblock]]
+// We declare the private variables of the problem.
+// [[codeblock]]
private:
- // We declare the private variables of the problem. They are initialized in the problems constructor.
- // We declare the variable for the analytic solution.
+ // variables for the analytic solution.
std::vector<Scalar> exactWaterDepth_;
std::vector<Scalar> exactVelocityX_;
- // constant friction value. An analytic solution is only available for const friction. If you want to run the simulation with a non constant friciton value (specified in the spatialParams) you have to remove the analytic solution.
+ // constant friction value (an analytic solution is only available for const friction).
Scalar constManningN_;
- // The constant bed slope.
+ // The constant channel bed slope.
Scalar bedSlope_;
// The water depth at the outflow boundary.
Scalar hBoundary_;
// The discharge at the inflow boundary.
Scalar discharge_;
- // eps is used as a small value for the definition of the boundry conditions
+ // We assign a private global variable for the epsilon:
static constexpr Scalar eps_ = 1.0e-6;
- std::string name_;
-};
+}; // end class definition RoughChannelProblem
} // end namespace Dumux
-
+// [[/codeblock]]
+// [[/content]]
#endif
diff --git a/examples/shallowwaterfriction/properties.hh b/examples/shallowwaterfriction/properties.hh
index d1e87e9aa4655ca2fa5ce396d4a6f518e40a7d95..86e03e69c14ea38809f30a05fc9e1686b0ab1ee7 100644
--- a/examples/shallowwaterfriction/properties.hh
+++ b/examples/shallowwaterfriction/properties.hh
@@ -20,34 +20,53 @@
#ifndef DUMUX_EXAMPLE_SHALLOWWATER_FRICTION_PROPERTIES_HH
#define DUMUX_EXAMPLE_SHALLOWWATER_FRICTION_PROPERTIES_HH
-// ## The file `properties.hh`
+// ## Compile-time settings (`properties.hh`)
//
-// The header includes will be mentioned in the text below.
-// <details><summary>Click to show the header includes</summary>
+// In this file, the type tag used for the shallow water flow simulation is defined,
+// for which we then specialize `properties` to the needs of the desired setup.
+//
+// [[content]]
+//
+// ### Includes
+// [[details]] include files
+//
+// The `ShallowWater` type tag specializes most of the `properties` required for a
+// shallow water flow simulation in DuMu<sup>x</sup>. We will use this in the following to inherit the
+// respective properties and subsequently specialize those `properties` for our
+// type tag, which we want to modify or for which no meaningful default can be set.
+#include <dumux/freeflow/shallowwater/model.hh>
+
+// We want to use `YaspGrid`, an implementation of the dune grid interface for structured grids:
#include <dune/grid/yaspgrid.hh>
-#include <dumux/common/properties.hh>
+// In this example, we want to discretize the equations with the cell centered finite volume
+// scheme using two-point-flux approximation:
#include <dumux/discretization/cctpfa.hh>
-#include <dumux/freeflow/shallowwater/model.hh>
-#include "spatialparams.hh"
+// We include the problem and spatial parameters headers used for this simulation.
#include "problem.hh"
-// </details>
+#include "spatialparams.hh"
+// [[/details]]
//
-
-// Let's define the properties for our simulation
-namespace Dumux::Properties {
-
-// First, a so-called TypeTag is created. Properties are traits specialized for this TypeTag (a simple `struct`).
-// The properties of two other TypeTags are inherited by adding the alias `InheritsFrom`.
+// ### Type tag definition
+//
+// First, a so-called type tag is created. Properties are traits specialized for this type tag (a simple `struct`).
+// The properties of two other type tags are inherited by adding the alias `InheritsFrom`.
// Here, properties from the shallow water model (`TTag::ShallowWater`) 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
+// are inherited. These other type tag definitions can be found in the included
// headers `dumux/freeflow/shallowwater/model.hh` and `dumux/discretization/cctpfa.hh`.
+// [[codeblock]]
+// We enter the namespace Dumux::Properties in order to import the entire Dumux namespace for general use:
+namespace Dumux::Properties {
+
namespace TTag {
struct RoughChannel { using InheritsFrom = std::tuple<ShallowWater, CCTpfaModel>; };
}
+// [[/codeblock]]
+// ### Property specializations
+//
// We use a structured Cartesian grid with tensor product structure.
// `Dune::YaspGrid` (Yet Another Structure Parallel Grid) is defined in `dune/grid/yaspgrid.hh`
// in the Dune module `dune-grid`.
@@ -55,7 +74,7 @@ template<class TypeTag>
struct Grid<TypeTag, TTag::RoughChannel>
{ using type = Dune::YaspGrid<2, Dune::TensorProductCoordinates<GetPropType<TypeTag, Properties::Scalar>, 2> >; };
-// Next, we specialize the properties `Problem` and `SpatialParams` for our new TypeTag and
+// Next, we specialize the properties `Problem` and `SpatialParams` for our new type tag and
// set the type to our problem and spatial parameter classes implemented
// in `problem.hh` and `spatialparams.hh`.
// [[codeblock]]
@@ -83,5 +102,5 @@ struct EnableGridGeometryCache<TypeTag, TTag::RoughChannel>
{ static constexpr bool value = true; };
} // end namespace Dumux::Properties
-
+// [[/content]]
#endif
diff --git a/examples/shallowwaterfriction/spatialparams.hh b/examples/shallowwaterfriction/spatialparams.hh
index 1638ce5c314c6f16b55a3a88ec138287414fc35d..68f188781db5b2c97e14888a3f4625ec418232c5 100644
--- a/examples/shallowwaterfriction/spatialparams.hh
+++ b/examples/shallowwaterfriction/spatialparams.hh
@@ -20,28 +20,40 @@
#ifndef DUMUX_ROUGH_CHANNEL_SPATIAL_PARAMETERS_HH
#define DUMUX_ROUGH_CHANNEL_SPATIAL_PARAMETERS_HH
-// ## The file `spatialparams.hh`
+// ## Parameter distributions (`spatialparams.hh`)
//
+// This file contains the __spatial parameters class__ which defines the
+// the friction law, including it's friction parameter, the acceleration
+// due to gravity and the altitude of the channel bed surface. In this example only the bed
+// surface has a non constant distribution.
//
-// We include the basic spatial parameters for finite volumes file from which we will inherit
+// [[content]]
+//
+// ### Include files
+// We include the basic spatial parameters file for finite volumes, from which we will inherit.
#include <dumux/material/spatialparams/fv.hh>
-// The parameters header is needed to retrieve run-time parameters.
-#include <dumux/common/parameters.hh>
-// We include all friction laws, between we can choose for the calculation of the bottom friction source.
+// We include all friction laws.
#include <dumux/material/fluidmatrixinteractions/frictionlaws/frictionlaw.hh>
#include <dumux/material/fluidmatrixinteractions/frictionlaws/manning.hh>
#include <dumux/material/fluidmatrixinteractions/frictionlaws/nikuradse.hh>
+#include <dumux/material/fluidmatrixinteractions/frictionlaws/nofriction.hh>
-// We enter the namespace Dumux. All Dumux functions and classes are in a namespace Dumux, to make sure they don`t clash with symbols from other libraries you may want to use in conjunction with Dumux.
+// ### The spatial parameters class
+//
+// In the `RoughChannelSpatialParams` class, we define all functions needed to describe
+// the rough channel for the shallow water problem.
+// We inherit from the `FVSpatialParams` class, which is the base class
+// for spatial parameters in the context of
+// applications using finite volume discretization schemes.
+// [[codeblock]]
namespace Dumux {
-//In the RoughChannelSpatialParams class we define all functions needed to describe the spatial distributed parameters.
template<class GridGeometry, class Scalar, class VolumeVariables>
class RoughChannelSpatialParams
: public FVSpatialParams<GridGeometry, Scalar,
RoughChannelSpatialParams<GridGeometry, Scalar, VolumeVariables>>
{
- // We introduce using declarations that are derived from the property system which we need in this class
+ // This convenience aliases will be used throughout this class
using ThisType = RoughChannelSpatialParams<GridGeometry, Scalar, VolumeVariables>;
using ParentType = FVSpatialParams<GridGeometry, Scalar, ThisType>;
using GridView = typename GridGeometry::GridView;
@@ -49,9 +61,14 @@ class RoughChannelSpatialParams
using SubControlVolume = typename FVElementGeometry::SubControlVolume;
using Element = typename GridView::template Codim<0>::Entity;
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+ // [[/codeblock]]
+ // In the following, the properties of the the rough channel are set. Namely, these are
+ // the friction law, including it's friction parameter, the acceleration
+ // due to gravity and the altitude of the channel bed surface.
+ // [[codeblock]]
public:
- // In the constructor be read some values from the `params.input` and initialize the friciton law.
+ // In the constructor we read some values from the `params.input` and initialize the friciton law.
RoughChannelSpatialParams(std::shared_ptr<const GridGeometry> gridGeometry)
: ParentType(gridGeometry)
{
@@ -61,7 +78,8 @@ public:
initFrictionLaw();
}
- // We initialize the friction law based on the law specified in `params.input`.
+ // This function handles the initialization of the friction law based on the settings
+ // specified in `params.input`.
void initFrictionLaw()
{
if (frictionLawType_ == "Manning")
@@ -74,45 +92,45 @@ public:
Scalar ks = getParam<Scalar>("Problem.Ks");
frictionLaw_ = std::make_unique<FrictionLawNikuradse<VolumeVariables>>(ks);
}
+ else if (frictionLawType_ == "None")
+ {
+ frictionLaw_ = std::make_unique<FrictionLawNoFriction<VolumeVariables>>();
+ }
else
{
- std::cout<<"The FrictionLaw in params.input is unknown. Valid entries are `Manning` and `Nikuradse`!"<<std::endl;
+ std::cout<<"The FrictionLaw in params.input is unknown. Valid entries are `Manning`,"
+ " `Nikuradse` and `None`!"<<std::endl;
}
}
- // Use this function, if you want to vary the value for the gravity.
- Scalar gravity(const GlobalPosition& globalPos) const
+ // This function returns an object of the friction law class, already initialized with a friction value.
+ const FrictionLaw<VolumeVariables>& frictionLaw(const Element& element,
+ const SubControlVolume& scv) const
{
- return gravity_;
+ return *frictionLaw_;
}
- // Use this function for a constant gravity.
- Scalar gravity() const
+ // This function returns the acceleration due to gravity.
+ Scalar gravity(const GlobalPosition& globalPos) const
{
return gravity_;
}
- // This function returns an object of the friction law class, which is initialized with the appropriate friction values. If you want to use different friciton values or laws, you have to use a vector of unique_ptr for `frictionLaw_` and pick the right friction law instances via the `element` argument.
- const FrictionLaw<VolumeVariables>& frictionLaw(const Element& element,
- const SubControlVolume& scv) const
- {
- return *frictionLaw_;
- }
-
- // Define the bed surface based on the `bedSlope_`.
+ // Define the bed surface based on the bed slope and the bed level at the inflow (10 m).
Scalar bedSurface(const Element& element,
const SubControlVolume& scv) const
{
return 10.0 - element.geometry().center()[0] * bedSlope_;
}
+// We declare the private variables of the problem.
private:
Scalar gravity_;
Scalar bedSlope_;
std::string frictionLawType_;
std::unique_ptr<FrictionLaw<VolumeVariables>> frictionLaw_;
-};
-// end of namespace Dumux.
-}
-
+}; // end class definition of RoughChannelSpatialParams
+} // end of namespace Dumux.
+// [[/codeblock]]
+// [[/content]]
#endif