Commit 0c6965d6 authored by Timo Koch's avatar Timo Koch
Browse files

Merge branch 'feature/improve-tracer-example' into 'master'

Feature/improve tracer example

See merge request !1899
parents a38e8857 56b9a826
{
"README.md" : [
"doc/intro.md",
"spatialparams_1p.hh",
"doc/intro.md"
],
"doc/1pmodel.md" : [
"doc/1pintro.md",
"properties_1p.hh",
"problem_1p.hh",
"spatialparams_tracer.hh",
"spatialparams_1p.hh"
],
"doc/tracermodel.md" : [
"doc/tracerintro.md",
"properties_tracer.hh",
"problem_tracer.hh",
"main.cc",
"doc/results.md"
]
"spatialparams_tracer.hh"
],
"doc/main.md" : [
"main.cc"
],
"navigation" : {
"mainpage" : "README.md",
"subpages" : [
"doc/1pmodel.md",
"doc/tracermodel.md",
"doc/main.md"
],
"subtitles" : [
"Single-phase flow simulation setup",
"Tracer transport simulation setup",
"Main program flow"
]
}
}
This diff is collapsed.
# Part 1: Implementation of the single-phase flow simulation setup
The single-phase flow setup is implemented in the files `properties_1p.hh`,
`problem_1p.hh` and `spatialparams_1p.hh`. In the first of these files, a new
type tag is declared for this problem. This then allows the specialization
of DuMuX `properties` for this type tag, which can be used to customize compile-time
settings for the simulation. Two exemplary `properties`, that are mandatory to be
specialized, are `Problem` and `SpatialParams`. With the first, one sets the
`Problem` class to be used, in which users can define initial and boundary conditions.
Similarly, in the `SpatialParams` class one implements the parameter distributions
(e.g. porosity and permeability) that should be used by the model.
The documentation provided in the sequel is structured as follows:
[[_TOC_]]
<!-- Important: This file has been automatically generated by generate_example_docs.py. Do not edit this file directly! -->
| [:arrow_left: Back to the main documentation](../README.md) | [:arrow_right: Continue with part 2](tracermodel.md) |
|---|---:|
# Part 1: Implementation of the single-phase flow simulation setup
The single-phase flow setup is implemented in the files `properties_1p.hh`,
`problem_1p.hh` and `spatialparams_1p.hh`. In the first of these files, a new
type tag is declared for this problem. This then allows the specialization
of DuMuX `properties` for this type tag, which can be used to customize compile-time
settings for the simulation. Two exemplary `properties`, that are mandatory to be
specialized, are `Problem` and `SpatialParams`. With the first, one sets the
`Problem` class to be used, in which users can define initial and boundary conditions.
Similarly, in the `SpatialParams` class one implements the parameter distributions
(e.g. porosity and permeability) that should be used by the model.
The documentation provided in the sequel is structured as follows:
[[_TOC_]]
## Compile-time settings (`properties_1p.hh`)
In this file, the type tag used for the single-phase flow simulation is defined,
for which we then specialize `properties` to the needs of the desired setup.
<details open>
<summary><b>Click to hide/show the file documentation</b> (or inspect the [source code](../properties_1p.hh))</summary>
### Includes
<details><summary> Click to show includes</summary>
The `OneP` type tag specializes most of the `properties` required for single-
phase flow simulations in DuMuX. We will use this in the following to inherit the
respective properties and subsequently specialize those `properties` for our
type tag, which we want to modify or for which no meaningful default can be set.
```cpp
#include <dumux/porousmediumflow/1p/model.hh>
```
We want to use `YaspGrid`, an implementation of the dune grid interface for structured grids:
```cpp
#include <dune/grid/yaspgrid.hh>
```
In this example, we want to discretize the equations with the cell centered finite volume
scheme using two-point-flux approximation:
```cpp
#include <dumux/discretization/cctpfa.hh>
```
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 problem and spatial parameters headers used for this simulation.
```cpp
#include "problem_1p.hh"
#include "spatialparams_1p.hh"
```
</details>
### Type tag definition
We define a type tag for our simulation with the name `IncompressibleTest`
and inherit the `properties` specialized for the type tags `OneP` and `CCTpfaModel`.
This way, most of the `properties` required for single-phase flow simulations
using the cell centered finite volume scheme with two-point-flux approximation
are conveniently specialized for our new type tag.
However, some properties depend on user choices and no meaningful default value can be set.
Those properties will be adressed later in this file.
```cpp
// We enter the namespace Dumux::Properties in order to import the entire Dumux namespace for general use:
namespace Dumux::Properties {
// declaration of the `IncompressibleTest` type tag for the single-phase flow problem
namespace TTag {
struct IncompressibleTest { using InheritsFrom = std::tuple<OneP, CCTpfaModel>; };
}
```
### Property specializations
In the following piece of code, mandatory `properties` for which no meaningful
default can be set, are specialized for our type tag `IncompressibleTest`.
```cpp
// This sets the grid type used for the simulation. Here, we use a structured 2D grid.
template<class TypeTag>
struct Grid<TypeTag, TTag::IncompressibleTest> { using type = Dune::YaspGrid<2>; };
// This sets our problem class (see problem_1p.hh) containing initial and boundary conditions.
template<class TypeTag>
struct Problem<TypeTag, TTag::IncompressibleTest> { using type = OnePTestProblem<TypeTag>; };
// This sets our spatial parameters class (see spatialparams_1p.hh) in which we define
// the distributions for porosity and permeability to be used.
template<class TypeTag>
struct SpatialParams<TypeTag, TTag::IncompressibleTest>
{
private:
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
public:
using type = OnePTestSpatialParams<GridGeometry, Scalar>;
};
// This sets the fluid system to be used. Here, we use liquid water as fluid phase.
template<class TypeTag>
struct FluidSystem<TypeTag, TTag::IncompressibleTest>
{
private:
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
public:
using type = FluidSystems::OnePLiquid<Scalar, Components::SimpleH2O<Scalar> >;
};
```
These are all mandatory `properties`. However, we also want to specialize the
`LocalResidual` property, as we want to use analytic differentation (see `main.cc`).
The default for the `LocalResidual` property, specialized by the type tag `OneP`,
is the implementation of the local residual for compressible single-phase flow.
Therein, the functionality required for analytic derivatives is not implemented
as for compressible fluids this would require the derivatives of the fluid properties
(eg density and viscosity) with respect to pressure. For the case of an incompressible
fluid phase, a specialized local residual is available that contains the required
functionality, and we set this with the following piece of code:
```cpp
template<class TypeTag>
struct LocalResidual<TypeTag, TTag::IncompressibleTest> { using type = OnePIncompressibleLocalResidual<TypeTag>; };
```
Apart from that, we also set some properties related to memory management
throughout the simulation.
<details><summary> Click to show caching properties</summary>
In Dumux, one has the option to activate/deactive the grid-wide caching of
geometries and variables. If active, the CPU time can be significantly reduced
as less dynamic memory allocation procedures are necessary. Per default, grid-wide
caching is disabled to ensure minimal memory requirements, however, in this example we
want to active all available caches, which significanlty increases the memory
demand but makes the simulation faster.
```cpp
// This enables grid-wide caching of the volume variables.
template<class TypeTag>
struct EnableGridVolumeVariablesCache<TypeTag, TTag::IncompressibleTest> { static constexpr bool value = true; };
//This enables grid wide caching for the flux variables.
template<class TypeTag>
struct EnableGridFluxVariablesCache<TypeTag, TTag::IncompressibleTest> { static constexpr bool value = true; };
// This enables grid-wide caching for the finite volume grid geometry
template<class TypeTag>
struct EnableGridGeometryCache<TypeTag, TTag::IncompressibleTest> { static constexpr bool value = true; };
} // end namespace Dumux::Properties
```
</details>
</details>
## Initial and boundary conditions (`problem_1p.hh`)
This file contains the __problem class__ which defines the initial and boundary
conditions for the single-phase flow simulation.
<details open>
<summary><b>Click to hide/show the file documentation</b> (or inspect the [source code](../problem_1p.hh))</summary>
### Include files
The only include we need here is the `PorousMediumFlowProblem` class, the base
class from which we will derive.
```cpp
#include <dumux/porousmediumflow/problem.hh>
```
### The problem class
We enter the problem class where all necessary boundary conditions and initial conditions are set for our simulation.
As we are solving a problem related to flow in porous media, we inherit from the base class `PorousMediumFlowProblem`.
```cpp
namespace Dumux {
template<class TypeTag>
class OnePTestProblem : public PorousMediumFlowProblem<TypeTag>
{
// A few convenience aliases used throughout this class.
using ParentType = PorousMediumFlowProblem<TypeTag>;
using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
using Element = typename GridView::template Codim<0>::Entity;
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
static constexpr int dimWorld = GridView::dimensionworld;
public:
// This is the constructor of our problem class:
OnePTestProblem(std::shared_ptr<const GridGeometry> gridGeometry)
: ParentType(gridGeometry) {}
```
#### Boundary conditions
With the following function we define the __type of boundary conditions__ depending on the location.
Two types of boundary conditions can be specified: Dirichlet or Neumann boundary conditions. On
Dirichlet boundaries, the values of the primary variables need to be fixed. On a Neumann boundaries,
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 boundaryTypesAtPos(const GlobalPosition& globalPos) const
{
// we define a small epsilon value
Scalar eps = 1.0e-6;
// Initially, set Neumann boundary conditions for all equations
BoundaryTypes values;
values.setAllNeumann();
// On the top and bottom, use Dirichlet boundary conditions to prescribe pressures later.
const auto yMax = this->gridGeometry().bBoxMax()[dimWorld-1];
if (globalPos[dimWorld-1] < eps || globalPos[dimWorld-1] > yMax - eps)
values.setAllDirichlet();
return values;
}
```
The following function specifies the __values on Dirichlet boundaries__.
We need to define values for the primary variable (pressure), for which we
set a pressure of 1.1 bar and 1 bar at the bottom and top boundaries, respectively.
```cpp
PrimaryVariables dirichletAtPos(const GlobalPosition& globalPos) const
{
// instantiate a primary variables object
PrimaryVariables values;
// and assign a pressure value in [Pa] such that at the bottom boundary
// a pressure of 1.1 bar is set, and on the top boundary a pressure of 1 bar.
values[0] = 1.0e5*(1.1 - globalPos[dimWorld-1]*0.1);
return values;
}
```
#### Temperature distribution
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*/ }
}; // end class definition of OnePTestProblem
} // end namespace Dumux
```
</details>
## Parameter distributions (`spatialparams_1p.hh`)
This file contains the __spatial parameters class__ which defines the
distributions for the porous medium parameters permeability and porosity
over the computational grid
<details open>
<summary><b>Click to hide/show the file documentation</b> (or inspect the [source code](../spatialparams_1p.hh))</summary>
### Include files
In this example, we use a randomly generated and element-wise distributed
permeability field. For this, we use the random number generation facilities
provided by the C++ standard library.
```cpp
#include <random>
```
We include the spatial parameters class for single-phase models discretized
by finite volume schemes, from which the spatial parameters defined for this
example will inherit.
```cpp
#include <dumux/material/spatialparams/fv1p.hh>
```
### The spatial parameters class
In the `OnePTestSpatialParams` class, we define all functions needed to describe
the porous medium, e.g. porosity and permeability, for the single-phase problem.
We inherit from the `FVSpatialParamsOneP` class here, which is the base class
for spatial parameters in the context of single-phase porous medium flow
applications using finite volume discretization schemes.
```cpp
namespace Dumux {
template<class GridGeometry, class Scalar>
class OnePTestSpatialParams
: public FVSpatialParamsOneP<GridGeometry, Scalar,
OnePTestSpatialParams<GridGeometry, Scalar>>
{
// The following convenience aliases will be used throughout this class
using GridView = typename GridGeometry::GridView;
using FVElementGeometry = typename GridGeometry::LocalView;
using SubControlVolume = typename FVElementGeometry::SubControlVolume;
using Element = typename GridView::template Codim<0>::Entity;
using ParentType = FVSpatialParamsOneP<GridGeometry, Scalar,
OnePTestSpatialParams<GridGeometry, Scalar>>;
static constexpr int dimWorld = GridView::dimensionworld;
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;
```
#### Generation of the random permeability field
We generate a random permeability field upon construction of the spatial parameters class
using lognormal distributions. The mean values for the permeability inside and outside of a
low-permeable lens (given by the coorinates `lensLowerLeft_` and `lensUpperRight_`) are defined
in the variables `permeabilityLens_` and `permeability_`. The respective values are obtained
from the input file making use of the free function `getParam`. We use a standard deviarion
of 10% here and compute permeabily values for all elements of the computational grid.
```cpp
OnePTestSpatialParams(std::shared_ptr<const GridGeometry> gridGeometry)
: ParentType(gridGeometry), K_(gridGeometry->gridView().size(0), 0.0)
{
// The permeability of the domain and the lens are obtained from the `params.input` file.
permeability_ = getParam<Scalar>("SpatialParams.Permeability");
permeabilityLens_ = getParam<Scalar>("SpatialParams.PermeabilityLens");
// Furthermore, the position of the lens, which is defined by the position of the lower left and the upper right corners, are obtained from the input file.
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_`.
std::mt19937 rand(0);
std::lognormal_distribution<Scalar> K(std::log(permeability_), std::log(permeability_)*0.1);
std::lognormal_distribution<Scalar> KLens(std::log(permeabilityLens_), std::log(permeabilityLens_)*0.1);
// loop over all elements and compute a permeability value
for (const auto& element : elements(gridGeometry->gridView()))
{
const auto eIdx = gridGeometry->elementMapper().index(element);
const auto globalPos = element.geometry().center();
K_[eIdx] = isInLens_(globalPos) ? KLens(rand) : K(rand);
}
}
```
#### 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
instance of `ElementSolution`. Here, we use element-wise distributed permeabilities that were randomly generated in
the constructor and stored in the vector `K_`(see above).
```cpp
template<class ElementSolution>
const PermeabilityType& permeability(const Element& element,
const SubControlVolume& scv,
const ElementSolution& elemSol) const
{
return K_[scv.elementIndex()];
}
```
We set the porosity $`[-]`$ for the whole domain to a value of $`20 \%`$.
Note that in case you want to use solution-dependent porosities, you can
use the overload
`porosity(const Element&, const SubControlVolume&, const ElementSolution&)`
that is defined in the base class `FVSpatialParamsOneP`. Per default, this
fowards to the `porosityAtPos` function per default, which we overload here.
```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_; }
```
The remainder of this class contains a convenient function to determine if
a position is inside the lens and defines the data members.
<details><summary> Click to show private data members and member functions</summary>
```cpp
private:
// The following function returns true if a given position is inside the lens.
// We use an epsilon of 1.5e-7 here for floating point comparisons.
bool isInLens_(const GlobalPosition& globalPos) const
{
for (int i = 0; i < dimWorld; ++i)
{
if (globalPos[i] < lensLowerLeft_[i] + 1.5e-7
|| globalPos[i] > lensUpperRight_[i] - 1.5e-7)
return false;
}
return true;
}
GlobalPosition lensLowerLeft_, lensUpperRight_;
Scalar permeability_, permeabilityLens_;
std::vector<Scalar> K_;
}; // end class definition of OnePTestSpatialParams
} // end namespace Dumux
```
</details>
</details>
| [:arrow_left: Back to the main documentation](../README.md) | [:arrow_right: Continue with part 2](tracermodel.md) |
|---|---:|
# One-phase flow with random permeability distribution and a tracer model
# Single-phase flow and tracer transport
In this example, single-phase flow and tracer transport through a domain with a
heterogeneous permeability distribution is considered. A velocity distribution
is obtained from the solution of a stationary single-phase problem, and subsequently,
this velocity field is used for the simulation of the transport of tracer through the
domain.
__The main points illustrated in this example are__
* setting up and solving a stationary single-phase flow problem
* setting up and solving a tracer transport problem
* solving two problems sequentially and realizing the data transfer
* using a simple method to generate a random permeability field
__Table of contents__. This description is structured as follows:
[[_TOC_]]
## Problem set-up
This example contains a contaminant transported by a base groundwater flow in a randomly distributed permeability field. The figure below shows the simulation set-up. The permeability values range between 6.12e-15 and 1.5 e-7 $`m^2`$. A pressure gradient between the top and the bottom boundary leads to a groundwater flux from the bottom to the top. Neumann no-flow boundaries are assigned to the left and right boundary. Initially, there is a contaminant concentration at the bottom of the domain.
A domain with an extent of $`10 \, \mathrm{m} \times 10 \, \mathrm{m}`$ is considered,
in which a heterogeneous permeability distribution is generated randomly. The problem
set-up is shown in the figure below. In the stationary single-phase simulation,
a pressure difference between the bottom and the top boundaries is prescribed (see left figure),
resulting in a non-uniform velocity distribution due to the heterogeneous medium.
Neumann no-flow conditions are used on the lateral sides. On the basis of the resulting
velocity field, the transport of an initial tracer concentration distribution through
the domain is simulated (see right figure). Initially, non-zero tracer concentrations
are prescribed on a small strip close to the bottom boundary.
<figure>
<center>
<img src="img/1p_setup.png" alt="Single-phase setup" width="47%"/> <img src="img/xtracer.gif" alt="Tracer result" width="47%"/>
<figcaption> <b> Fig.1 </b> - Setup for the single-phase problem (left) and tracer mass fraction over time as computed with the tracer model (right).</figcaption>
</center>
</figure>
![](./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 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.
As mentioned above, two models are solved sequentially in this example. A single-phase
model (_1p model_) is used to solve for the stationary velocity distribution of a fluid phase
in the domain. The tracer transport is solved with the _tracer model_, which solves an advection-diffusion
equation for a tracer component, which is assumed not to affect the density and viscosity
of the fluid phase.
### 1p Model
The single phase model uses Darcy's law as the equation for the momentum conservation:
......@@ -21,29 +55,60 @@ with the darcy velocity $` \textbf v `$, the permeability $` \textbf K`$, the dy
Darcy's law is inserted into the mass balance equation:
```math
\phi \frac{\partial \varrho}{\partial t} + \text{div} \textbf v = 0
\phi \frac{\partial \varrho}{\partial t} + \text{div} \textbf v = 0,
```
where $`\phi`$ is the porosity.
The equation is discretized using a cell-centered 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).
where $`\phi`$ is the porosity. The primary variable used in this model is the pressure $`p`$.
### Tracer Model
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:
The tracer model solves the mass conservation equation of a tracer component $`\kappa`$,
in which both advective and diffusive transport mechanisms are considered:
```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.
```
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 (see [dumux/material/fluidmatrixinteractions/diffusivityconstanttortuosity.hh](https://git.iws.uni-stuttgart.de/dumux-repositories/dumux/-/blob/master/dumux/material/fluidmatrixinteractions/diffusivityconstanttortuosity.hh)):
Here, $`\textbf v`$ is a velocity field, which in this example is computed using the _1p model_ (see above). Moreover, $`X^\kappa`$ is the tracer mass fraction and $` D^\kappa_\text{pm} `$ is the
effective diffusivity. In this example, the effective diffusivity is a function of the diffusion
coefficient of the tracer component $`D^\kappa`$ and the porosity and tortuosity $`\tau`$ of the porous
medium (see [dumux/material/fluidmatrixinteractions/diffusivityconstanttortuosity.hh](https://git.iws.uni-stuttgart.de/dumux-repositories/dumux/-/blob/master/dumux/material/fluidmatrixinteractions/diffusivityconstanttortuosity.hh)):
```math
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](https://dumux.org/handbook).
The primary variable used in this model is the tracer mass fraction $`X^\kappa`$.