diff --git a/examples/1ptracer/README.md b/examples/1ptracer/README.md
index 864c1f55ee1e2b0624a1319186588544d5c11af9..68e153ff58d6792c1b84c5bb3d9710a591dca2dd 100644
--- a/examples/1ptracer/README.md
+++ b/examples/1ptracer/README.md
@@ -54,21 +54,28 @@ For the tracer model, this is done in the files `problem_tracer.hh` and `spatial
In this file, we generate a random permeability field in the constructor of the `OnePTestSpatialParams` class.
For this, we use the random number generation facilities provided by the C++ standard library.
+
```cpp
#include <random>
```
+
We use the properties for porous medium flow models, declared in the file `properties.hh`.
+
```cpp
#include <dumux/porousmediumflow/properties.hh>
```
+
We include the spatial parameters class for single-phase models discretized by finite volume schemes.
The spatial parameters defined for this example will inherit from those.
+
```cpp
#include <dumux/material/spatialparams/fv1p.hh>
namespace Dumux {
```
+
In the `OnePTestSpatialParams` class, we define all functions needed to describe the porous medium, e.g. porosity and permeability for the 1p_problem.
+
```cpp
template<class GridGeometry, class Scalar>
class OnePTestSpatialParams
@@ -76,7 +83,9 @@ class OnePTestSpatialParams
OnePTestSpatialParams<GridGeometry, Scalar>>
{
```
+
We declare aliases for types that we are going to need in this class.
+
```cpp
using GridView = typename GridGeometry::GridView;
using FVElementGeometry = typename GridGeometry::LocalView;
@@ -90,27 +99,35 @@ We declare aliases for types that we are going to need in this class.
public:
```
+
The spatial parameters must export the type used to define permeabilities.
Here, we are using scalar permeabilities, but tensors are also supported.
+
```cpp
using PermeabilityType = Scalar;
OnePTestSpatialParams(std::shared_ptr<const GridGeometry> gridGeometry)
: ParentType(gridGeometry), K_(gridGeometry->gridView().size(0), 0.0)
{
```
+
### Generation of the random permeability field
We get the permeability of the domain and the lens from the `params.input` file.
+
```cpp
permeability_ = getParam<Scalar>("SpatialParams.Permeability");
permeabilityLens_ = getParam<Scalar>("SpatialParams.PermeabilityLens");
```
+
Furthermore, we get the position of the lens, which is defined by the position of the lower left and the upper right corner.
+
```cpp
lensLowerLeft_ = getParam<GlobalPosition>("SpatialParams.LensLowerLeft");
lensUpperRight_ =getParam<GlobalPosition>("SpatialParams.LensUpperRight");
```
+
We generate random fields for the permeability using lognormal distributions, with `permeability_` as mean value and 10 % of it as standard deviation.
A separate distribution is used for the lens using `permeabilityLens_`. A permeability value is created for each element of the grid and is stored in the vector `K_`.
+
```cpp
std::mt19937 rand(0);
std::lognormal_distribution<Scalar> K(std::log(permeability_), std::log(permeability_)*0.1);
@@ -123,10 +140,12 @@ A separate distribution is used for the lens using `permeabilityLens_`. A permea
}
}
```
+
### Properties of the porous matrix
This function returns the permeability $`[m^2]`$ to be used within a sub-control volume (`scv`) inside the element `element`.
One can define the permeability as function of the primary variables on the element, which are given in the provided `ElementSolution`.
Here, we use element-wise distributed permeabilities that were randomly generated in the constructor (see above).
+
```cpp
template<class ElementSolution>
const PermeabilityType& permeability(const Element& element,
@@ -135,21 +154,26 @@ Here, we use element-wise distributed permeabilities that were randomly generate
{
return K_[scv.dofIndex()];
}
-
```
+
We set the porosity $`[-]`$ for the whole domain to a value of $`20 \%`$.
+
```cpp
Scalar porosityAtPos(const GlobalPosition &globalPos) const
{ return 0.2; }
```
+
We reference to the permeability field. This is used in the main function to write an output for the permeability field.
+
```cpp
const std::vector<Scalar>& getKField() const
{ return K_; }
private:
```
+
We have a convenient definition of the position of the lens.
+
```cpp
bool isInLens_(const GlobalPosition &globalPos) const
{
@@ -171,141 +195,192 @@ We have a convenient definition of the position of the lens.
};
} // end namespace Dumux
-
```
+
## The file `problem_1p.hh`
Before we enter the problem class containing initial and boundary conditions, we include necessary files and introduce properties.
### Include files
We use `YaspGrid`, an implementation of the dune grid interface for structured grids:
+
```cpp
#include <dune/grid/yaspgrid.hh>
```
+
The cell centered, two-point-flux discretization scheme is included:
+
```cpp
#include <dumux/discretization/cctpfa.hh>
```
+
The one-phase flow model is included:
+
```cpp
#include <dumux/porousmediumflow/1p/model.hh>
```
+
This is the porous medium problem class that this class is derived from:
+
```cpp
#include <dumux/porousmediumflow/problem.hh>
```
+
The fluid properties are specified in the following headers (we use liquid water as the fluid phase):
+
```cpp
#include <dumux/material/components/simpleh2o.hh>
#include <dumux/material/fluidsystems/1pliquid.hh>
```
+
The local residual for incompressible flow is included:
+
```cpp
#include <dumux/porousmediumflow/1p/incompressiblelocalresidual.hh>
```
+
We include the header that specifies all spatially variable parameters:
+
```cpp
#include "spatialparams_1p.hh"
```
+
### Define basic properties for our simulation
We enter the namespace Dumux in order to import the entire Dumux namespace for general use
+
```cpp
namespace Dumux {
```
+
The problem class is forward declared:
+
```cpp
template<class TypeTag>
class OnePTestProblem;
```
+
We enter the namespace Properties, which is a sub-namespace of the namespace Dumux:
+
```cpp
namespace Properties {
```
+
A `TypeTag` for our simulation is created which inherits from the one-phase flow model and the cell centered, two-point-flux discretization scheme.
+
```cpp
namespace TTag {
struct IncompressibleTest { using InheritsFrom = std::tuple<OneP, CCTpfaModel>; };
}
```
+
We use a structured 2D grid:
+
```cpp
template<class TypeTag>
struct Grid<TypeTag, TTag::IncompressibleTest> { using type = Dune::YaspGrid<2>; };
```
+
The problem class specifies initial and boundary conditions:
+
```cpp
template<class TypeTag>
struct Problem<TypeTag, TTag::IncompressibleTest> { using type = OnePTestProblem<TypeTag>; };
```
+
We define the spatial parameters for our simulation:
+
```cpp
template<class TypeTag>
struct SpatialParams<TypeTag, TTag::IncompressibleTest>
{
```
+
We define convenient shortcuts to the properties `GridGeometry` and `Scalar`:
+
```cpp
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
```
+
Finally, we set the spatial parameters:
+
```cpp
using type = OnePTestSpatialParams<GridGeometry, Scalar>;
};
```
+
We use the local residual that contains analytic derivative methods for incompressible flow:
+
```cpp
template<class TypeTag>
struct LocalResidual<TypeTag, TTag::IncompressibleTest> { using type = OnePIncompressibleLocalResidual<TypeTag>; };
```
+
In the following we define the fluid properties.
+
```cpp
template<class TypeTag>
struct FluidSystem<TypeTag, TTag::IncompressibleTest>
{
```
+
We define a convenient shortcut to the property Scalar:
+
```cpp
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
```
+
We create a fluid system that consists of one liquid water phase. We use the simple
description of water, which means we do not use tabulated values but more general equations of state.
+
```cpp
using type = FluidSystems::OnePLiquid<Scalar, Components::SimpleH2O<Scalar> >;
};
```
+
We enable caching for the grid volume variables.
+
```cpp
template<class TypeTag>
struct EnableGridVolumeVariablesCache<TypeTag, TTag::IncompressibleTest> { static constexpr bool value = true; };
```
+
We enable caching for the grid flux variables.
+
```cpp
template<class TypeTag>
struct EnableGridFluxVariablesCache<TypeTag, TTag::IncompressibleTest> { static constexpr bool value = true; };
```
+
We enable caching for the FV grid geometry
+
```cpp
template<class TypeTag>
struct EnableGridGeometryCache<TypeTag, TTag::IncompressibleTest> { static constexpr bool value = true; };
```
+
The cache stores values that were already calculated for later usage. This increases the memory demand but makes the simulation faster.
We leave the namespace Properties.
+
```cpp
}
```
+
### The problem class
We enter the problem class where all necessary boundary conditions and initial conditions are set for our simulation.
As this is a porous medium problem, we inherit from the base class `PorousMediumFlowProblem`.
+
```cpp
template<class TypeTag>
class OnePTestProblem : public PorousMediumFlowProblem<TypeTag>
{
```
+
We use convenient declarations that we derive from the property system.
+
```cpp
using ParentType = PorousMediumFlowProblem<TypeTag>;
using GridView = GetPropType<TypeTag, Properties::GridView>;
@@ -321,37 +396,49 @@ We use convenient declarations that we derive from the property system.
public:
```
+
This is the constructor of our problem class:
+
```cpp
OnePTestProblem(std::shared_ptr<const GridGeometry> gridGeometry)
: ParentType(gridGeometry) {}
```
+
First, we define the type of boundary conditions depending on the location. Two types of boundary conditions
can be specified: Dirichlet or Neumann boundary condition. On a Dirichlet boundary, the values of the
primary variables need to be fixed. On a Neumann boundary condition, values for derivatives need to be fixed.
Mixed boundary conditions (different types for different equations on the same boundary) are not accepted for
cell-centered finite volume schemes.
+
```cpp
BoundaryTypes boundaryTypes(const Element &element,
const SubControlVolumeFace &scvf) const
{
BoundaryTypes values;
```
+
we retrieve the global position, i.e. the vector with the global coordinates,
of the integration point on the boundary sub-control volume face `scvf`
+
```cpp
const auto globalPos = scvf.ipGlobal();
```
+
we define a small epsilon value
+
```cpp
Scalar eps = 1.0e-6;
```
+
We specify Dirichlet boundaries on the top and bottom of our domain:
+
```cpp
if (globalPos[dimWorld-1] < eps || globalPos[dimWorld-1] > this->gridGeometry().bBoxMax()[dimWorld-1] - eps)
values.setAllDirichlet();
```
+
The top and bottom of our domain are Neumann boundaries:
+
```cpp
else
values.setAllNeumann();
@@ -359,60 +446,80 @@ The top and bottom of our domain are Neumann boundaries:
return values;
}
```
+
Second, we specify the values for the Dirichlet boundaries. We need to fix values of our primary variable
+
```cpp
PrimaryVariables dirichlet(const Element &element,
const SubControlVolumeFace &scvf) const
{
```
+
we retreive again the global position
+
```cpp
const auto& pos = scvf.ipGlobal();
PrimaryVariables values(0);
```
+
and assign pressure values in [Pa] according to a pressure gradient to 1e5 Pa at the top and 1.1e5 Pa at the bottom.
+
```cpp
values[0] = 1.0e+5*(1.1 - pos[dimWorld-1]*0.1);
return values;
}
```
+
We need to specify a constant temperature for our isothermal problem.
Fluid properties that depend on temperature will be calculated with this value.
+
```cpp
Scalar temperature() const
{
return 283.15; // 10°C
}
```
+
This is everything the one phase problem class contains.
+
```cpp
};
```
+
We leave the namespace Dumux.
+
```cpp
} // end namespace Dumux
```
+
## The file `spatialparams_tracer.hh`
In this file, we define spatial properties of the porous medium such as the permeability and the porosity in various functions for the tracer problem.
Furthermore, spatial dependent properties of the tracer fluid system are defined and in the end two functions handle the calculated volume fluxes from the solution of the 1p problem.
We use the properties for porous medium flow models, declared in the file `properties.hh`.
+
```cpp
#include <dumux/porousmediumflow/properties.hh>
```
+
As in the 1p spatialparams, we inherit from the spatial parameters for single-phase models using finite volumes, which we include here.
+
```cpp
#include <dumux/material/spatialparams/fv1p.hh>
```
+
We enter the namespace Dumux
+
```cpp
namespace Dumux {
```
+
In the `TracerTestSpatialParams` class, we define all functions needed to describe spatially dependent parameters for the `tracer_problem`.
+
```cpp
template<class GridGeometry, class Scalar>
class TracerTestSpatialParams
@@ -435,13 +542,17 @@ public:
TracerTestSpatialParams(std::shared_ptr<const GridGeometry> gridGeometry)
: ParentType(gridGeometry) {}
```
+
### Properties of the porous matrix
We define the same porosity for the whole domain as in the 1p spatialparams.
+
```cpp
Scalar porosityAtPos(const GlobalPosition& globalPos) const
{ return 0.2; }
```
+
We do not consider dispersivity for the tracer transport. So we set the dispersivity coefficient to zero.
+
```cpp
template<class ElementSolution>
Scalar dispersivity(const Element &element,
@@ -449,31 +560,39 @@ We do not consider dispersivity for the tracer transport. So we set the dispersi
const ElementSolution& elemSol) const
{ return 0; }
```
+
### Properties of the fluid system
In the following, we define fluid properties that are spatial parameters in the tracer model.
They can possible vary in space but are usually constants.
Furthermore, spatially constant values of the fluid system are defined in the `TracerFluidSystem` class in `problem.hh`.
We define the fluid density to a constant value of 1000 $`\frac{kg}{m^3}`$.
+
```cpp
Scalar fluidDensity(const Element &element,
const SubControlVolume& scv) const
{ return 1000; }
```
+
This interface defines the fluid molar mass within the sub-control volume `scv`.
+
```cpp
Scalar fluidMolarMass(const Element &element,
const SubControlVolume& scv) const
{ return fluidMolarMassAtPos(scv.dofPosition()); }
```
+
This interface defines the fluid molar mass depending on the position in the domain.
+
```cpp
Scalar fluidMolarMassAtPos(const GlobalPosition &globalPos) const
{ return 18.0; }
```
+
### The volume fluxes
We define a function which returns the volume flux across the given sub-control volume face `scvf`.
This flux is obtained from the vector `volumeFlux_` that contains the fluxes across al sub-control volume faces of the discretization.
This vector can be set using the `setVolumeFlux` function.
+
```cpp
template<class ElementVolumeVariables>
Scalar volumeFlux(const Element &element,
@@ -484,8 +603,10 @@ This vector can be set using the `setVolumeFlux` function.
return volumeFlux_[scvf.index()];
}
```
+
We define a function that allows setting the volume fluxes for all sub-control volume faces of the discretization.
This is used in the main function after these fluxes have been based on the pressure solution obtained with the single-phase model.
+
```cpp
void setVolumeFlux(const std::vector<Scalar>& f)
{ volumeFlux_ = f; }
@@ -495,63 +616,84 @@ private:
};
} // end namespace Dumux
-
```
+
## The file `problem_tracer.hh`
Before we enter the problem class containing initial and boundary conditions, we include necessary files and introduce properties.
### Include files
Again, we use YaspGrid, the implementation of the dune grid interface for structured grids:
+
```cpp
#include <dune/grid/yaspgrid.hh>
```
+
and the cell centered, two-point-flux discretization.
+
```cpp
#include <dumux/discretization/cctpfa.hh>
```
+
Here, we include the tracer model:
+
```cpp
#include <dumux/porousmediumflow/tracer/model.hh>
```
+
We include again the porous medium problem class that this class is derived from:
+
```cpp
#include <dumux/porousmediumflow/problem.hh>
```
+
and the base fluidsystem. We will define a custom fluid system that inherits from that class.
+
```cpp
#include <dumux/material/fluidsystems/base.hh>
```
+
We include the header that specifies all spatially variable parameters for the tracer problem:
+
```cpp
#include "spatialparams_tracer.hh"
```
+
### Define basic properties for our simulation
We enter the namespace Dumux
+
```cpp
namespace Dumux {
```
+
The problem class is forward declared:
+
```cpp
template <class TypeTag>
class TracerTestProblem;
```
+
We enter the namespace Properties,
+
```cpp
namespace Properties {
```
+
A `TypeTag` for our simulation is created which inherits from the tracer model and the
cell centered, two-point-flux discretization scheme.
+
```cpp
namespace TTag {
struct TracerTest { using InheritsFrom = std::tuple<Tracer>; };
struct TracerTestCC { using InheritsFrom = std::tuple<TracerTest, CCTpfaModel>; };
}
```
+
We enable caching for the grid volume variables, the flux variables and the FV grid geometry.
+
```cpp
template<class TypeTag>
struct EnableGridVolumeVariablesCache<TypeTag, TTag::TracerTest> { static constexpr bool value = true; };
@@ -560,17 +702,23 @@ struct EnableGridFluxVariablesCache<TypeTag, TTag::TracerTest> { static constexp
template<class TypeTag>
struct EnableGridGeometryCache<TypeTag, TTag::TracerTest> { static constexpr bool value = true; };
```
+
We use the same grid as in the stationary one-phase model, a structured 2D grid:
+
```cpp
template<class TypeTag>
struct Grid<TypeTag, TTag::TracerTest> { using type = Dune::YaspGrid<2>; };
```
+
The problem class specifies initial and boundary conditions:
+
```cpp
template<class TypeTag>
struct Problem<TypeTag, TTag::TracerTest> { using type = TracerTestProblem<TypeTag>; };
```
+
We define the spatial parameters for our tracer simulation:
+
```cpp
template<class TypeTag>
struct SpatialParams<TypeTag, TTag::TracerTest>
@@ -582,27 +730,35 @@ public:
using type = TracerTestSpatialParams<GridGeometry, Scalar>;
};
```
+
One can choose between a formulation in terms of mass or mole fractions. Here, we are using mass fractions.
+
```cpp
template<class TypeTag>
struct UseMoles<TypeTag, TTag::TracerTest> { static constexpr bool value = false; };
```
+
We use solution-independent molecular diffusion coefficients. Per default, solution-dependent
diffusion coefficients are assumed during the computation of the jacobian matrix entries. Specifying
solution-independent diffusion coefficients can speed up computations:
+
```cpp
template<class TypeTag>
struct SolutionDependentMolecularDiffusion<TypeTag, TTag::TracerTestCC>
{ static constexpr bool value = false; };
```
+
In the following, we create a new tracer fluid system and derive from the base fluid system.
+
```cpp
template<class TypeTag>
class TracerFluidSystem : public FluidSystems::Base<GetPropType<TypeTag, Properties::Scalar>,
TracerFluidSystem<TypeTag>>
{
```
+
We define some convenience aliases:
+
```cpp
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using Problem = GetPropType<TypeTag, Properties::Problem>;
@@ -613,41 +769,48 @@ We define some convenience aliases:
public:
```
+
We specify that the fluid system only contains tracer components,
+
```cpp
static constexpr bool isTracerFluidSystem()
{ return true; }
```
-and that no component is the main component
-```cpp
- static constexpr int getMainComponent(int phaseIdx)
- { return -1; }
-```
+
We define the number of components of this fluid system (one single tracer component)
+
```cpp
static constexpr int numComponents = 1;
```
+
This interface is designed to define the names of the components of the fluid system.
Here, we only have a single component, so `compIdx` should always be 0.
The component name is used for the vtk output.
+
```cpp
static std::string componentName(int compIdx = 0)
{ return "tracer_" + std::to_string(compIdx); }
```
+
We set the phase name for the phase index (`phaseIdx`) for velocity vtk output:
Here, we only have a single phase, so `phaseIdx` should always be zero.
+
```cpp
static std::string phaseName(int phaseIdx = 0)
{ return "Groundwater"; }
```
+
We set the molar mass of the tracer component with index `compIdx` (should again always be zero here).
+
```cpp
static Scalar molarMass(unsigned int compIdx = 0)
{ return 0.300; }
```
+
We set the value for the binary diffusion coefficient. This
might depend on spatial parameters like pressure / temperature.
But, in this case we neglect diffusion and return 0.0:
+
```cpp
static Scalar binaryDiffusionCoefficient(unsigned int compIdx,
const Problem& problem,
@@ -656,25 +819,32 @@ But, in this case we neglect diffusion and return 0.0:
{ return 0.0; }
};
```
+
We set the above created tracer fluid system:
+
```cpp
template<class TypeTag>
struct FluidSystem<TypeTag, TTag::TracerTest> { using type = TracerFluidSystem<TypeTag>; };
```
+
We leave the namespace Properties.
+
```cpp
}
-
```
+
### The problem class
We enter the problem class where all necessary boundary conditions and initial conditions are set for our simulation.
As this is a porous medium problem, we inherit from the base class `PorousMediumFlowProblem`.
+
```cpp
template <class TypeTag>
class TracerTestProblem : public PorousMediumFlowProblem<TypeTag>
{
```
+
We use convenient declarations that we derive from the property system.
+
```cpp
using ParentType = PorousMediumFlowProblem<TypeTag>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
@@ -694,24 +864,32 @@ We use convenient declarations that we derive from the property system.
using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
```
+
We create a convenience bool stating whether mole or mass fractions are used
+
```cpp
static constexpr bool useMoles = getPropValue<TypeTag, Properties::UseMoles>();
```
+
We create additional convenience integers to make dimWorld and numComponents available in the problem
+
```cpp
static constexpr int dimWorld = GridView::dimensionworld;
static const int numComponents = FluidSystem::numComponents;
public:
```
+
This is the constructor of our problem class:
+
```cpp
TracerTestProblem(std::shared_ptr<const GridGeometry> gridGeometry)
: ParentType(gridGeometry)
{
```
+
We print to the terminal whether mole or mass fractions are used
+
```cpp
if(useMoles)
std::cout<<"problem uses mole fractions" << '\n';
@@ -719,8 +897,10 @@ We print to the terminal whether mole or mass fractions are used
std::cout<<"problem uses mass fractions" << '\n';
}
```
+
We define the type of boundary conditions depending on the location.
All boundaries are set to a neumann-type flow boundary condition.
+
```cpp
BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
{
@@ -729,19 +909,25 @@ All boundaries are set to a neumann-type flow boundary condition.
return values;
}
```
+
We specify the initial conditions for the primary variable (tracer concentration) depending on the location.
+
```cpp
PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
{
PrimaryVariables initialValues(0.0);
```
+
The initial contamination is located at the bottom of the domain:
+
```cpp
if (globalPos[1] < 0.1 + eps_)
{
```
+
We chose a mole fraction of $`1e-9`$, but in case the mass fractions
are used by the model, we have to convert this value:
+
```cpp
if (useMoles)
initialValues = 1e-9;
@@ -752,7 +938,9 @@ are used by the model, we have to convert this value:
return initialValues;
}
```
+
We implement an outflow boundary on the top of the domain and prescribe zero-flux Neumann boundary conditions on all other boundaries.
+
```cpp
NumEqVector neumann(const Element& element,
const FVElementGeometry& fvGeometry,
@@ -764,7 +952,9 @@ We implement an outflow boundary on the top of the domain and prescribe zero-flu
const auto& volVars = elemVolVars[scvf.insideScvIdx()];
const auto& globalPos = scvf.center();
```
+
This is the outflow boundary, where tracer is transported by advection with the given flux field.
+
```cpp
if (globalPos[dimWorld-1] > this->gridGeometry().bBoxMax()[dimWorld-1] - eps_)
{
@@ -774,7 +964,9 @@ This is the outflow boundary, where tracer is transported by advection with the
assert(values>=0.0 && "Volume flux at outflow boundary is expected to have a positive sign");
}
```
+
Prescribe zero-flux Neumann boundary conditions elsewhere
+
```cpp
else
values = 0.0;
@@ -784,140 +976,194 @@ Prescribe zero-flux Neumann boundary conditions elsewhere
private:
```
+
We assign a private global variable for the epsilon:
+
```cpp
static constexpr Scalar eps_ = 1e-6;
```
+
This is everything the tracer problem class contains.
+
```cpp
};
```
+
We leave the namespace Dumux here.
+
```cpp
} // end namespace Dumux
-
```
-## The file `main.cc`
+
+## Main program flow (`main.cc`)
-We look now at the main file for the tracer problem. We set up two problems in this file and solve them sequentially, first the 1p problem and afterwards the tracer problem. The result of the 1p problem is the pressure distribution in the problem domain. We use it to calculate the volume fluxes, which act as an input for the tracer problem. Based on this volume fluxes, we calculate the transport of a tracer in the following tracer problem.
-### Includes
+This file contains the main program flow. In this example, we solve a single-phase flow problem
+to obtain a pressure distribution on the domain. Subsequently, the distribution of volume fluxes
+is computed from that pressure distribution, which is then passed to a tracer problem to solve
+the transport of an initial contamination through the model domain.
+### Included header files
+
```cpp
#include <config.h>
```
+
We include both problems in the main file, the `problem_1p.hh` and the `problem_tracer.hh`.
+
```cpp
#include "problem_1p.hh"
#include "problem_tracer.hh"
```
+
Further, we include a standard header file for C++, to get time and date information
+
```cpp
#include <ctime>
```
+
and another one 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.
Here, we include 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>
```
+
In Dumux, the property system is used to specify classes and compile-time options to be used by 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 and retrieval 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>
```
+
The following file contains the class, which defines the sequential linear solver backends.
+
```cpp
#include <dumux/linear/seqsolverbackend.hh>
```
+
Further we include the assembler, which assembles the linear systems for finite volume schemes (box-scheme, tpfa-approximation, mpfa-approximation).
+
```cpp
#include <dumux/assembly/fvassembler.hh>
```
+
The containing class in the following file defines the different differentiation methods used to compute the derivatives of the residual.
+
```cpp
#include <dumux/assembly/diffmethod.hh>
```
+
We need the following class to simplify the writing of dumux simulation data to VTK format.
+
```cpp
#include <dumux/io/vtkoutputmodule.hh>
```
+
The gridmanager constructs a grid from the information in the input or grid file. There is a specification for the different supported grid managers.
+
```cpp
#include <dumux/io/grid/gridmanager.hh>
```
+
### Beginning of the main function
+
```cpp
int main(int argc, char** argv) try
{
using namespace Dumux;
```
+
Convenience aliases for the type tags of the two problems, which are defined in the individual problem files.
+
```cpp
using OnePTypeTag = Properties::TTag::IncompressibleTest;
using TracerTypeTag = Properties::TTag::TracerTestCC;
```
+
We initialize MPI. Finalization is done automatically on exit.
+
```cpp
const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
```
+
We print the dumux start message.
+
```cpp
if (mpiHelper.rank() == 0)
DumuxMessage::print(/*firstCall=*/true);
```
+
We parse the command line arguments.
+
```cpp
Parameters::init(argc, argv);
```
+
### 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, by specifying the coordinates
of the corners of the grid and the number of cells to be used to discretize each spatial direction.
Here, we solve both the single-phase and the tracer problem on the same grid.
Hence, the grid is only created once using the grid type defined by the type tag of the 1p problem.
+
```cpp
GridManager<GetPropType<OnePTypeTag, Properties::Grid>> gridManager;
gridManager.init();
```
+
We compute on the leaf grid view.
+
```cpp
const auto& leafGridView = gridManager.grid().leafGridView();
```
+
### Set-up and solving of the 1p problem
In the following section, we set up and solve the 1p problem. As the result of this problem, we obtain the pressure distribution in the domain.
#### Set-up
We create and initialize the finite volume grid geometry, the problem, the linear system, including the jacobian matrix, the residual and the solution vector and the gridvariables.
We need the finite volume geometry to build up the subcontrolvolumes (scv) and subcontrolvolume faces (scvf) for each element of the grid partition.
+
```cpp
using GridGeometry = GetPropType<OnePTypeTag, Properties::GridGeometry>;
auto gridGeometry = std::make_shared<GridGeometry>(leafGridView);
gridGeometry->update();
```
+
In the problem, we define the boundary and initial conditions.
+
```cpp
using OnePProblem = GetPropType<OnePTypeTag, Properties::Problem>;
auto problemOneP = std::make_shared<OnePProblem>(gridGeometry);
```
+
The jacobian matrix (`A`), the solution vector (`p`) and the residual (`r`) are parts of the linear system.
+
```cpp
using JacobianMatrix = GetPropType<OnePTypeTag, Properties::JacobianMatrix>;
using SolutionVector = GetPropType<OnePTypeTag, Properties::SolutionVector>;
@@ -926,32 +1172,42 @@ The jacobian matrix (`A`), the solution vector (`p`) and the residual (`r`) are
auto A = std::make_shared<JacobianMatrix>();
auto r = std::make_shared<SolutionVector>();
```
+
The grid variables store variables (primary and secondary variables) on sub-control volumes and faces (volume and flux variables).
+
```cpp
using OnePGridVariables = GetPropType<OnePTypeTag, Properties::GridVariables>;
auto onePGridVariables = std::make_shared<OnePGridVariables>(problemOneP, gridGeometry);
onePGridVariables->init(p);
```
+
#### Assembling the linear system
We create and inizialize the assembler.
+
```cpp
using OnePAssembler = FVAssembler<OnePTypeTag, DiffMethod::analytic>;
auto assemblerOneP = std::make_shared<OnePAssembler>(problemOneP, gridGeometry, onePGridVariables);
assemblerOneP->setLinearSystem(A, r);
```
+
We assemble the local jacobian and the residual and stop the time needed, which is displayed in the terminal output, using the `assemblyTimer`. Further, we start the timer to evaluate the total time of the assembly, solving and updating.
+
```cpp
Dune::Timer timer;
Dune::Timer assemblyTimer; std::cout << "Assembling linear system ..." << std::flush;
assemblerOneP->assembleJacobianAndResidual(p);
assemblyTimer.stop(); std::cout << " took " << assemblyTimer.elapsed() << " seconds." << std::endl;
```
+
We want to solve `Ax = -r`.
+
```cpp
(*r) *= -1.0;
```
+
#### Solution
We set the linear solver "UMFPack" as the linear solver. Afterwards we solve the linear system. The time needed to solve the system is recorded by the `solverTimer` and displayed in the terminal output.
+
```cpp
using LinearSolver = UMFPackBackend;
Dune::Timer solverTimer; std::cout << "Solving linear system ..." << std::flush;
@@ -959,15 +1215,18 @@ We set the linear solver "UMFPack" as the linear solver. Afterwards we solve the
linearSolver->solve(*A, p, *r);
solverTimer.stop(); std::cout << " took " << solverTimer.elapsed() << " seconds." << std::endl;
```
+
#### Update and output
We update the grid variables with the new solution.
+
```cpp
Dune::Timer updateTimer; std::cout << "Updating variables ..." << std::flush;
onePGridVariables->update(p);
updateTimer.elapsed(); std::cout << " took " << updateTimer.elapsed() << std::endl;
-
```
+
We initialize the vtkoutput. Each model has a predefined model specific output with relevant parameters for that model. We add the pressure data from the solution vector (`p`) and the permeability field as output data.
+
```cpp
using GridView = GetPropType<OnePTypeTag, Properties::GridView>;
Dune::VTKWriter<GridView> onepWriter(leafGridView);
@@ -976,7 +1235,9 @@ We initialize the vtkoutput. Each model has a predefined model specific output w
onepWriter.addCellData(k, "permeability");
onepWriter.write("1p");
```
+
We stop the timer and display the total time of the simulation as well as the cumulative CPU time.
+
```cpp
timer.stop();
@@ -984,10 +1245,11 @@ We stop the timer and display the total time of the simulation as well as the cu
std::cout << "Simulation took " << timer.elapsed() << " seconds on "
<< comm.size() << " processes.\n"
<< "The cumulative CPU time was " << timer.elapsed()*comm.size() << " seconds.\n";
-
```
+
### Computation of the volume fluxes
We use the results of the 1p problem to calculate the volume fluxes in the model domain.
+
```cpp
using Scalar = GetPropType<OnePTypeTag, Properties::Scalar>;
std::vector<Scalar> volumeFlux(gridGeometry->numScvf(), 0.0);
@@ -995,13 +1257,17 @@ We use the results of the 1p problem to calculate the volume fluxes in the model
using FluxVariables = GetPropType<OnePTypeTag, Properties::FluxVariables>;
auto upwindTerm = [](const auto& volVars) { return volVars.mobility(0); };
```
+
We iterate over all elements
+
```cpp
for (const auto& element : elements(leafGridView))
{
```
+
Compute the element-local views on geometry, primary and secondary variables
as well as variables needed for flux computations
+
```cpp
auto fvGeometry = localView(*gridGeometry);
fvGeometry.bind(element);
@@ -1012,66 +1278,87 @@ as well as variables needed for flux computations
auto elemFluxVars = localView(onePGridVariables->gridFluxVarsCache());
elemFluxVars.bind(element, fvGeometry, elemVolVars);
```
+
We calculate the volume fluxes for all sub-control volume faces except for Neumann boundary faces
+
```cpp
for (const auto& scvf : scvfs(fvGeometry))
{
```
+
skip Neumann boundary faces
+
```cpp
if (scvf.boundary() && problemOneP->boundaryTypes(element, scvf).hasNeumann())
continue;
```
+
let the `FluxVariables` class do the flux computation.
+
```cpp
FluxVariables fluxVars;
fluxVars.init(*problemOneP, element, fvGeometry, elemVolVars, scvf, elemFluxVars);
volumeFlux[scvf.index()] = fluxVars.advectiveFlux(0, upwindTerm);
}
}
-
```
+
### Set-up and solving of the tracer problem
#### Set-up
Similar to the 1p problem, we first create and initialize the problem.
+
```cpp
using TracerProblem = GetPropType<TracerTypeTag, Properties::Problem>;
auto tracerProblem = std::make_shared<TracerProblem>(gridGeometry);
```
+
We use the volume fluxes calculated in the previous section as input for the tracer model.
+
```cpp
tracerProblem->spatialParams().setVolumeFlux(volumeFlux);
```
+
We create and initialize the solution vector. As the tracer problem is transient, the initial solution defined in the problem is applied to the solution vector.
+
```cpp
SolutionVector x(leafGridView.size(0));
tracerProblem->applyInitialSolution(x);
auto xOld = x;
```
+
We create and initialize the grid variables.
+
```cpp
using GridVariables = GetPropType<TracerTypeTag, Properties::GridVariables>;
auto gridVariables = std::make_shared<GridVariables>(tracerProblem, gridGeometry);
gridVariables->init(x);
```
+
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.
+
```cpp
const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
auto dt = getParam<Scalar>("TimeLoop.DtInitial");
const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
```
+
We instantiate the time loop.
+
```cpp
auto timeLoop = std::make_shared<CheckPointTimeLoop<Scalar>>(0.0, dt, tEnd);
timeLoop->setMaxTimeStepSize(maxDt);
```
+
We create and inizialize the assembler with time loop for the instationary problem.
+
```cpp
using TracerAssembler = FVAssembler<TracerTypeTag, DiffMethod::analytic, /*implicit=*/false>;
auto assembler = std::make_shared<TracerAssembler>(tracerProblem, gridGeometry, gridVariables, timeLoop);
assembler->setLinearSystem(A, r);
```
+
We initialize the vtk output module and add a velocity output.
+
```cpp
VtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables, x, tracerProblem->name());
using IOFields = GetPropType<TracerTypeTag, Properties::IOFields>;
@@ -1079,43 +1366,56 @@ We initialize the vtk output module and add a velocity output.
using VelocityOutput = GetPropType<TracerTypeTag, Properties::VelocityOutput>;
vtkWriter.addVelocityOutput(std::make_shared<VelocityOutput>(*gridVariables));
vtkWriter.write(0.0);
-
```
+
We define 10 check points in the time loop at which we will write the solution to vtk files.
+
```cpp
timeLoop->setPeriodicCheckPoint(tEnd/10.0);
```
+
#### The time loop
We start the time loop and calculate a new time step as long as `tEnd` is not reached. In every single time step, the problem is assembled and solved.
+
```cpp
timeLoop->start(); do
{
```
+
First we define the old solution as the solution of the previous time step for storage evaluations.
+
```cpp
assembler->setPreviousSolution(xOld);
```
+
Then the linear system is assembled.
+
```cpp
Dune::Timer assembleTimer;
assembler->assembleJacobianAndResidual(x);
assembleTimer.stop();
```
+
We solve the linear system `A(xOld-xNew) = r`.
+
```cpp
Dune::Timer solveTimer;
SolutionVector xDelta(x);
linearSolver->solve(*A, xDelta, *r);
solveTimer.stop();
```
+
We calculate the actual solution and update it in the grid variables.
+
```cpp
updateTimer.reset();
x -= xDelta;
gridVariables->update(x);
updateTimer.stop();
```
+
We display the statistics of the actual time step.
+
```cpp
const auto elapsedTot = assembleTimer.elapsed() + solveTimer.elapsed() + updateTimer.elapsed();
std::cout << "Assemble/solve/update time: "
@@ -1124,34 +1424,45 @@ We display the statistics of the actual time step.
<< updateTimer.elapsed() << "(" << 100*updateTimer.elapsed()/elapsedTot << "%)"
<< std::endl;
```
+
The new solution is defined as the old solution.
+
```cpp
xOld = x;
gridVariables->advanceTimeStep();
```
+
We advance the time loop to the next time step.
+
```cpp
timeLoop->advanceTimeStep();
```
+
We write the Vtk output on check points.
+
```cpp
if (timeLoop->isCheckPoint())
vtkWriter.write(timeLoop->time());
```
+
We report the statistics of this time step.
+
```cpp
timeLoop->reportTimeStep();
```
+
We set the time step size dt of the next time step.
+
```cpp
timeLoop->setTimeStepSize(dt);
} while (!timeLoop->finished());
timeLoop->finalize(leafGridView.comm());
-
```
+
### Final Output
+
```cpp
if (mpiHelper.rank() == 0)
DumuxMessage::print(/*firstCall=*/false);
@@ -1160,9 +1471,11 @@ We set the time step size dt of the next time step.
}
```
+
### Exception handling
In this part of the main file we catch and print possible exceptions that could
occur during the simulation.
+
```cpp
catch (Dumux::ParameterException &e)
{
@@ -1190,6 +1503,7 @@ catch (...)
}
```
+
## Results
The 1p-model calculated a stationary pressure distribution. It is shown in the following figure:
diff --git a/examples/2pinfiltration/README.md b/examples/2pinfiltration/README.md
index d068fbcda1c377e242140a6241bb829766da3188..32fe01bdd88dedca5bfe1769ff5d197d7029a25f 100644
--- a/examples/2pinfiltration/README.md
+++ b/examples/2pinfiltration/README.md
@@ -50,17 +50,22 @@ The grid is adapitvely refined around the injection. The adaptive behaviour can
we include the basic spatial parameters for finite volumes file from which we will inherit
+
```cpp
#include <dumux/material/spatialparams/fv.hh>
```
+
we include all laws which are needed to define the interaction between the solid matrix and the fluids, e.g. laws for capillary pressure saturation relationships.
+
```cpp
#include <dumux/material/fluidmatrixinteractions/2p/regularizedvangenuchten.hh>
#include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh>
namespace Dumux {
```
+
In the TwoPTestSpatialParams class we define all functions needed to describe the porous matrix, e.g. porosity and permeability
+
```cpp
template<class GridGeometry, class Scalar>
@@ -68,7 +73,9 @@ class TwoPTestSpatialParams
: public FVSpatialParams<GridGeometry, Scalar, TwoPTestSpatialParams<GridGeometry, Scalar>>
{
```
+
we introduce using declarations that are derived from the property system which we need in this class
+
```cpp
using GridView = typename GridGeometry::GridView;
using Element = typename GridView::template Codim<0>::Entity;
@@ -91,32 +98,42 @@ public:
: ParentType(gridGeometry)
{
```
+
we get the position of the lens from the params.input file. The lens is defined by the position of the lower left and the upper right corner
+
```cpp
lensLowerLeft_ = getParam<GlobalPosition>("SpatialParams.LensLowerLeft");
lensUpperRight_ = getParam<GlobalPosition>("SpatialParams.LensUpperRight");
```
+
we set the parameters for the material law (here Van-Genuchten Law). First we set the residual saturations for the wetting phase and the non-wetting phase. lensMaterialParams_ define the material parameters for the lens while outerMaterialParams_ define material params for the rest of the domain.
+
```cpp
lensMaterialParams_.setSwr(0.18);
lensMaterialParams_.setSnr(0.0);
outerMaterialParams_.setSwr(0.05);
outerMaterialParams_.setSnr(0.0);
```
+
we set the parameters for the Van Genuchten law alpha and n
+
```cpp
lensMaterialParams_.setVgAlpha(0.00045);
lensMaterialParams_.setVgn(7.3);
outerMaterialParams_.setVgAlpha(0.0037);
outerMaterialParams_.setVgn(4.7);
```
+
here we get the permeabilities from the params.input file. In case that no parameter is set, the default parameters (9.05e-12 and 4.6e-10) are used
+
```cpp
lensK_ = getParam<Scalar>("SpatialParams.lensK", 9.05e-12);
outerK_ = getParam<Scalar>("SpatialParams.outerK", 4.6e-10);
}
```
+
We define the (intrinsic) permeability $`[m^2]`$. In this test, we use element-wise distributed permeabilities.
+
```cpp
template<class ElementSolution>
PermeabilityType permeability(const Element& element,
@@ -128,7 +145,9 @@ We define the (intrinsic) permeability $`[m^2]`$. In this test, we use element-w
return outerK_;
}
```
+
We set the porosity $`[-]`$ depending on the position
+
```cpp
Scalar porosityAtPos(const GlobalPosition& globalPos) const
{
@@ -137,7 +156,9 @@ We set the porosity $`[-]`$ depending on the position
return 0.4;
}
```
+
We set the parameter object for the Van Genuchten material law.
+
```cpp
template<class ElementSolution>
const MaterialLawParams& materialLawParams(const Element& element,
@@ -148,9 +169,10 @@ We set the parameter object for the Van Genuchten material law.
return lensMaterialParams_;
return outerMaterialParams_;
}
-
```
+
Here we can define which phase is to be considered as the wetting phase. Here the wetting phase is the the first phase of the fluidsystem. In this case that is water.
+
```cpp
template<class FluidSystem>
int wettingPhaseAtPos(const GlobalPosition& globalPos) const
@@ -158,7 +180,9 @@ Here we can define which phase is to be considered as the wetting phase. Here th
private:
```
+
we have a convenience definition of the position of the lens
+
```cpp
bool isInLens_(const GlobalPosition &globalPos) const
{
@@ -181,140 +205,189 @@ we have a convenience definition of the position of the lens
};
} // end namespace Dumux
-
```
+
## The file `problem.hh`
-## Include files
+### Include files
The grid we use
+
```cpp
#include <dune/alugrid/grid.hh>
```
+
The cell centered, two-point-flux discretization scheme is included:
+
```cpp
#include <dumux/discretization/cctpfa.hh>
```
+
The fluid properties are specified in the following headers:
+
```cpp
#include <dumux/material/components/trichloroethene.hh>
#include <dumux/material/components/simpleh2o.hh>
#include <dumux/material/fluidsystems/1pliquid.hh>
#include <dumux/material/fluidsystems/2pimmiscible.hh>
```
+
This is the porous medium problem class that this class is derived from:
+
```cpp
#include <dumux/porousmediumflow/problem.hh>
```
+
The two-phase flow model is included:
+
```cpp
#include <dumux/porousmediumflow/2p/model.hh>
```
+
The local residual for incompressible flow is included:
+
```cpp
#include <dumux/porousmediumflow/2p/incompressiblelocalresidual.hh>
```
+
We include the header that specifies all spatially variable parameters:
+
```cpp
#include "spatialparams.hh"
```
+
A container to read values for the initial condition is included:
+
```cpp
#include <dumux/io/container.hh>
```
+
## Define basic properties for our simulation
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.
One could use these functions and classes by prefixing every use of these names by ::,
but that would quickly become cumbersome and annoying.
Rather, we simply import the entire Dumux namespace for general use:
+
```cpp
namespace Dumux {
```
+
The problem class is forward declared:
+
```cpp
template<class TypeTag>
class PointSourceProblem;
```
+
We enter the namespace Properties, which is a sub-namespace of the namespace Dumux:
+
```cpp
namespace Properties {
```
+
A TypeTag for our simulation is created which inherits from the two-phase flow model and the
cell centered, two-point-flux discretization scheme.
+
```cpp
namespace TTag {
struct PointSourceExample { using InheritsFrom = std::tuple<TwoP, CCTpfaModel>; };
}
```
+
We use non-conforming refinement in our simulation:
+
```cpp
template<class TypeTag>
struct Grid<TypeTag, TTag::PointSourceExample> { using type = Dune::ALUGrid<2, 2, Dune::cube, Dune::nonconforming>; };
```
+
The problem class specifies initial and boundary conditions:
+
```cpp
template<class TypeTag>
struct Problem<TypeTag, TTag::PointSourceExample> { using type = PointSourceProblem<TypeTag>; };
```
+
The local residual contains analytic derivative methods for incompressible flow:
+
```cpp
template<class TypeTag>
struct LocalResidual<TypeTag, TTag::PointSourceExample> { using type = TwoPIncompressibleLocalResidual<TypeTag>; };
```
+
In the following we define our fluid properties.
+
```cpp
template<class TypeTag>
struct FluidSystem<TypeTag, TTag::PointSourceExample>
{
```
+
We define a convenient shortcut to the property Scalar:
+
```cpp
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
```
+
First, we create a fluid system that consists of one liquid water phase. We use the simple
description of water, which means we do not use tabulated values but more general equations of state.
+
```cpp
using WettingPhase = FluidSystems::OnePLiquid<Scalar, Components::SimpleH2O<Scalar> >;
```
+
Second, we create another fluid system consisting of a liquid phase as well, the Trichlorethene (DNAPL) phase:
+
```cpp
using NonwettingPhase = FluidSystems::OnePLiquid<Scalar, Components::Trichloroethene<Scalar> >;
```
+
Third, we combine both fluid systems in our final fluid system which consist of two immiscible liquid phases:
+
```cpp
using type = FluidSystems::TwoPImmiscible<Scalar, WettingPhase, NonwettingPhase>;
};
```
+
we set the formulation for the primary variables to p0s1. In this case that means that the water pressure and the DNAPL saturation are our primary variables.
+
```cpp
template<class TypeTag>
struct Formulation<TypeTag, TTag::PointSourceExample>
{ static constexpr auto value = TwoPFormulation::p0s1; };
```
+
We define the spatial parameters for our simulation:
+
```cpp
template<class TypeTag>
struct SpatialParams<TypeTag, TTag::PointSourceExample>
{
```
+
We define convenient shortcuts to the properties GridGeometry and Scalar:
+
```cpp
private:
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
```
+
Finally we set the spatial parameters:
+
```cpp
public:
using type = TwoPTestSpatialParams<GridGeometry, Scalar>;
};
```
+
We enable caching for the grid volume variables, the grid flux variables and the FV grid geometry. The cache
stores values that were already calculated for later usage. This makes the simulation faster.
+
```cpp
template<class TypeTag>
struct EnableGridVolumeVariablesCache<TypeTag, TTag::PointSourceExample> { static constexpr bool value = false; };
@@ -323,19 +396,25 @@ stores values that were already calculated for later usage. This makes the simul
template<class TypeTag>
struct EnableGridGeometryCache<TypeTag, TTag::PointSourceExample> { static constexpr bool value = false; };
```
+
We leave the namespace Properties.
+
```cpp
}
```
-## The problem class
+
+### The problem class
We enter the problem class where all necessary boundary conditions and initial conditions are set for our simulation.
As this is a porous medium problem, we inherit from the basic PorousMediumFlowProblem.
+
```cpp
template <class TypeTag >
class PointSourceProblem : public PorousMediumFlowProblem<TypeTag>
{
```
+
We use convenient declarations that we derive from the property system.
+
```cpp
using ParentType = PorousMediumFlowProblem<TypeTag>;
using GridView = GetPropType<TypeTag, Properties::GridView>;
@@ -351,7 +430,9 @@ We use convenient declarations that we derive from the property system.
using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
```
+
We define some indices for convenient use in the problem class:
+
```cpp
enum {
pressureH2OIdx = Indices::pressureIdx,
@@ -363,46 +444,60 @@ We define some indices for convenient use in the problem class:
public:
```
+
This is the constructor of our problem class:
+
```cpp
PointSourceProblem(std::shared_ptr<const GridGeometry> gridGeometry)
: ParentType(gridGeometry)
{
```
+
We read in the values for the initial condition of our simulation:
+
```cpp
initialValues_ = readFileToContainer<std::vector<PrimaryVariables>>("initialsolutioncc.txt");
}
```
+
First, we define the type of boundary conditions depending on location. Two types of boundary conditions
can be specified: Dirichlet or Neumann boundary condition. On a Dirichlet boundary, the values of the
primary variables need to be fixed. On a Neumann boundary condition, values for derivatives need to be fixed.
Mixed boundary conditions (different types for different equations on the same boundary) are not accepted.
+
```cpp
BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
{
BoundaryTypes values;
```
+
We specify Dirichlet boundaries on the left and right hand side of our domain:
+
```cpp
if (onLeftBoundary_(globalPos) || onRightBoundary_(globalPos))
values.setAllDirichlet();
else
```
+
The top and bottom of our domain are Neumann boundaries:
+
```cpp
values.setAllNeumann();
return values;
}
```
+
Second, we specify the values for the Dirichlet boundaries, depending on location. As mentioned,
we need to fix values of our two primary variables: the water pressure
and the Trichlorethene saturation.
+
```cpp
PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
{
```
+
To determine the density of water for a given state, we build a fluid state with the given conditions:
+
```cpp
PrimaryVariables values;
GetPropType<TypeTag, Properties::FluidState> fluidState;
@@ -410,11 +505,15 @@ To determine the density of water for a given state, we build a fluid state with
fluidState.setPressure(waterPhaseIdx, /*pressure=*/1e5);
fluidState.setPressure(dnaplPhaseIdx, /*pressure=*/1e5);
```
+
The density is then calculated by the fluid system:
+
```cpp
Scalar densityW = FluidSystem::density(fluidState, waterPhaseIdx);
```
+
The water phase pressure is the hydrostatic pressure, scaled with a factor:
+
```cpp
Scalar height = this->gridGeometry().bBoxMax()[1] - this->gridGeometry().bBoxMin()[1];
Scalar depth = this->gridGeometry().bBoxMax()[1] - globalPos[1];
@@ -424,26 +523,34 @@ The water phase pressure is the hydrostatic pressure, scaled with a factor:
values[pressureH2OIdx] = 1e5 - factor*densityW*this->spatialParams().gravity(globalPos)[1]*depth;
```
+
The saturation of the DNAPL Trichlorethene is zero on our Dirichlet boundary:
+
```cpp
values[saturationDNAPLIdx] = 0.0;
return values;
}
```
+
Third, we specify the values for the Neumann boundaries.
In our case, we need to specify mass fluxes for our two liquid phases.
The inflow is denoted by a negative sign, and the outflow by a positive sign.
+
```cpp
NumEqVector neumannAtPos(const GlobalPosition &globalPos) const
{
```
+
We initialize the fluxes with zero:
+
```cpp
NumEqVector values(0.0);
```
+
At the inlet, we specify an inflow for our DNAPL Trichlorethene.
The units are $`kg/(m^2 s)`$.
+
```cpp
if (onInlet_(globalPos))
values[contiDNAPLEqIdx] = -0.04;
@@ -451,15 +558,19 @@ The units are $`kg/(m^2 s)`$.
return values;
}
```
+
Last, we specify the initial conditions. The initial condition needs to be set for all primary variables.
Here, we take the data from the file that we read in previously.
+
```cpp
PrimaryVariables initial(const Element& element) const
{
```
+
The input data is written for a uniform grid with discretization length delta.
Accordingly, we need to find the index of our cells, depending on the x and y coordinates,
that corresponds to the indices of the input data set.
+
```cpp
const auto delta = 0.0625;
unsigned int cellsX = this->gridGeometry().bBoxMax()[0]/delta;
@@ -469,26 +580,32 @@ that corresponds to the indices of the input data set.
return initialValues_[dataIdx];
}
```
+
We need to specify a constant temperature for our isothermal problem.
Fluid properties that depend on temperature will be calculated with this value.
+
```cpp
Scalar temperature() const
{
return 293.15; // 10°C
}
```
+
Additionally, we set a point source. The point source can be solution dependent.
It is specified in form of a vector that contains source values for alle phases and positions in space.
The first entry is a tuple containing the position in space, the second entry contains a tuple with the source (with the unit of $`kg/s`$)
for the phases (first phase is the water phase, the second phase is the DNAPL Trichlorethene phase).
+
```cpp
void addPointSources(std::vector<PointSource>& pointSources) const
{
pointSources.push_back(PointSource({0.502, 3.02}, {0, 0.1}));
}
```
+
We define private global functions that are used to determine if a point in space is on the left, right or upper boundary, or
at the inlet.
+
```cpp
private:
bool onLeftBoundary_(const GlobalPosition &globalPos) const
@@ -513,41 +630,53 @@ at the inlet.
return onUpperBoundary_(globalPos) && 0.5 < lambda && lambda < 2.0/3.0;
}
```
+
Our private global variables are the epsilon value and the vector containing the initial values read from file.
+
```cpp
static constexpr Scalar eps_ = 1e-6;
std::vector<PrimaryVariables> initialValues_;
```
+
This is everything the problem class contains.
+
```cpp
};
```
+
We leave the namespace Dumux here, too.
+
```cpp
}
-
```
+
## The file `main.cc`
-## The main file
+
This is the main file for the 2pinfiltration example. Here we can see the programme sequence and how the system is solved using Newton's method
### Includes
+
```cpp
#include <config.h>
-
```
+
Standard header file for C++, to get time and date information.
+
```cpp
#include <ctime>
```
+
Standard header file for C++, for in- and output.
+
```cpp
#include <iostream>
```
+
Dumux is based on DUNE, the Distributed and Unified Numerics Environment, which provides several grid managers and linear solvers. So we need some includes from that.
+
```cpp
#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/timer.hh>
@@ -555,49 +684,71 @@ Dumux is based on DUNE, the Distributed and Unified Numerics Environment, which
#include <dune/grid/io/file/vtk.hh>
#include <dune/istl/io.hh>
```
+
In Dumux, a property system is used to specify the model. For this, different properties are defined containing type definitions, values and methods. All properties are declared in the file properties.hh.
+
```cpp
#include <dumux/common/properties.hh>
```
+
The following file contains the parameter class, which manages the definition of input parameters by a default value, the inputfile or the command line.
+
```cpp
#include <dumux/common/parameters.hh>
```
+
The file dumuxmessage.hh contains the class defining the start and end message of the simulation.
+
```cpp
#include <dumux/common/dumuxmessage.hh>
#include <dumux/common/defaultusagemessage.hh>
```
+
we include the linear solver to be used to solve the linear system
+
```cpp
#include <dumux/linear/amgbackend.hh>
#include <dumux/linear/linearsolvertraits.hh>
```
+
we include the nonlinear Newton's method
+
```cpp
#include <dumux/nonlinear/newtonsolver.hh>
```
+
Further, we include assembler, which assembles the linear systems for finite volume schemes (box-scheme, tpfa-approximation, mpfa-approximation).
+
```cpp
#include <dumux/assembly/fvassembler.hh>
```
+
The containing class in the following file defines the different differentiation methods used to compute the derivatives of the residual.
+
```cpp
#include <dumux/assembly/diffmethod.hh>
```
+
we include the spatial discretization methods available in Dumux
+
```cpp
#include <dumux/discretization/method.hh>
```
+
We need the following class to simplify the writing of dumux simulation data to VTK format.
+
```cpp
#include <dumux/io/vtkoutputmodule.hh>
```
+
The gridmanager constructs a grid from the information in the input or grid file. There is a specification for the different supported grid managers.
+
```cpp
#include <dumux/io/grid/gridmanager.hh>
```
+
We include several files which are needed for the adaptive grid
+
```cpp
#include <dumux/adaptive/adapt.hh>
#include <dumux/adaptive/markelements.hh>
@@ -605,167 +756,231 @@ We include several files which are needed for the adaptive grid
#include <dumux/porousmediumflow/2p/griddatatransfer.hh>
#include <dumux/porousmediumflow/2p/gridadaptindicator.hh>
```
+
We include the problem file which defines initial and boundary conditions to describe our example problem
+
```cpp
#include "problem.hh"
```
+
### Beginning of the main function
+
```cpp
int main(int argc, char** argv) try
{
using namespace Dumux;
```
+
we define the type tag for this problem
+
```cpp
using TypeTag = Properties::TTag::PointSourceExample;
```
+
We initialize MPI, finalize is done automatically on exit
+
```cpp
const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
```
+
We print dumux start message
+
```cpp
if (mpiHelper.rank() == 0)
DumuxMessage::print(/*firstCall=*/true);
```
+
We parse command line arguments and input file
+
```cpp
Parameters::init(argc, argv);
```
+
### Create the grid
A gridmanager tries to create the grid either from a grid file or the input file.
+
```cpp
GridManager<GetPropType<TypeTag, Properties::Grid>> gridManager;
gridManager.init();
```
+
The instationary non-linear problem is run on this grid.
we compute on the leaf grid view
+
```cpp
const auto& leafGridView = gridManager.grid().leafGridView();
```
+
### Set-up and solving of the problem
#### Set-up
We create and initialize the finite volume grid geometry, the problem, the linear system, including the jacobian matrix, the residual and the solution vector and the gridvariables.
We need the finite volume geometry to build up the subcontrolvolumes (scv) and subcontrolvolume faces (scvf) for each element of the grid partition.
+
```cpp
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
auto gridGeometry = std::make_shared<GridGeometry>(leafGridView);
gridGeometry->update();
```
+
In the problem, we define the boundary and initial conditions.
+
```cpp
using Problem = GetPropType<TypeTag, Properties::Problem>;
auto problem = std::make_shared<Problem>(gridGeometry);
```
+
We call the `computePointSourceMap` method to compute the point sources. The `computePointSourceMap` method is inherited from the fvproblem and therefore specified in the `dumux/common/fvproblem.hh`. It calls the `addPointSources` method specified in the `problem.hh` file.
+
```cpp
problem->computePointSourceMap();
```
+
We initialize the solution vector,
+
```cpp
using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
SolutionVector x(gridGeometry->numDofs());
problem->applyInitialSolution(x);
auto xOld = x;
```
+
and then use the solution vector to intialize the `gridVariables`.
+
```cpp
using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
auto gridVariables = std::make_shared<GridVariables>(problem, gridGeometry);
gridVariables->init(x);
```
+
We instantiate the indicator for grid adaption & the data transfer, we read some parameters for indicator from the input file.
+
```cpp
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
const Scalar refineTol = getParam<Scalar>("Adaptive.RefineTolerance");
const Scalar coarsenTol = getParam<Scalar>("Adaptive.CoarsenTolerance");
```
+
We use an indicator for a two-phase flow problem that is saturation-dependent and defined in the file `dumux/porousmediumflow/2p/gridadaptindicator.hh.` It allows to set the minimum and maximum allowed refinement levels via the input parameters.
+
```cpp
TwoPGridAdaptIndicator<TypeTag> indicator(gridGeometry);
```
+
The data transfer performs the transfer of data on a grid from before to after adaptation and is defined in the file `dumux/porousmediumflow/2p/griddatatransfer.hh`. Its main functions are to store and reconstruct the primary variables.
+
```cpp
TwoPGridDataTransfer<TypeTag> dataTransfer(problem, gridGeometry, gridVariables, x);
```
+
We do an initial refinement around sources/BCs. We use the `GridAdaptInitializationIndicator` defined in `dumux/adaptive/initializationindicator.hh` for that.
+
```cpp
GridAdaptInitializationIndicator<TypeTag> initIndicator(problem, gridGeometry, gridVariables);
```
+
We refine up to the maximum level. For every level, the indicator used for the refinement/coarsening is calculated. If any grid cells have to be adapted, the gridvariables and the pointsourcemap are updated.
+
```cpp
const auto maxLevel = getParam<std::size_t>("Adaptive.MaxLevel", 0);
for (std::size_t i = 0; i < maxLevel; ++i)
{
```
+
we calculate the initial indicator for adaption for each grid cell using the initial solution x
+
```cpp
initIndicator.calculate(x);
```
+
and then we mark the elements that were adapted.
+
```cpp
bool wasAdapted = false;
if (markElements(gridManager.grid(), initIndicator))
wasAdapted = adapt(gridManager.grid(), dataTransfer);
```
+
In case of a grid adaptation, the gridvariables and the pointsourcemap are updated.
+
```cpp
if (wasAdapted)
{
```
+
We overwrite the old solution with the new (resized & interpolated) one
+
```cpp
xOld = x;
```
+
We initialize the secondary variables to the new (and "new old") solution
+
```cpp
gridVariables->updateAfterGridAdaption(x);
```
+
we update the point source map after adaption
+
```cpp
problem->computePointSourceMap();
}
}
```
+
Depending on the initial conditions, another grid adaptation might be necessary. The gridadaptindicator uses the input parameters `Adaptive.RefineTolerance` and `Adaptive.CoarsenTolerance` for this step.
+
```cpp
indicator.calculate(x, refineTol, coarsenTol);
```
+
we mark the elements that were adapted
+
```cpp
bool wasAdapted = false;
if (markElements(gridManager.grid(), indicator))
wasAdapted = adapt(gridManager.grid(), dataTransfer);
```
+
In case of an additional grid adaptation, the gridvariables and the pointsourcemap are updated again.
+
```cpp
if (wasAdapted)
{
```
+
Overwrite the old solution with the new (resized & interpolated) one
+
```cpp
xOld = x;
```
+
Initialize the secondary variables to the new (and "new old") solution
+
```cpp
gridVariables->updateAfterGridAdaption(x);
```
+
Update the point source map
+
```cpp
problem->computePointSourceMap();
}
```
+
We get some time loop parameters from the input file params.input
+
```cpp
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
auto dt = getParam<Scalar>("TimeLoop.DtInitial");
```
+
and initialize the vtkoutput. Each model has a predefined model specific output with relevant parameters for that model.
+
```cpp
using IOFields = GetPropType<TypeTag, Properties::IOFields>;
VtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables, x, problem->name());
@@ -774,102 +989,146 @@ and initialize the vtkoutput. Each model has a predefined model specific output
IOFields::initOutputModule(vtkWriter); // Add model specific output fields
vtkWriter.write(0.0);
```
+
we instantiate the time loop
+
```cpp
auto timeLoop = std::make_shared<TimeLoop<Scalar>>(0, dt, tEnd);
timeLoop->setMaxTimeStepSize(maxDt);
```
+
we set the assembler with the time loop because we have an instationary problem
+
```cpp
using Assembler = FVAssembler<TypeTag, DiffMethod::numeric>;
auto assembler = std::make_shared<Assembler>(problem, gridGeometry, gridVariables, timeLoop);
```
+
we set the linear solver
+
```cpp
using LinearSolver = AMGBiCGSTABBackend<LinearSolverTraits<GridGeometry>>;
auto linearSolver = std::make_shared<LinearSolver>(leafGridView, gridGeometry->dofMapper());
```
+
additionally we set the non-linear solver.
+
```cpp
using NewtonSolver = Dumux::NewtonSolver<Assembler, LinearSolver>;
NewtonSolver nonLinearSolver(assembler, linearSolver);
```
+
we start the time loop. In each time step before we start calculating a new solution we check if we have to refine the mesh again based on the solution of the previous time step.
+
```cpp
timeLoop->start(); do
{
```
+
We only want to refine/coarsen after first time step is finished, not before. The initial refinement was already done before the start of the time loop. This means we only refine when the time is greater than 0.
+
```cpp
if (timeLoop->time() > 0)
{
```
+
again we compute the refinement indicator with the `TwoPGridAdaptIndicator`
+
```cpp
indicator.calculate(x, refineTol, coarsenTol);
```
+
we mark elements and adapt grid if necessary
+
```cpp
wasAdapted = false;
if (markElements(gridManager.grid(), indicator))
wasAdapted = adapt(gridManager.grid(), dataTransfer);
```
+
In case of a grid adaptation, the gridvariables and the pointsourcemap are updated again.
+
```cpp
if (wasAdapted)
{
```
+
We overwrite the old solution with the new (resized & interpolated) one
+
```cpp
xOld = x;
```
+
We tell the assembler to resize the matrix and set pattern
+
```cpp
assembler->setJacobianPattern();
```
+
We tell the assembler to resize the residual
+
```cpp
assembler->setResidualSize();
```
+
We initialize the secondary variables to the new (and "new old") solution
+
```cpp
gridVariables->updateAfterGridAdaption(x);
```
+
We update the point source map
+
```cpp
problem->computePointSourceMap();
}
```
+
we leave the refinement step
+
```cpp
}
```
+
Now, we start to calculate the new solution of that time step. First, we define the old solution as the solution of the previous time step for storage evaluations.
+
```cpp
assembler->setPreviousSolution(xOld);
```
+
We solve the non-linear system with time step control.
+
```cpp
nonLinearSolver.solve(x, *timeLoop);
```
+
We make the new solution the old solution.
+
```cpp
xOld = x;
gridVariables->advanceTimeStep();
```
+
We advance to the time loop to the next step.
+
```cpp
timeLoop->advanceTimeStep();
```
+
We write vtk output for each time step
+
```cpp
vtkWriter.write(timeLoop->time());
```
+
We report statistics of this time step
+
```cpp
timeLoop->reportTimeStep();
```
+
We set a new dt as suggested by the newton solver for the next time step
+
```cpp
timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()));
@@ -877,9 +1136,11 @@ We set a new dt as suggested by the newton solver for the next time step
timeLoop->finalize(leafGridView.comm());
```
+
### Final Output
print dumux end message
+
```cpp
if (mpiHelper.rank() == 0)
{
@@ -915,6 +1176,7 @@ catch (...)
}
```
+
## Results
The 2p-infiltration model computes the water saturation distribution after 500 seconds according tthe following figure:
diff --git a/examples/freeflowchannel/README.md b/examples/freeflowchannel/README.md
index 2159c8619bf64768e3e3c1bbc669a99d340ab7e6..c06fd91176746329c21569312356b9ad5c42a2ae 100644
--- a/examples/freeflowchannel/README.md
+++ b/examples/freeflowchannel/README.md
@@ -24,12 +24,12 @@ In the following, we take a close look at the files containing the set-up: At fi
## The file `problem.hh`
-
### The problem class
We enter the problem class where all necessary initial and boundary conditions are set for our simulation.
As this is a Stokes problem, we inherit from the basic <code>NavierStokesProblem</code>.
<details><summary>Toggle to expand code:</summary>
+
```cpp
#include <dumux/freeflow/navierstokes/problem.hh>
@@ -40,12 +40,14 @@ template <class TypeTag>
class ChannelExampleProblem : public NavierStokesProblem<TypeTag>
{
```
+
</details>
We use convenient declarations that we derive from the property system.
<details>
<summary>Toggle to expand code (convenient declarations)</summary>
+
```cpp
using ParentType = NavierStokesProblem<TypeTag>;
using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
@@ -62,6 +64,7 @@ We use convenient declarations that we derive from the property system.
public:
```
+
</details>
There follows the constructor of our problem class:
@@ -70,6 +73,7 @@ As no run-time value is specified, we set the outlet pressure to 1.1e5 Pa.
<details>
<summary>Toggle to expand code (constructor)</summary>
+
```cpp
ChannelExampleProblem(std::shared_ptr<const GridGeometry> gridGeometry)
: ParentType(gridGeometry)
@@ -78,6 +82,7 @@ As no run-time value is specified, we set the outlet pressure to 1.1e5 Pa.
outletPressure_ = getParam<Scalar>("Problem.OutletPressure", 1.1e5);
}
```
+
</details>
Now, we define the type of initial and boundary conditions depending on location.
@@ -94,6 +99,7 @@ of our domain else.
<details>
<summary>Toggle to expand code (<code>boundaryTypesAtPos</code>)</summary>
+
```cpp
BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
{
@@ -117,6 +123,7 @@ of our domain else.
return values;
}
```
+
</details>
Second, we specify the values for the Dirichlet boundaries. We need to fix the values of our primary variables.
@@ -125,6 +132,7 @@ in x-direction is set to zero if not at the inlet.
<details>
<summary>Toggle to expand code (<code>dirichletAtPos</code>)</summary>
+
```cpp
PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
{
@@ -138,6 +146,7 @@ in x-direction is set to zero if not at the inlet.
return values;
}
```
+
</details>
We specify the values for the initial conditions.
@@ -145,6 +154,7 @@ We assign constant values for pressure and velocity components.
<details>
<summary>Toggle to expand code (<code>initialAtPos</code>)</summary>
+
```cpp
PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
{
@@ -157,6 +167,7 @@ We assign constant values for pressure and velocity components.
return values;
}
```
+
</details>
We need to specify a constant temperature for our isothermal problem.
@@ -164,53 +175,61 @@ We set it to 10°C.
<details>
<summary>Toggle to expand code (<code>temperature</code>)</summary>
+
```cpp
Scalar temperature() const
{ return 273.15 + 10; }
private:
```
+
</details>
The inlet is at the left side of the physical domain.
<details>
<summary>Toggle to expand code (<code>isInlet_</code>)</summary>
+
```cpp
bool isInlet_(const GlobalPosition& globalPos) const
{
return globalPos[0] < eps_;
}
```
+
</details>
The outlet is at the right side of the physical domain.
<details>
<summary>Toggle to expand code (<code>isOutlet_</code>)</summary>
+
```cpp
bool isOutlet_(const GlobalPosition& globalPos) const
{
return globalPos[0] > this->gridGeometry().bBoxMax()[0] - eps_;
}
```
+
</details>
Finally, private variables are declared:
<details>
<summary>Toggle to expand code (private variables)</summary>
+
```cpp
static constexpr Scalar eps_=1e-6;
Scalar inletVelocity_;
Scalar outletPressure_;
};
}
-#endif
```
+
</details>
+
## The file `main.cc`
@@ -226,12 +245,14 @@ and another standard header for in- and output.
<details>
<summary>Toggle to expand code (includes of problem file and of standard headers)</summary>
+
```cpp
#include <config.h>
#include <ctime>
#include <iostream>
```
+
</details>
Dumux is based on DUNE, the Distributed and Unified Numerics Environment, which provides several grid managers and
@@ -239,6 +260,7 @@ linear solvers. So we need some includes from that.
<details>
<summary>Toggle to expand code (dune includes)</summary>
+
```cpp
#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/timer.hh>
@@ -246,6 +268,7 @@ linear solvers. So we need some includes from that.
#include <dune/grid/io/file/vtk.hh>
#include <dune/istl/io.hh>
```
+
</details>
In Dumux, a property system is used to specify the model. For this, different properties are defined containing
@@ -267,6 +290,7 @@ The following class contains functionality for additional flux output to the con
<details>
<summary>Toggle to expand code (dumux includes)</summary>
+
```cpp
#include <dumux/common/properties.hh>
#include <dumux/common/parameters.hh>
@@ -290,6 +314,7 @@ The following class contains functionality for additional flux output to the con
#include "problem.hh"
```
+
</details>
</details>
@@ -302,6 +327,7 @@ We setup the DuMux properties for our simulation (click [here](https://git.iws.u
4. The problem class `ChannelExampleProblem`, which is forward declared before we enter `namespace Dumux` and defined later in this file, is defined to be the problem used in this test problem (charaterized by the TypeTag `ChannelExample`). The fluid system, which contains information about the properties such as density, viscosity or diffusion coefficient of the fluid we're simulating, is set to a constant one phase liquid.
5. We enable caching for the following classes (which stores values that were already calculated for later usage and thus results in higher memory usage but improved CPU speed): the grid volume variables, the grid flux variables, the finite volume grid geometry.
+
```cpp
namespace Dumux::Properties {
@@ -334,11 +360,13 @@ struct EnableGridGeometryCache<TypeTag, TTag::ChannelExample> { static constexpr
```
+
### Beginning of the main function
We begin the main function by making the type tag `ChannelExample`, that we defined in `problem.hh` for this test problem available here.
Then we initializing the message passing interface (MPI), even if we do not plan to run the application in parallel. Finalizing of the MPI is done automatically on exit.
We continue by printing the dumux start message and parsing the command line arguments and runtimeparameters from the input file in the init function.
+
```cpp
int main(int argc, char** argv) try
{
@@ -354,6 +382,7 @@ int main(int argc, char** argv) try
Parameters::init(argc, argv);
```
+
### Set-up and solving of the problem
A gridmanager tries to create the grid either from a grid file or the input file. You can learn more about grids in
@@ -370,6 +399,7 @@ primary variables (velocities, pressures) as well as secondary variables (densit
We then initialize the vtkoutput. Each model has a predefined model-specific output with relevant parameters
for that model. Here, it is pressure, velocity, density and process rank (relevant in the case of parallelisation).
+
```cpp
GridManager<GetPropType<TypeTag, Properties::Grid>> gridManager;
gridManager.init();
@@ -399,6 +429,7 @@ for that model. Here, it is pressure, velocity, density and process rank (releva
vtkWriter.write(0.0);
```
+
We set up two surfaces over which fluxes are calculated.
We determine the extensions [xMin,xMax]x[yMin,yMax] of the physical domain.
The first surface (added by the first call of addSurface) shall be placed at the middle of the channel.
@@ -408,6 +439,7 @@ In this case, we add half a cell-width to the x-position in order to make sure t
the cell faces lie on the surface. This assumes a regular cartesian grid.
The second surface (second call of addSurface) is placed at the outlet of the channel.
+
```cpp
FluxOverSurface<GridVariables,
SolutionVector,
@@ -442,6 +474,7 @@ The second surface (second call of addSurface) is placed at the outlet of the ch
flux.addSurface("outlet", p0outlet, p1outlet);
```
+
The incompressible Stokes equation depends only linearly on the velocity, so we could use a linear solver to solve the problem.
Here, we use the show the more general case which would also work for incompressible fluids or the
Navier-Stokes equation. We use non-linear Newton solver for the solution.
@@ -451,6 +484,7 @@ As a postprocessing, we calculate mass and volume fluxes over the planes specifi
<details>
<summary>Toggle to expand code (assembly, solution process, postprocessing)</summary>
+
```cpp
using Assembler = StaggeredFVAssembler<TypeTag, DiffMethod::numeric>;
auto assembler = std::make_shared<Assembler>(problem, gridGeometry, gridVariables);
@@ -466,6 +500,7 @@ As a postprocessing, we calculate mass and volume fluxes over the planes specifi
flux.calculateMassOrMoleFluxes();
flux.calculateVolumeFluxes();
```
+
</details>
### Final Output
@@ -475,6 +510,7 @@ possibly catched error messages are printed.
<details>
<summary>Toggle to expand code (final output)</summary>
+
```cpp
vtkWriter.write(1.0);
@@ -525,6 +561,7 @@ catch (...)
return 4;
}
```
+
</details>
diff --git a/examples/shallowwaterfriction/README.md b/examples/shallowwaterfriction/README.md
index 6d4c51c877bd6eac834cada02c9758d69462ead2..c59eded330c24598808e562bbe5410e2c04cc264 100644
--- a/examples/shallowwaterfriction/README.md
+++ b/examples/shallowwaterfriction/README.md
@@ -73,24 +73,33 @@ for the simulation are given in the file *params.input*.
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
@@ -98,7 +107,9 @@ class RoughChannelSpatialParams
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>;
@@ -110,7 +121,9 @@ We introduce using declarations that are derived from the property system which
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)
@@ -121,7 +134,9 @@ In the constructor be read some values from the `params.input` and initialize th
initFrictionLaw();
}
```
+
We initialize the friction law based on the law specified in `params.input`.
+
```cpp
void initFrictionLaw()
{
@@ -141,21 +156,27 @@ We initialize the friction law based on the law specified in `params.input`.
}
}
```
+
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
@@ -163,7 +184,9 @@ This function returns an object of the friction law class, which is initialized
return *frictionLaw_;
}
```
+
Define the bed surface based on the `bedSlope_`.
+
```cpp
Scalar bedSurface(const Element& element,
const SubControlVolume& scv) const
@@ -178,95 +201,126 @@ private:
std::unique_ptr<FrictionLaw<VolumeVariables>> frictionLaw_;
};
```
+
end of namespace Dumux.
+
```cpp
}
-
```
+
## The file `problem.hh`
## Include files
We use the dune yasp grid.
+
```cpp
#include <dune/grid/yaspgrid.hh>
```
+
We include the cell centered, two-point-flux discretization scheme.
+
```cpp
#include <dumux/discretization/cctpfa.hh>
```
+
The parameters header is needed to retrieve run-time parameters.
+
```cpp
#include <dumux/common/parameters.hh>
```
+
We include the header which are needed for shallow water models.
+
```cpp
#include <dumux/freeflow/shallowwater/model.hh>
#include <dumux/freeflow/shallowwater/problem.hh>
#include <dumux/freeflow/shallowwater/boundaryfluxes.hh>
```
+
We include the header that specifies all spatially variable parameters.
+
```cpp
#include "spatialparams.hh"
```
+
## Define basic properties for our simulation
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. One could use these functions and classes by prefixing every use of these names by ::, but that would quickly become cumbersome and annoying. Rather, we simply import the entire Dumux namespace for general use.
+
```cpp
namespace Dumux {
```
+
The problem class is forward declared.
+
```cpp
template <class TypeTag>
class RoughChannelProblem;
```
+
We enter the namespace Properties, which is a sub-namespace of the namespace Dumux.
+
```cpp
namespace Properties {
```
+
A TypeTag for our simulation is created which inherits from the shallow water model and the
cell centered, two-point-flux discretization scheme.
+
```cpp
namespace TTag {
struct RoughChannel { using InheritsFrom = std::tuple<ShallowWater, CCTpfaModel>; };
}
```
+
We define the grid of our simulation. We use a two-dimensional Yasp Grid.
+
```cpp
template<class TypeTag>
struct Grid<TypeTag, TTag::RoughChannel>
{ using type = Dune::YaspGrid<2, Dune::TensorProductCoordinates<GetPropType<TypeTag, Properties::Scalar>, 2> >; };
```
+
We set the problem. The problem class specifies initial and boundary conditions and is defined below.
+
```cpp
template<class TypeTag>
struct Problem<TypeTag, TTag::RoughChannel>
{ using type = Dumux::RoughChannelProblem<TypeTag>; };
```
+
We define the spatial parameters for our simulation. The values are specified in the corresponding spatialparameters header file, which is included above.
+
```cpp
template<class TypeTag>
struct SpatialParams<TypeTag, TTag::RoughChannel>
{
private:
```
+
We define convenient shortcuts to the properties GridGeometry, Scalar, ElementVolumeVariables and VolumeVariables:
+
```cpp
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;
```
+
Finally we set the spatial parameters:
+
```cpp
public:
using type = RoughChannelSpatialParams<GridGeometry, Scalar, VolumeVariables>;
};
```
+
We enable caching for the FV grid geometry and the grid volume variables. The cache
stores values that were already calculated for later usage. This makes the simulation faster.
+
```cpp
template<class TypeTag>
struct EnableGridGeometryCache<TypeTag, TTag::RoughChannel>
@@ -276,19 +330,25 @@ template<class TypeTag>
struct EnableGridVolumeVariablesCache<TypeTag, TTag::RoughChannel>
{ static constexpr bool value = false; };
```
+
We leave the namespace Properties.
+
```cpp
}
```
+
## The problem class
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
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>;
@@ -311,44 +371,58 @@ We use convenient declarations that we derive from the property system.
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)
{
@@ -359,7 +433,9 @@ Get the water depth with Gauckler-Manning-Strickler
return pow(abs(discharge)*manningN/sqrt(bedSlope), 0.6);
}
```
+
Get the analytical solution
+
```cpp
void analyticalSolution()
{
@@ -376,14 +452,18 @@ Get the analytical solution
}
}
```
+
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,
@@ -393,14 +473,18 @@ Get the source term.
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,
@@ -410,11 +494,15 @@ Get the source term due to bottom friction.
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];
@@ -423,19 +511,29 @@ The bottom shear stress causes a pure loss of momentum. Thus the first entry of
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).
+
+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,
@@ -451,7 +549,9 @@ We specify the neumann boundary. Due to the weak imposition we calculate the flu
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_)
{
@@ -463,7 +563,9 @@ Calculate the rieman invariants for imposed discharge at the left side.
nxy);
}
```
+
Calculate the rieman invariants for impose water depth at the right side.
+
```cpp
else if (scvf.center()[0] > 100.0 - eps_)
{
@@ -475,7 +577,9 @@ Calculate the rieman invariants for impose water depth at the right side.
nxy);
}
```
+
Calculate the rieman invarianty for the no-flow boundary.
+
```cpp
else
{
@@ -484,7 +588,9 @@ Calculate the rieman invarianty for the no-flow boundary.
boundaryStateVariables[2] = -insideVolVars.velocity(1);
}
```
+
We calculate the boundary fluxes based on a Riemann problem.
+
```cpp
auto riemannFlux = ShallowWater::riemannProblem(insideVolVars.waterDepth(),
boundaryStateVariables[0],
@@ -504,268 +610,373 @@ We calculate the boundary fluxes based on a Riemann problem.
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.
+
```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_;
};
```
+
We leave the namespace Dumux.
+
```cpp
}
-
```
+
## The file `main.cc`
-## The main file
+
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 problem file which defines initial and boundary conditions to describe our example problem
+
```cpp
#include "problem.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);
```
+
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()));
@@ -774,8 +985,10 @@ We set new dt as suggested by newton controller for the next time step.
timeLoop->finalize(leafGridView.comm());
```
+
### Final Output
We print dumux end message.
+
```cpp
if (mpiHelper.rank() == 0)
{
@@ -812,5 +1025,6 @@ catch (...)
}
```
+
## Results
The solution and the analytical result for the problem is shown in .