From a770a49e19d859592d57f1a75eb46f0952f519a0 Mon Sep 17 00:00:00 2001
From: "Dennis.Glaeser"
Date: Thu, 19 Mar 2020 14:06:16 +0100
Subject: [PATCH 1/9] [examples][tracer] edits on intro.md

examples/1ptracer/doc/intro.md  27 ++++++++++++++
1 file changed, 14 insertions(+), 13 deletions()
diff git a/examples/1ptracer/doc/intro.md b/examples/1ptracer/doc/intro.md
index 3d95a6063c..22461959bc 100644
 a/examples/1ptracer/doc/intro.md
+++ b/examples/1ptracer/doc/intro.md
@@ 8,43 +8,44 @@ This example contains a contaminant transported by a base groundwater flow in a
![](./img/setup.png)
## Model description
Two different models are applied to simulate the system: In a first step, the groundwater velocity is evaluated under stationary conditions. Therefore the single phase model is applied.
In a second step, the contaminant gets transported based on the groundwater velocity field. It is assumed, that the dissolved contaminant does not affect density and viscosity of the groundwater and thus, it is handled as a tracer by the tracer model. The tracer model is then solved instationarily.
+Two different models are applied to simulate the system: In a first step, the groundwater velocity is evaluated under stationary conditions using the single phase model.
+In a second step, the contaminant is transported with the groundwater velocity field. It is assumed, that the dissolved contaminant does not affect density and viscosity of the groundwater, and thus, it is handled as a tracer by the tracer model. The tracer model is then solved instationarily.
### 1p Model
The single phase model uses Darcy's law as the equation for the momentum conservation:
```math
\textbf v =  \frac{\textbf K}{\mu} \left(\textbf{grad}\, p  \varrho {\textbf g} \right)
+\textbf v =  \frac{\textbf K}{\mu} \left(\textbf{grad}\, p  \varrho {\textbf g} \right),
```
With the darcy velocity $` \textbf v `$, the permeability $` \textbf K`$, the dynamic viscosity $` \mu`$, the pressure $`p`$, the density $`\rho`$ and the gravity $`\textbf g`$.
+with the darcy velocity $` \textbf v `$, the permeability $` \textbf K`$, the dynamic viscosity $` \mu`$, the pressure $`p`$, the density $`\rho`$ and the gravity $`\textbf g`$.
Darcy's law is inserted into the continuity equation:
+Darcy's law is inserted into the mass balance equation:
```math
\phi \frac{\partial \varrho}{\partial t} + \text{div} \textbf v = 0
```
with porosity $`\phi`$.
+where $`\phi`$ is the porosity.
The equation is discretized using a cellcentered finite volume scheme as spatial discretization for the pressure as primary variable. For details on the discretization scheme, have a look at the dumux [handbook](https://dumux.org/handbook).
### Tracer Model
The transport of the contaminant component $`\kappa`$ is based on the previously evaluated velocity field $`\textbf v`$ with the help the following mass balance equation:
+The transport of the contaminant component $`\kappa`$ is based on the previously evaluated velocity field $`\textbf v`$ with the help of the following mass balance equation:
```math
\phi \frac{ \partial \varrho X^\kappa}{\partial t}  \text{div} \left\lbrace \varrho X^\kappa {\textbf v} + \varrho D^\kappa_\text{pm} \textbf{grad} X^\kappa \right\rbrace = 0
+\phi \frac{ \partial \varrho X^\kappa}{\partial t}  \text{div} \left\lbrace \varrho X^\kappa {\textbf v} + \varrho D^\kappa_\text{pm} \textbf{grad} X^\kappa \right\rbrace = 0,
```
With the porosity $`\phi`$, the mass fraction of the contaminant component $`\kappa`$: $`X^\kappa`$, the porous medium diffusivity $` D^\kappa_\text{pm} `$ and the density of the fluid phase $`\varrho`$.
+where $`X^\kappa`$ is the mass fraction of the contaminant component $`\kappa`$ and $` D^\kappa_\text{pm} `$ is the effective diffusivity.
The porous medium diffusivity is a function of the diffusion coefficient of the component $`D^\kappa`$, the porosity $`\phi`$ and the porous medium tortuosity $`\tau`$ by the following equation:
+The effective diffusivity is a function of the diffusion coefficient of the component $`D^\kappa`$ and the porosity and tortuosity $`\tau`$ of the porous medium:
```math
D^\kappa_\text{pm}= \phi \tau D^\kappa
+D^\kappa_\text{pm}= \phi \tau D^\kappa.
```
The primary variable of this model is the mass fraction $`X^\kappa`$. We apply the same spatial discretization as in the single phase model and use the implicit Euler method for time discretization. For more information, have a look at the dumux handbook.
+The primary variable of this model is the mass fraction $`X^\kappa`$. We apply the same spatial discretization as in the single phase model and use the implicit Euler method for time discretization. For more information, have a look at the dumux [handbook](https://dumux.org/handbook).
In the following, we take a close look at the files containing the setup: At first, boundary conditions and spatially distributed parameters are set in `problem_1p.hh` and `spatialparams_1p.hh`, respectively, for the single phase model and subsequently in `problem_tracer.hh` and `spatialparams_tracer.hh` for the tracer model. Afterwards, we show the different steps for solving the model in the source file `main.cc`. At the end, we show some simulation results.
+In the following, we take a close look at the files containing the setup: The boundary conditions and spatially distributed parameters for the single phase model are set in `problem_1p.hh` and `spatialparams_1p.hh`.
+For the tracer model, this is done in the files `problem_tracer.hh` and `spatialparams_tracer.hh`, respectively. Afterwards, we show the different steps for solving the model in the source file `main.cc`. Finally, some simulation results are shown.

GitLab
From c14cb62351f2251a324a0c59f980ac2b5e7cb661 Mon Sep 17 00:00:00 2001
From: "Dennis.Glaeser"
Date: Thu, 19 Mar 2020 14:06:47 +0100
Subject: [PATCH 2/9] [examples][tracer] edits on spatialparams_1p

examples/1ptracer/spatialparams_1p.hh  29 +++++++++++++++++
1 file changed, 18 insertions(+), 11 deletions()
diff git a/examples/1ptracer/spatialparams_1p.hh b/examples/1ptracer/spatialparams_1p.hh
index 38d22f0fba..191f144383 100644
 a/examples/1ptracer/spatialparams_1p.hh
+++ b/examples/1ptracer/spatialparams_1p.hh
@@ 21,25 +21,26 @@
#ifndef DUMUX_INCOMPRESSIBLE_ONEP_TEST_SPATIAL_PARAMS_HH
#define DUMUX_INCOMPRESSIBLE_ONEP_TEST_SPATIAL_PARAMS_HH
// In this file, we generate the random permeability field in the constructor of the `OnePTestSpatialParams` class. Thereafter, spatial properties of the porous medium such as the permeability and the porosity are defined in various functions for the 1p problem.
// We want to generate a random permeability field. For this, we use a random number generation of the C++ standard library.
+// 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.
#include
// In the file `properties.hh` all properties are declared.
+
+// We use the properties for porous medium flow models, declared in the file `properties.hh`.
#include
// We include the spatial parameters for singlephase, finite volumes from which we will inherit.
+// We include the spatial parameters class for singlephase models discretized by finite volume schemes.
+// The spatial parameters defined for this example will inherit from those.
#include
namespace Dumux {
// In the `OnePTestSpatialParams` class, we define all functions needed to describe the porous matrix, e.g. porosity and permeability for the 1p_problem.

+// In the `OnePTestSpatialParams` class, we define all functions needed to describe the porous medium, e.g. porosity and permeability for the 1p_problem.
template
class OnePTestSpatialParams
: public FVSpatialParamsOneP>
{
 // We introduce `using` declarations that are derived from the property system, which we need in this class.
+ // We declare aliases for types that we are going to need in this class.
using GridView = typename GridGeometry::GridView;
using FVElementGeometry = typename GridGeometry::LocalView;
using SubControlVolume = typename FVElementGeometry::SubControlVolume;
@@ 51,6 +52,9 @@ class OnePTestSpatialParams
using GlobalPosition = typename SubControlVolume::GlobalPosition;
public:
+
+ // The spatial parameters must export the type used to define permeabilities.
+ // Here, we are using scalar permeabilities, but tensors are also supported.
using PermeabilityType = Scalar;
OnePTestSpatialParams(std::shared_ptr gridGeometry)
: ParentType(gridGeometry), K_(gridGeometry>gridView().size(0), 0.0)
@@ 60,11 +64,12 @@ public:
permeability_ = getParam("SpatialParams.Permeability");
permeabilityLens_ = getParam("SpatialParams.PermeabilityLens");
 // Further, we get the position of the lens, which is defined by the position of the lower left and the upper right corner.
+ // Furthermore, we get the position of the lens, which is defined by the position of the lower left and the upper right corner.
lensLowerLeft_ = getParam("SpatialParams.LensLowerLeft");
lensUpperRight_ =getParam("SpatialParams.LensUpperRight");
 // We generate random fields for the permeability using a lognormal distribution, with the `permeability_` as mean value and 10 % of it as standard deviation. A seperate distribution is used for the lens using `permeabilityLens_`.
+ // 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_`.
std::mt19937 rand(0);
std::lognormal_distribution K(std::log(permeability_), std::log(permeability_)*0.1);
std::lognormal_distribution KLens(std::log(permeabilityLens_), std::log(permeabilityLens_)*0.1);
@@ 77,7 +82,9 @@ public:
}
// ### Properties of the porous matrix
 // We define the (intrinsic) permeability $`[m^2]`$ using the generated random permeability field. In this test, we use elementwise distributed permeabilities.
+ // This function returns the permeability $`[m^2]`$ to be used within a subcontrol 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 elementwise distributed permeabilities that were randomly generated in the constructor (see above).
template
const PermeabilityType& permeability(const Element& element,
const SubControlVolume& scv,
@@ 87,7 +94,7 @@ public:
}
 // We set the porosity $`[]`$ for the whole domain.
+ // We set the porosity $`[]`$ for the whole domain to a value of $`20 \%`$.
Scalar porosityAtPos(const GlobalPosition &globalPos) const
{ return 0.2; }

GitLab
From 193ae32e13ef22fed577b9d56961dbbfd63b100a Mon Sep 17 00:00:00 2001
From: "Dennis.Glaeser"
Date: Thu, 19 Mar 2020 14:07:04 +0100
Subject: [PATCH 3/9] [examples][tracer] edits on problem_1p

examples/1ptracer/problem_1p.hh  31 ++++++++++++++++
1 file changed, 16 insertions(+), 15 deletions()
diff git a/examples/1ptracer/problem_1p.hh b/examples/1ptracer/problem_1p.hh
index b3c8d67478..e9fae345a1 100644
 a/examples/1ptracer/problem_1p.hh
+++ b/examples/1ptracer/problem_1p.hh
@@ 19,12 +19,12 @@
// ### Header guard
#ifndef DUMUX_ONEP_TRACER_TEST_PROBLEM_HH
#define DUMUX_ONEP_TRACER_TEST_PROBLEM_HH
+#ifndef DUMUX_ONEP_TEST_PROBLEM_HH
+#define DUMUX_ONEP_TEST_PROBLEM_HH
//Before we enter the problem class containing initial and boundary conditions, we include necessary files and introduce properties.
+// Before we enter the problem class containing initial and boundary conditions, we include necessary files and introduce properties.
// ### Include files
// The dune grid interphase is included here:
+// We use `YaspGrid`, an implementation of the dune grid interface for structured grids:
#include
// The cell centered, twopointflux discretization scheme is included:
@@ 34,7 +34,7 @@
// This is the porous medium problem class that this class is derived from:
#include
// The fluid properties are specified in the following headers:
+// The fluid properties are specified in the following headers (we use liquid water as the fluid phase):
#include
#include
@@ 78,11 +78,11 @@ struct SpatialParams
using type = OnePTestSpatialParams;
};
// The local residual contains analytic derivative methods for incompressible flow:
+// We use the local residual that contains analytic derivative methods for incompressible flow:
template
struct LocalResidual { using type = OnePIncompressibleLocalResidual; };
 // In the following we define our fluid properties.
+// In the following we define the fluid properties.
template
struct FluidSystem
{
@@ 102,13 +102,13 @@ struct EnableGridFluxVariablesCache { static
// We enable caching for the FV grid geometry
template
struct EnableGridGeometryCache { static constexpr bool value = true; };
//The cache stores values that were already calculated for later usage. This makes the simulation faster.
+// 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.
}
// ### 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`.
+// As this is a porous medium problem, we inherit from the base class `PorousMediumFlowProblem`.
template
class OnePTestProblem : public PorousMediumFlowProblem
{
@@ 130,24 +130,25 @@ public:
OnePTestProblem(std::shared_ptr gridGeometry)
: ParentType(gridGeometry) {}
 // First, we define the type of boundary conditions depending on location. Two types of boundary conditions
+ // 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.
+ // Mixed boundary conditions (different types for different equations on the same boundary) are not accepted for
+ // cellcentered finite volume schemes.
BoundaryTypes boundaryTypes(const Element &element,
const SubControlVolumeFace &scvf) const
{
BoundaryTypes values;
 // we retreive the global position, i.e. the vector including the global coordinates
 // of the finite volume
+ // we retrieve the global position, i.e. the vector with the global coordinates,
+ // of the integration point on the boundary subcontrol volume face `scvf`
const auto globalPos = scvf.ipGlobal();
// we define a small epsilon value
Scalar eps = 1.0e6;
// We specify Dirichlet boundaries on the top and bottom of our domain:
if (globalPos[dimWorld1] < eps  globalPos[dimWorld1] > this>gridGeometry().bBoxMax()[dimWorld1]  eps)
values.setAllDirichlet();
+ // The top and bottom of our domain are Neumann boundaries:
else
 // The top and bottom of our domain are Neumann boundaries:
values.setAllNeumann();
return values;
@@ 160,7 +161,7 @@ public:
// we retreive again the global position
const auto& pos = scvf.ipGlobal();
PrimaryVariables values(0);
 // we assign pressure values in [Pa] according to a pressure gradient to 1e5 Pa at the top and 1.1e5 Pa at the bottom.
+ // and assign pressure values in [Pa] according to a pressure gradient to 1e5 Pa at the top and 1.1e5 Pa at the bottom.
values[0] = 1.0e+5*(1.1  pos[dimWorld1]*0.1);
return values;
}

GitLab
From 43945255a285d81501739d8c38c7c9f72880c1f2 Mon Sep 17 00:00:00 2001
From: "Dennis.Glaeser"
Date: Thu, 19 Mar 2020 14:07:28 +0100
Subject: [PATCH 4/9] [examples][tracer] edits on spatialparams_tracer

examples/1ptracer/problem_tracer.hh  3 ++
examples/1ptracer/spatialparams_tracer.hh  26 ++++++++++++++
2 files changed, 18 insertions(+), 11 deletions()
diff git a/examples/1ptracer/problem_tracer.hh b/examples/1ptracer/problem_tracer.hh
index 62f03c09a0..7f95bf8a33 100644
 a/examples/1ptracer/problem_tracer.hh
+++ b/examples/1ptracer/problem_tracer.hh
@@ 202,7 +202,8 @@ public:
if (useMoles)
initialValues = 1e9;
else
 initialValues = 1e9*FluidSystem::molarMass(0)/this>spatialParams().fluidMolarMass(globalPos);
+ initialValues = 1e9*FluidSystem::molarMass(0)
+ /this>spatialParams().fluidMolarMassAtPos(globalPos);
}
return initialValues;
}
diff git a/examples/1ptracer/spatialparams_tracer.hh b/examples/1ptracer/spatialparams_tracer.hh
index 00bedcae4c..7159034943 100644
 a/examples/1ptracer/spatialparams_tracer.hh
+++ b/examples/1ptracer/spatialparams_tracer.hh
@@ 23,16 +23,16 @@
#ifndef DUMUX_TRACER_TEST_SPATIAL_PARAMS_HH
#define DUMUX_TRACER_TEST_SPATIAL_PARAMS_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. Further, spatial dependent properties of the tracer fluid system are defined and in the end two functions handel the calculated volume fluxes from the solution of the 1p problem.
// In the file `properties.hh`, all properties are declared.
+// 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`.
#include
// As in the 1p spatialparams, we inherit from the spatial parameters for singlephase, finite volumes, which we include here.
+// As in the 1p spatialparams, we inherit from the spatial parameters for singlephase models using finite volumes, which we include here.
#include
// We enter the namespace Dumux
namespace Dumux {
// In the `TracerTestSpatialParams` class, we define all functions needed to describe spatially dependent parameters for the `tracer_problem`.

template
class TracerTestSpatialParams
: public FVSpatialParamsOneP
Scalar volumeFlux(const Element &element,
const FVElementGeometry& fvGeometry,
@@ 92,7 +97,8 @@ public:
return volumeFlux_[scvf.index()];
}
 // We define a function to set the volume flux. This is used in the main function to set the volume flux to the calculated value based on the solution of the 1p problem.
+ // We define a function that allows setting the volume fluxes for all subcontrol 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 singlephase model.
void setVolumeFlux(const std::vector& f)
{ volumeFlux_ = f; }

GitLab
From 2c62c05935afe2201c83f8dd6dc1a6d9022d2e0b Mon Sep 17 00:00:00 2001
From: "Dennis.Glaeser"
Date: Thu, 19 Mar 2020 14:08:27 +0100
Subject: [PATCH 5/9] [examples][tracer] edits on problem_tracer

examples/1ptracer/problem_tracer.hh  104 +++++++++++++++
1 file changed, 57 insertions(+), 47 deletions()
diff git a/examples/1ptracer/problem_tracer.hh b/examples/1ptracer/problem_tracer.hh
index 7f95bf8a33..804d5ae6bf 100644
 a/examples/1ptracer/problem_tracer.hh
+++ b/examples/1ptracer/problem_tracer.hh
@@ 23,7 +23,7 @@
//Before we enter the problem class containing initial and boundary conditions, we include necessary files and introduce properties.
// ### Include files
// Again, we have to include the dune grid interface:
+// Again, we use YaspGrid, the implementation of the dune grid interface for structured grids:
#include
// and the cell centered, twopointflux discretization.
#include
@@ 31,7 +31,7 @@
#include
// We include again the porous medium problem class that this class is derived from:
#include
// and the base fluidsystem
+// and the base fluidsystem. We will define a custom fluid system that inherits from that class.
#include
// We include the header that specifies all spatially variable parameters for the tracer problem:
#include "spatialparams_tracer.hh"
@@ 73,25 +73,29 @@ struct Problem { using type = TracerTestProblem
struct SpatialParams
{
+private:
using GridGeometry = GetPropType;
using Scalar = GetPropType;
+public:
using type = TracerTestSpatialParams;
};
// We define that mass fractions are used to define the concentrations
+// One can choose between a formulation in terms of mass or mole fractions. Here, we are using mass fractions.
template
struct UseMoles { static constexpr bool value = false; };
// We do not use a solution dependent molecular diffusion coefficient:
+// We use solutionindependent molecular diffusion coefficients. Per default, solutiondependent
+// diffusion coefficients are assumed during the computation of the jacobian matrix entries. Specifying
+// solutionindependent diffusion coefficients can speed up computations:
template
struct SolutionDependentMolecularDiffusion { static constexpr bool value = false; };
+struct SolutionDependentMolecularDiffusion
+{ static constexpr bool value = false; };
// In the following, we create a new tracer fluid system and derive it from the base fluid system.
+// In the following, we create a new tracer fluid system and derive from the base fluid system.
template
class TracerFluidSystem : public FluidSystems::Base,
TracerFluidSystem>
{
 // We define convenient shortcuts to the properties `Scalar`, `Problem`, `GridView`,
 // `Element`, `FVElementGeometry` and `SubControlVolume`:
+ // We define some convenience aliases:
using Scalar = GetPropType;
using Problem = GetPropType;
using GridView = GetPropType;
@@ 100,31 +104,35 @@ class TracerFluidSystem : public FluidSystems::Base { using type = TracerFluidSystem
class TracerTestProblem : public PorousMediumFlowProblem
{
@@ 164,9 +172,9 @@ class TracerTestProblem : public PorousMediumFlowProblem
using ElementFluxVariablesCache = typename GridVariables::GridFluxVariablesCache::LocalView;
using FVElementGeometry = typename GetPropType::LocalView;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
 //We create a bool saying whether mole or mass fractions are used
+ // We create a convenience bool stating whether mole or mass fractions are used
static constexpr bool useMoles = getPropValue();
 //We create additional int to make dimWorld and numComponents available in the problem
+ // We create additional convenience integers to make dimWorld and numComponents available in the problem
static constexpr int dimWorld = GridView::dimensionworld;
static const int numComponents = FluidSystem::numComponents;
@@ 175,7 +183,7 @@ public:
TracerTestProblem(std::shared_ptr gridGeometry)
: ParentType(gridGeometry)
{
 // We print out whether mole or mass fractions are used
+ // We print to the terminal whether mole or mass fractions are used
if(useMoles)
std::cout<<"problem uses mole fractions" << '\n';
else
@@ 195,10 +203,12 @@ public:
PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
{
PrimaryVariables initialValues(0.0);
 // The tracer concentration is located on the domain bottom:
+
+ // The initial contamination is located at the bottom of the domain:
if (globalPos[1] < 0.1 + eps_)
{
 // We assign different values, depending on wether mole concentrations or mass concentrations are used:
+ // We chose a mole fraction of $`1e9`$, but in case the mass fractions
+ // are used by the model, we have to convert this value:
if (useMoles)
initialValues = 1e9;
else
@@ 208,33 +218,33 @@ public:
return initialValues;
}
 //We implement an outflow boundary on the top of the domain and prescribe zeroflux Neumann boundary conditions on all other boundaries.
 NumEqVector neumann(const Element& element,
 const FVElementGeometry& fvGeometry,
 const ElementVolumeVariables& elemVolVars,
 const ElementFluxVariablesCache& elemFluxVarsCache,
 const SubControlVolumeFace& scvf) const
 {
 NumEqVector values(0.0);
 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 and diffusive flux is enforced to be zero
 if (globalPos[dimWorld1] > this>gridGeometry().bBoxMax()[dimWorld1]  eps_)
 {
 values = this>spatialParams().volumeFlux(element, fvGeometry, elemVolVars, scvf)
 * volVars.massFraction(0, 0) * volVars.density(0)
 / scvf.area();
 assert(values>=0.0 && "Volume flux at outflow boundary is expected to have a positive sign");
 }
 //prescribe zeroflux Neumann boundary conditions elsewhere
 else
 values = 0.0;
+ // We implement an outflow boundary on the top of the domain and prescribe zeroflux Neumann boundary conditions on all other boundaries.
+ NumEqVector neumann(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const ElementFluxVariablesCache& elemFluxVarsCache,
+ const SubControlVolumeFace& scvf) const
+ {
+ NumEqVector values(0.0);
+ const auto& volVars = elemVolVars[scvf.insideScvIdx()];
+ const auto& globalPos = scvf.center();
 return values;
+ // This is the outflow boundary, where tracer is transported by advection with the given flux field.
+ if (globalPos[dimWorld1] > this>gridGeometry().bBoxMax()[dimWorld1]  eps_)
+ {
+ values = this>spatialParams().volumeFlux(element, fvGeometry, elemVolVars, scvf)
+ * volVars.massFraction(0, 0) * volVars.density(0)
+ / scvf.area();
+ assert(values>=0.0 && "Volume flux at outflow boundary is expected to have a positive sign");
}
+ // Prescribe zeroflux Neumann boundary conditions elsewhere
+ else
+ values = 0.0;
+
+ return values;
+ }
+
private:
// We assign a private global variable for the epsilon:
static constexpr Scalar eps_ = 1e6;

GitLab
From 7f091910b4a7fd6058d98c159b384dde3b7f89b1 Mon Sep 17 00:00:00 2001
From: "Dennis.Glaeser"
Date: Thu, 19 Mar 2020 14:08:42 +0100
Subject: [PATCH 6/9] [examples][tracer] edits on main file

examples/1ptracer/main.cc  64 +++++++++++++++++++
1 file changed, 32 insertions(+), 32 deletions()
diff git a/examples/1ptracer/main.cc b/examples/1ptracer/main.cc
index 0325eeac46..5ee3a9152d 100644
 a/examples/1ptracer/main.cc
+++ b/examples/1ptracer/main.cc
@@ 30,15 +30,19 @@
// and another one for in and output.
#include
// 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.
+// 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.
#include
#include
#include
#include
// In Dumux, a property system is used to specify the model. For this, different properties are defined containing type definitions, values and methods. All properties are declared in the file `properties.hh`.
+// In Dumux, the property system is used to specify classes and compiletime 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`.
#include
// The following file contains the parameter class, which manages the definition of input parameters by a default value, the inputfile or the command line.
+// 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.
#include
// The file `dumuxmessage.hh` contains the class defining the start and end message of the simulation.
#include
@@ 60,7 +64,7 @@ int main(int argc, char** argv) try
{
using namespace Dumux;
 // We define the type tags for the two problems. They are created in the individual problem files.
+ // Convenience aliases for the type tags of the two problems, which are defined in the individual problem files.
using OnePTypeTag = Properties::TTag::IncompressibleTest;
using TracerTypeTag = Properties::TTag::TracerTestCC;
@@ 75,7 +79,11 @@ int main(int argc, char** argv) try
Parameters::init(argc, argv);
// ### Create the grid
 // A gridmanager tries to create the grid either from a grid file or the input file. We solve both problems on the same grid. Hence, the grid is only created once using the type tag of the 1p problem.
+ // 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 singlephase 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.
GridManager> gridManager;
gridManager.init();
@@ 83,7 +91,7 @@ int main(int argc, char** argv) try
const auto& leafGridView = gridManager.grid().leafGridView();
// ### Setup 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 problem domain.
+ // 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.
// #### 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.
@@ 103,13 +111,13 @@ int main(int argc, char** argv) try
auto A = std::make_shared();
auto r = std::make_shared();
 // The grid variables store variables on scv and scvf (volume and flux variables).
+ // The grid variables store variables (primary and secondary variables) on subcontrol volumes and faces (volume and flux variables).
using OnePGridVariables = GetPropType;
auto onePGridVariables = std::make_shared(problemOneP, gridGeometry);
onePGridVariables>init(p);
// #### Assembling the linear system
 // We created and inizialize the assembler.
+ // We create and inizialize the assembler.
using OnePAssembler = FVAssembler;
auto assemblerOneP = std::make_shared(problemOneP, gridGeometry, onePGridVariables);
assemblerOneP>setLinearSystem(A, r);
@@ 157,16 +165,17 @@ int main(int argc, char** argv) try
// ### Computation of the volume fluxes
// We use the results of the 1p problem to calculate the volume fluxes in the model domain.

using Scalar = GetPropType;
std::vector volumeFlux(gridGeometry>numScvf(), 0.0);
using FluxVariables = GetPropType;
auto upwindTerm = [](const auto& volVars) { return volVars.mobility(0); };
 // We iterate over all elements.
+ // We iterate over all elements
for (const auto& element : elements(leafGridView))
{
+ // Compute the elementlocal views on geometry, primary and secondary variables
+ // as well as variables needed for flux computations
auto fvGeometry = localView(*gridGeometry);
fvGeometry.bind(element);
@@ 176,28 +185,17 @@ int main(int argc, char** argv) try
auto elemFluxVars = localView(onePGridVariables>gridFluxVarsCache());
elemFluxVars.bind(element, fvGeometry, elemVolVars);
 // We calculate the volume flux for every subcontrolvolume face, which is not on a Neumann boundary (is not on the boundary or is on a Dirichlet boundary).

+ // We calculate the volume fluxes for all subcontrol volume faces except for Neumann boundary faces
for (const auto& scvf : scvfs(fvGeometry))
{
 const auto idx = scvf.index();

 if (!scvf.boundary())
 {
 FluxVariables fluxVars;
 fluxVars.init(*problemOneP, element, fvGeometry, elemVolVars, scvf, elemFluxVars);
 volumeFlux[idx] = fluxVars.advectiveFlux(0, upwindTerm);
 }
 else
 {
 const auto bcTypes = problemOneP>boundaryTypes(element, scvf);
 if (bcTypes.hasOnlyDirichlet())
 {
 FluxVariables fluxVars;
 fluxVars.init(*problemOneP, element, fvGeometry, elemVolVars, scvf, elemFluxVars);
 volumeFlux[idx] = fluxVars.advectiveFlux(0, upwindTerm);
 }
 }
+ // skip Neumann boundary faces
+ if (scvf.boundary() && problemOneP>boundaryTypes(element, scvf).hasNeumann())
+ continue;
+
+ // let the `FluxVariables` class do the flux computation.
+ FluxVariables fluxVars;
+ fluxVars.init(*problemOneP, element, fvGeometry, elemVolVars, scvf, elemFluxVars);
+ volumeFlux[scvf.index()] = fluxVars.advectiveFlux(0, upwindTerm);
}
}
@@ 244,7 +242,7 @@ int main(int argc, char** argv) try
vtkWriter.write(0.0);
 // For the time loop we set a check point.
+ // We define 10 check points in the time loop at which we will write the solution to vtk files.
timeLoop>setPeriodicCheckPoint(tEnd/10.0);
// #### The time loop
@@ 302,13 +300,15 @@ int main(int argc, char** argv) try
// ### Final Output

if (mpiHelper.rank() == 0)
DumuxMessage::print(/*firstCall=*/false);
return 0;
}
+// ### Exception handling
+// In this part of the main file we catch and print possible exceptions that could
+// occur during the simulation.
catch (Dumux::ParameterException &e)
{
std::cerr << std::endl << e << " > Abort!" << std::endl;

GitLab
From af65683c2df34c4619af933cabcba17b549fb3fe Mon Sep 17 00:00:00 2001
From: "Dennis.Glaeser"
Date: Thu, 19 Mar 2020 14:19:01 +0100
Subject: [PATCH 7/9] [examples][tracer] generate new readme from changes

examples/1ptracer/README.md  282 ++++++++++++++++++++
1 file changed, 157 insertions(+), 125 deletions()
diff git a/examples/1ptracer/README.md b/examples/1ptracer/README.md
index a2e1d56de2..c85a34d1ec 100644
 a/examples/1ptracer/README.md
+++ b/examples/1ptracer/README.md
@@ 8,76 +8,77 @@ This example contains a contaminant transported by a base groundwater flow in a
![](./img/setup.png)
## Model description
Two different models are applied to simulate the system: In a first step, the groundwater velocity is evaluated under stationary conditions. Therefore the single phase model is applied.
In a second step, the contaminant gets transported based on the groundwater velocity field. It is assumed, that the dissolved contaminant does not affect density and viscosity of the groundwater and thus, it is handled as a tracer by the tracer model. The tracer model is then solved instationarily.
+Two different models are applied to simulate the system: In a first step, the groundwater velocity is evaluated under stationary conditions using the single phase model.
+In a second step, the contaminant is transported with the groundwater velocity field. It is assumed, that the dissolved contaminant does not affect density and viscosity of the groundwater, and thus, it is handled as a tracer by the tracer model. The tracer model is then solved instationarily.
### 1p Model
The single phase model uses Darcy's law as the equation for the momentum conservation:
```math
\textbf v =  \frac{\textbf K}{\mu} \left(\textbf{grad}\, p  \varrho {\textbf g} \right)
+\textbf v =  \frac{\textbf K}{\mu} \left(\textbf{grad}\, p  \varrho {\textbf g} \right),
```
With the darcy velocity $` \textbf v `$, the permeability $` \textbf K`$, the dynamic viscosity $` \mu`$, the pressure $`p`$, the density $`\rho`$ and the gravity $`\textbf g`$.
+with the darcy velocity $` \textbf v `$, the permeability $` \textbf K`$, the dynamic viscosity $` \mu`$, the pressure $`p`$, the density $`\rho`$ and the gravity $`\textbf g`$.
Darcy's law is inserted into the continuity equation:
+Darcy's law is inserted into the mass balance equation:
```math
\phi \frac{\partial \varrho}{\partial t} + \text{div} \textbf v = 0
```
with porosity $`\phi`$.
+where $`\phi`$ is the porosity.
The equation is discretized using a cellcentered finite volume scheme as spatial discretization for the pressure as primary variable. For details on the discretization scheme, have a look at the dumux [handbook](https://dumux.org/handbook).
### Tracer Model
The transport of the contaminant component $`\kappa`$ is based on the previously evaluated velocity field $`\textbf v`$ with the help the following mass balance equation:
+The transport of the contaminant component $`\kappa`$ is based on the previously evaluated velocity field $`\textbf v`$ with the help of the following mass balance equation:
```math
\phi \frac{ \partial \varrho X^\kappa}{\partial t}  \text{div} \left\lbrace \varrho X^\kappa {\textbf v} + \varrho D^\kappa_\text{pm} \textbf{grad} X^\kappa \right\rbrace = 0
+\phi \frac{ \partial \varrho X^\kappa}{\partial t}  \text{div} \left\lbrace \varrho X^\kappa {\textbf v} + \varrho D^\kappa_\text{pm} \textbf{grad} X^\kappa \right\rbrace = 0,
```
With the porosity $`\phi`$, the mass fraction of the contaminant component $`\kappa`$: $`X^\kappa`$, the porous medium diffusivity $` D^\kappa_\text{pm} `$ and the density of the fluid phase $`\varrho`$.
+where $`X^\kappa`$ is the mass fraction of the contaminant component $`\kappa`$ and $` D^\kappa_\text{pm} `$ is the effective diffusivity.
The porous medium diffusivity is a function of the diffusion coefficient of the component $`D^\kappa`$, the porosity $`\phi`$ and the porous medium tortuosity $`\tau`$ by the following equation:
+The effective diffusivity is a function of the diffusion coefficient of the component $`D^\kappa`$ and the porosity and tortuosity $`\tau`$ of the porous medium:
```math
D^\kappa_\text{pm}= \phi \tau D^\kappa
+D^\kappa_\text{pm}= \phi \tau D^\kappa.
```
The primary variable of this model is the mass fraction $`X^\kappa`$. We apply the same spatial discretization as in the single phase model and use the implicit Euler method for time discretization. For more information, have a look at the dumux handbook.
+The primary variable of this model is the mass fraction $`X^\kappa`$. We apply the same spatial discretization as in the single phase model and use the implicit Euler method for time discretization. For more information, have a look at the dumux [handbook](https://dumux.org/handbook).
In the following, we take a close look at the files containing the setup: At first, boundary conditions and spatially distributed parameters are set in `problem_1p.hh` and `spatialparams_1p.hh`, respectively, for the single phase model and subsequently in `problem_tracer.hh` and `spatialparams_tracer.hh` for the tracer model. Afterwards, we show the different steps for solving the model in the source file `main.cc`. At the end, we show some simulation results.
+In the following, we take a close look at the files containing the setup: The boundary conditions and spatially distributed parameters for the single phase model are set in `problem_1p.hh` and `spatialparams_1p.hh`.
+For the tracer model, this is done in the files `problem_tracer.hh` and `spatialparams_tracer.hh`, respectively. Afterwards, we show the different steps for solving the model in the source file `main.cc`. Finally, some simulation results are shown.
## The file `spatialparams_1p.hh`
In this file, we generate the random permeability field in the constructor of the `OnePTestSpatialParams` class. Thereafter, spatial properties of the porous medium such as the permeability and the porosity are defined in various functions for the 1p problem.
We want to generate a random permeability field. For this, we use a random number generation of the C++ standard library.
+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
```
In the file `properties.hh` all properties are declared.
+We use the properties for porous medium flow models, declared in the file `properties.hh`.
```cpp
#include
```
We include the spatial parameters for singlephase, finite volumes from which we will inherit.
+We include the spatial parameters class for singlephase models discretized by finite volume schemes.
+The spatial parameters defined for this example will inherit from those.
```cpp
#include
namespace Dumux {
```
In the `OnePTestSpatialParams` class, we define all functions needed to describe the porous matrix, e.g. porosity and permeability for the 1p_problem.
+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 OnePTestSpatialParams
: public FVSpatialParamsOneP>
{
```
We introduce `using` declarations that are derived from the property system, which we need in this class.
+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,6 +91,10 @@ We introduce `using` declarations that are derived from the property system, whi
using GlobalPosition = typename SubControlVolume::GlobalPosition;
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 gridGeometry)
: ParentType(gridGeometry), K_(gridGeometry>gridView().size(0), 0.0)
@@ 101,12 +106,13 @@ We get the permeability of the domain and the lens from the `params.input` file.
permeability_ = getParam("SpatialParams.Permeability");
permeabilityLens_ = getParam("SpatialParams.PermeabilityLens");
```
Further, we get the position of the lens, which is defined by the position of the lower left and the upper right corner.
+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("SpatialParams.LensLowerLeft");
lensUpperRight_ =getParam("SpatialParams.LensUpperRight");
```
We generate random fields for the permeability using a lognormal distribution, with the `permeability_` as mean value and 10 % of it as standard deviation. A seperate distribution is used for the lens using `permeabilityLens_`.
+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 K(std::log(permeability_), std::log(permeability_)*0.1);
@@ 120,7 +126,9 @@ We generate random fields for the permeability using a lognormal distribution, w
}
```
### Properties of the porous matrix
We define the (intrinsic) permeability $`[m^2]`$ using the generated random permeability field. In this test, we use elementwise distributed permeabilities.
+This function returns the permeability $`[m^2]`$ to be used within a subcontrol 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 elementwise distributed permeabilities that were randomly generated in the constructor (see above).
```cpp
template
const PermeabilityType& permeability(const Element& element,
@@ 131,7 +139,7 @@ We define the (intrinsic) permeability $`[m^2]`$ using the generated random perm
}
```
We set the porosity $`[]`$ for the whole domain.
+We set the porosity $`[]`$ for the whole domain to a value of $`20 \%`$.
```cpp
Scalar porosityAtPos(const GlobalPosition &globalPos) const
{ return 0.2; }
@@ 175,7 +183,7 @@ We have a convenient definition of the position of the lens.
Before we enter the problem class containing initial and boundary conditions, we include necessary files and introduce properties.
### Include files
The dune grid interphase is included here:
+We use `YaspGrid`, an implementation of the dune grid interface for structured grids:
```cpp
#include
```
@@ 191,7 +199,7 @@ This is the porous medium problem class that this class is derived from:
```cpp
#include
```
The fluid properties are specified in the following headers:
+The fluid properties are specified in the following headers (we use liquid water as the fluid phase):
```cpp
#include
#include
@@ 250,12 +258,12 @@ Finally, we set the spatial parameters:
using type = OnePTestSpatialParams;
};
```
The local residual contains analytic derivative methods for incompressible flow:
+We use the local residual that contains analytic derivative methods for incompressible flow:
```cpp
template
struct LocalResidual { using type = OnePIncompressibleLocalResidual; };
```
In the following we define our fluid properties.
+In the following we define the fluid properties.
```cpp
template
struct FluidSystem
@@ 286,14 +294,14 @@ We enable caching for the FV grid geometry
template
struct EnableGridGeometryCache { static constexpr bool value = true; };
```
The cache stores values that were already calculated for later usage. This makes the simulation faster.
+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 basic `PorousMediumFlowProblem`.
+As this is a porous medium problem, we inherit from the base class `PorousMediumFlowProblem`.
```cpp
template
class OnePTestProblem : public PorousMediumFlowProblem
@@ 320,18 +328,19 @@ This is the constructor of our problem class:
OnePTestProblem(std::shared_ptr gridGeometry)
: ParentType(gridGeometry) {}
```
First, we define the type of boundary conditions depending on location. Two types of boundary conditions
+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.
+Mixed boundary conditions (different types for different equations on the same boundary) are not accepted for
+cellcentered finite volume schemes.
```cpp
BoundaryTypes boundaryTypes(const Element &element,
const SubControlVolumeFace &scvf) const
{
BoundaryTypes values;
```
we retreive the global position, i.e. the vector including the global coordinates
of the finite volume
+we retrieve the global position, i.e. the vector with the global coordinates,
+of the integration point on the boundary subcontrol volume face `scvf`
```cpp
const auto globalPos = scvf.ipGlobal();
```
@@ 343,10 +352,10 @@ We specify Dirichlet boundaries on the top and bottom of our domain:
```cpp
if (globalPos[dimWorld1] < eps  globalPos[dimWorld1] > this>gridGeometry().bBoxMax()[dimWorld1]  eps)
values.setAllDirichlet();
 else
```
The top and bottom of our domain are Neumann boundaries:
```cpp
+ else
values.setAllNeumann();
return values;
@@ 363,7 +372,7 @@ we retreive again the global position
const auto& pos = scvf.ipGlobal();
PrimaryVariables values(0);
```
we assign pressure values in [Pa] according to a pressure gradient to 1e5 Pa at the top and 1.1e5 Pa at the bottom.
+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[dimWorld1]*0.1);
return values;
@@ 391,12 +400,13 @@ We leave the 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. Further, spatial dependent properties of the tracer fluid system are defined and in the end two functions handel the calculated volume fluxes from the solution of the 1p problem.
In the file `properties.hh`, all properties are declared.
+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
```
As in the 1p spatialparams, we inherit from the spatial parameters for singlephase, finite volumes, which we include here.
+As in the 1p spatialparams, we inherit from the spatial parameters for singlephase models using finite volumes, which we include here.
```cpp
#include
```
@@ 406,7 +416,6 @@ namespace Dumux {
```
In the `TracerTestSpatialParams` class, we define all functions needed to describe spatially dependent parameters for the `tracer_problem`.
```cpp

template
class TracerTestSpatialParams
: public FVSpatialParamsOneP
Scalar volumeFlux(const Element &element,
@@ 471,7 +486,8 @@ We define a function which returns the field of volume fluxes. This is e.g. used
return volumeFlux_[scvf.index()];
}
```
We define a function to set the volume flux. This is used in the main function to set the volume flux to the calculated value based on the solution of the 1p problem.
+We define a function that allows setting the volume fluxes for all subcontrol 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 singlephase model.
```cpp
void setVolumeFlux(const std::vector& f)
{ volumeFlux_ = f; }
@@ 491,7 +507,7 @@ private:
Before we enter the problem class containing initial and boundary conditions, we include necessary files and introduce properties.
### Include files
Again, we have to include the dune grid interface:
+Again, we use YaspGrid, the implementation of the dune grid interface for structured grids:
```cpp
#include
```
@@ 507,7 +523,7 @@ We include again the porous medium problem class that this class is derived from
```cpp
#include
```
and the base fluidsystem
+and the base fluidsystem. We will define a custom fluid system that inherits from that class.
```cpp
#include
```
@@ 561,30 +577,34 @@ We define the spatial parameters for our tracer simulation:
template
struct SpatialParams
{
+private:
using GridGeometry = GetPropType;
using Scalar = GetPropType;
+public:
using type = TracerTestSpatialParams;
};
```
We define that mass fractions are used to define the concentrations
+One can choose between a formulation in terms of mass or mole fractions. Here, we are using mass fractions.
```cpp
template
struct UseMoles { static constexpr bool value = false; };
```
We do not use a solution dependent molecular diffusion coefficient:
+We use solutionindependent molecular diffusion coefficients. Per default, solutiondependent
+diffusion coefficients are assumed during the computation of the jacobian matrix entries. Specifying
+solutionindependent diffusion coefficients can speed up computations:
```cpp
template
struct SolutionDependentMolecularDiffusion { static constexpr bool value = false; };
+struct SolutionDependentMolecularDiffusion
+{ static constexpr bool value = false; };
```
In the following, we create a new tracer fluid system and derive it from the base fluid system.
+In the following, we create a new tracer fluid system and derive from the base fluid system.
```cpp
template
class TracerFluidSystem : public FluidSystems::Base,
TracerFluidSystem>
{
```
We define convenient shortcuts to the properties `Scalar`, `Problem`, `GridView`,
`Element`, `FVElementGeometry` and `SubControlVolume`:
+We define some convenience aliases:
```cpp
using Scalar = GetPropType;
using Problem = GetPropType;
@@ 595,37 +615,41 @@ We define convenient shortcuts to the properties `Scalar`, `Problem`, `GridView`
public:
```
We specify, that the fluid system only contains tracer components,
+We specify that the fluid system only contains tracer components,
```cpp
static constexpr bool isTracerFluidSystem()
{ return true; }
```
that no component is the main component
+and that no component is the main component
```cpp
static constexpr int getMainComponent(int phaseIdx)
{ return 1; }
```
and the number of components
+We define the number of components of this fluid system (one single tracer component)
```cpp
static constexpr int numComponents = 1;
```
We set the component name for the component index (`compIdx`) for the vtk output:
+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)
+ 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`.
+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)
+ 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 for our case it is 0.0:
+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,
@@ 646,7 +670,7 @@ We leave the namespace Properties.
```
### 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`.
+As this is a porous medium problem, we inherit from the base class `PorousMediumFlowProblem`.
```cpp
template
class TracerTestProblem : public PorousMediumFlowProblem
@@ 672,11 +696,11 @@ We use convenient declarations that we derive from the property system.
using FVElementGeometry = typename GetPropType::LocalView;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
```
We create a bool saying whether mole or mass fractions are used
+We create a convenience bool stating whether mole or mass fractions are used
```cpp
static constexpr bool useMoles = getPropValue();
```
We create additional int to make dimWorld and numComponents available in the problem
+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;
@@ 689,7 +713,7 @@ This is the constructor of our problem class:
: ParentType(gridGeometry)
{
```
We print out whether mole or mass fractions are used
+We print to the terminal whether mole or mass fractions are used
```cpp
if(useMoles)
std::cout<<"problem uses mole fractions" << '\n';
@@ 713,51 +737,52 @@ We specify the initial conditions for the primary variable (tracer concentration
{
PrimaryVariables initialValues(0.0);
```
The tracer concentration is located on the domain bottom:
+The initial contamination is located at the bottom of the domain:
```cpp
if (globalPos[1] < 0.1 + eps_)
{
```
We assign different values, depending on wether mole concentrations or mass concentrations are used:
+We chose a mole fraction of $`1e9`$, but in case the mass fractions
+are used by the model, we have to convert this value:
```cpp
if (useMoles)
initialValues = 1e9;
else
 initialValues = 1e9*FluidSystem::molarMass(0)/this>spatialParams().fluidMolarMass(globalPos);
+ initialValues = 1e9*FluidSystem::molarMass(0)
+ /this>spatialParams().fluidMolarMassAtPos(globalPos);
}
return initialValues;
}
```
We implement an outflow boundary on the top of the domain and prescribe zeroflux Neumann boundary conditions on all other boundaries.
```cpp
 NumEqVector neumann(const Element& element,
 const FVElementGeometry& fvGeometry,
 const ElementVolumeVariables& elemVolVars,
 const ElementFluxVariablesCache& elemFluxVarsCache,
 const SubControlVolumeFace& scvf) const
 {
 NumEqVector values(0.0);
 const auto& volVars = elemVolVars[scvf.insideScvIdx()];
 const auto& globalPos = scvf.center();
+ NumEqVector neumann(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const ElementFluxVariablesCache& elemFluxVarsCache,
+ const SubControlVolumeFace& scvf) const
+ {
+ NumEqVector values(0.0);
+ 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 and diffusive flux is enforced to be zero
+This is the outflow boundary, where tracer is transported by advection with the given flux field.
```cpp
 if (globalPos[dimWorld1] > this>gridGeometry().bBoxMax()[dimWorld1]  eps_)
 {
 values = this>spatialParams().volumeFlux(element, fvGeometry, elemVolVars, scvf)
 * volVars.massFraction(0, 0) * volVars.density(0)
 / scvf.area();
 assert(values>=0.0 && "Volume flux at outflow boundary is expected to have a positive sign");
 }
+ if (globalPos[dimWorld1] > this>gridGeometry().bBoxMax()[dimWorld1]  eps_)
+ {
+ values = this>spatialParams().volumeFlux(element, fvGeometry, elemVolVars, scvf)
+ * volVars.massFraction(0, 0) * volVars.density(0)
+ / scvf.area();
+ assert(values>=0.0 && "Volume flux at outflow boundary is expected to have a positive sign");
+ }
```
prescribe zeroflux Neumann boundary conditions elsewhere
+Prescribe zeroflux Neumann boundary conditions elsewhere
```cpp
 else
 values = 0.0;
+ else
+ values = 0.0;
 return values;
 }
+ return values;
+ }
private:
```
@@ 798,18 +823,22 @@ and another one for in and output.
```cpp
#include
```
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.
+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
#include
#include
#include
```
In Dumux, a property system is used to specify the model. For this, different properties are defined containing type definitions, values and methods. All properties are declared in the file `properties.hh`.
+In Dumux, the property system is used to specify classes and compiletime 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
```
The following file contains the parameter class, which manages the definition of input parameters by a default value, the inputfile or the command line.
+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
```
@@ 843,7 +872,7 @@ int main(int argc, char** argv) try
{
using namespace Dumux;
```
We define the type tags for the two problems. They are created in the individual problem files.
+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;
@@ 862,7 +891,11 @@ We parse the command line arguments.
Parameters::init(argc, argv);
```
### Create the grid
A gridmanager tries to create the grid either from a grid file or the input file. We solve both problems on the same grid. Hence, the grid is only created once using the type tag of the 1p problem.
+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 singlephase 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> gridManager;
gridManager.init();
@@ 872,7 +905,7 @@ We compute on the leaf grid view.
const auto& leafGridView = gridManager.grid().leafGridView();
```
### Setup 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 problem domain.
+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.
#### 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.
@@ 895,14 +928,14 @@ The jacobian matrix (`A`), the solution vector (`p`) and the residual (`r`) are
auto A = std::make_shared();
auto r = std::make_shared();
```
The grid variables store variables on scv and scvf (volume and flux variables).
+The grid variables store variables (primary and secondary variables) on subcontrol volumes and faces (volume and flux variables).
```cpp
using OnePGridVariables = GetPropType;
auto onePGridVariables = std::make_shared(problemOneP, gridGeometry);
onePGridVariables>init(p);
```
#### Assembling the linear system
We created and inizialize the assembler.
+We create and inizialize the assembler.
```cpp
using OnePAssembler = FVAssembler;
auto assemblerOneP = std::make_shared(problemOneP, gridGeometry, onePGridVariables);
@@ 958,17 +991,20 @@ We stop the timer and display the total time of the simulation as well as the cu
### 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;
std::vector volumeFlux(gridGeometry>numScvf(), 0.0);
using FluxVariables = GetPropType;
auto upwindTerm = [](const auto& volVars) { return volVars.mobility(0); };
```
We iterate over all elements.
+We iterate over all elements
```cpp
for (const auto& element : elements(leafGridView))
{
+```
+Compute the elementlocal views on geometry, primary and secondary variables
+as well as variables needed for flux computations
+```cpp
auto fvGeometry = localView(*gridGeometry);
fvGeometry.bind(element);
@@ 978,29 +1014,21 @@ We iterate over all elements.
auto elemFluxVars = localView(onePGridVariables>gridFluxVarsCache());
elemFluxVars.bind(element, fvGeometry, elemVolVars);
```
We calculate the volume flux for every subcontrolvolume face, which is not on a Neumann boundary (is not on the boundary or is on a Dirichlet boundary).
+We calculate the volume fluxes for all subcontrol volume faces except for Neumann boundary faces
```cpp

for (const auto& scvf : scvfs(fvGeometry))
{
 const auto idx = scvf.index();

 if (!scvf.boundary())
 {
 FluxVariables fluxVars;
 fluxVars.init(*problemOneP, element, fvGeometry, elemVolVars, scvf, elemFluxVars);
 volumeFlux[idx] = fluxVars.advectiveFlux(0, upwindTerm);
 }
 else
 {
 const auto bcTypes = problemOneP>boundaryTypes(element, scvf);
 if (bcTypes.hasOnlyDirichlet())
 {
 FluxVariables fluxVars;
 fluxVars.init(*problemOneP, element, fvGeometry, elemVolVars, scvf, elemFluxVars);
 volumeFlux[idx] = fluxVars.advectiveFlux(0, upwindTerm);
 }
 }
+```
+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);
}
}
@@ 1055,7 +1083,7 @@ We initialize the vtk output module and add a velocity output.
vtkWriter.write(0.0);
```
For the time loop we set a check point.
+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);
```
@@ 1127,13 +1155,17 @@ We set the time step size dt of the next time step.
```
### Final Output
```cpp

if (mpiHelper.rank() == 0)
DumuxMessage::print(/*firstCall=*/false);
return 0;
}
+```
+### 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)
{
std::cerr << std::endl << e << " > Abort!" << std::endl;

GitLab
From 9b0e50102227aec4e27cf2ea7de4ee5f931b9e09 Mon Sep 17 00:00:00 2001
From: "Dennis.Glaeser"
Date: Thu, 19 Mar 2020 14:21:15 +0100
Subject: [PATCH 8/9] [examples][README][tracer] update bullet point list

examples/README.md  7 ++++
1 file changed, 4 insertions(+), 3 deletions()
diff git a/examples/README.md b/examples/README.md
index 29cf75a805..222a5d7886 100644
 a/examples/README.md
+++ b/examples/README.md
@@ 10,9 +10,10 @@ We first solve the pressure field, compute the steady state flow field, and then
You learn how to
* generate a randomly distributed permeability field
* solve a onephase flow in porous media problem
* compute the flow field from a pressure solution to pass to a tracer problem
* sequentially solve two types of problems after each other
+* sequentially solve two types of problems after each other:
+ * solve a onephase flow in porous media problem
+ * compute the flow field from a pressure solution to pass to a tracer problem
+ * solve an instationary tracer transport problem with a given flow field
### [:open_file_folder: Example 2: Twophase flow with infiltration and adaptive grid](https://git.iws.unistuttgart.de/dumuxrepositories/dumux/tree/master/examples/2pinfiltration)

GitLab
From aec715f0345a135a73feca1b74a5634d7dab6169 Mon Sep 17 00:00:00 2001
From: "Dennis.Glaeser"
Date: Thu, 19 Mar 2020 15:37:26 +0100
Subject: [PATCH 9/9] [examples][tracer][intro] link effective diffusivity
model

examples/1ptracer/README.md  2 +
examples/1ptracer/doc/intro.md  2 +
2 files changed, 2 insertions(+), 2 deletions()
diff git a/examples/1ptracer/README.md b/examples/1ptracer/README.md
index c85a34d1ec..b0cd88f3bd 100644
 a/examples/1ptracer/README.md
+++ b/examples/1ptracer/README.md
@@ 39,7 +39,7 @@ The transport of the contaminant component $`\kappa`$ is based on the previously
where $`X^\kappa`$ is the mass fraction of the contaminant component $`\kappa`$ and $` D^\kappa_\text{pm} `$ is the effective diffusivity.
The effective diffusivity is a function of the diffusion coefficient of the component $`D^\kappa`$ and the porosity and tortuosity $`\tau`$ of the porous medium:
+The effective diffusivity is a function of the diffusion coefficient of the component $`D^\kappa`$ and the porosity and tortuosity $`\tau`$ of the porous medium (see [dumux/material/fluidmatrixinteractions/diffusivityconstanttortuosity.hh](https://git.iws.unistuttgart.de/dumuxrepositories/dumux//blob/master/dumux/material/fluidmatrixinteractions/diffusivityconstanttortuosity.hh)):
```math
D^\kappa_\text{pm}= \phi \tau D^\kappa.
diff git a/examples/1ptracer/doc/intro.md b/examples/1ptracer/doc/intro.md
index 22461959bc..fa0964231f 100644
 a/examples/1ptracer/doc/intro.md
+++ b/examples/1ptracer/doc/intro.md
@@ 39,7 +39,7 @@ The transport of the contaminant component $`\kappa`$ is based on the previously
where $`X^\kappa`$ is the mass fraction of the contaminant component $`\kappa`$ and $` D^\kappa_\text{pm} `$ is the effective diffusivity.
The effective diffusivity is a function of the diffusion coefficient of the component $`D^\kappa`$ and the porosity and tortuosity $`\tau`$ of the porous medium:
+The effective diffusivity is a function of the diffusion coefficient of the component $`D^\kappa`$ and the porosity and tortuosity $`\tau`$ of the porous medium (see [dumux/material/fluidmatrixinteractions/diffusivityconstanttortuosity.hh](https://git.iws.unistuttgart.de/dumuxrepositories/dumux//blob/master/dumux/material/fluidmatrixinteractions/diffusivityconstanttortuosity.hh)):
```math
D^\kappa_\text{pm}= \phi \tau D^\kappa.

GitLab