Guided Doxygen Documentation ¶

Contents

Guided Doxygen Documentation
- Introduction
  - Programming Philosophy
  - Programming Guide
- DORiE Routine

Introduction ¶

This documentation includes the documentation created by Doxygen, and enhances it by adding further overview features and explanations not suitable for inclusion in the source code.

DORiE consists of several thousand lines of C++ source code, not accounting for the DUNE functions and classes included in the program. There are two approaches on getting to know the code, first by following the programs course fairly linearly, and second by studying class functions and interactions via the modules overview. - Don’t worry! The modularity of the code enables you to work with portions of it while neglecting other parts.

Programming Philosophy ¶

DORiE makes high use of object oriented programming by constructing different classes and letting them interact quite independently.

In general, the constructor of every class reads the parameters important to the same class and immediately assembles all variables and functions needed for further use. Class constructors usually do not receive parameters explicitly, but read them from the parameter file independently. This improves readability of the code but can cause conflicts in parameter declarations. The dorie.parscraper makes sure that all parameter queries are consistent before DORiE is compiled.

Query functions are handled in a top-down approach, where requests for and the sending of data are initiated by the higher-order or owner class. Ownership itself, however, is intended to be as flat as possible, with the main program functions assembling most of the classes and handing them over to other classes with which they interact.

Programming Guide ¶

None so far.

DORiE Routine ¶

Main Function ¶

The Main Routine of DORiE creates a Parameter File Parser object, to read out the parameter file. It assembles a Dune::Grid object by calling the appropriate functions for a structured rectangular or an unstructured simplex grid. The latter is created from a Gmsh .msh Mesh data file.

The grid is handed over to the appropriate Solver functions, which build the Dune Grid Function Space for a rectangular or simplex grid, respectively. The Solver helper functions then call the actual code-solver for computing the solution.

Typedefs

using Sim = Dune::Dorie::Simulation<Traits>¶

using Simplex = Dune::Dorie::BaseTraits<Dune::UGGrid, Dune::GeometryType::BasicType::simplex, dim, order, true, false>¶

using SimplexAdaptive = Dune::Dorie::BaseTraits<Dune::UGGrid, Dune::GeometryType::BasicType::simplex, dim, order, true, true>¶

using Cube = Dune::Dorie::BaseTraits<Dune::Dorie::YaspGrid, Dune::GeometryType::BasicType::cube, dim, order, true, false>¶

using CubeAdaptive = Dune::Dorie::BaseTraits<Dune::UGGrid, Dune::GeometryType::BasicType::cube, dim, order, true, true>¶

Functions

int main(int argc, char **argv)¶

Main Program Function: Initialize parameter file reader, build grid and call Richards Solver.

As simplex and rectangular grids demand different FiniteElementMaps, the program calls different functions for these tasks. The objects and types are then passed to the generic solver function.

The UG Grid and FEM templates need the dimension as

constexpr at compile time, so we need the tedious dim and FEorder queries.

See: RichardsSolver Main functions for assembling operators and solving the problem
See: RichardsSolverSimplex Helper for building simplex Grid Function Spaces
See: RichardsSolverRectangular Helper for building rectangular Grid Function Spaces
See: Dune::Dorie::Traits Variable type definitions Parameter file is specified Output directory is writable Finite Element Order is supported Grid type is supported Dimensions are supported

namespace Dune¶

namespace Dorie¶

Typedefs

using Dune::Dorie::YaspGrid = typedef Dune::YaspGrid<dim>: Resolve the second (default) template parameter of YaspGrid.

template <class GridType> std::shared_ptr<GridType> Dune::Dorie::build_grid_gmsh(const Dune::ParameterTree &inifile, const Dune::MPIHelper &helper)¶

Read simplex Gmsh file in 2D or 3D, create a grid, and return a shared pointer to the grid.

Refine the grid globally and call default load balancing for multi-processing support.

Return

Shared pointer to grid

Template Parameters

GridType: Type of the grid. Must be supported by Dune::GridFactory

Parameters

inifile: Parameter file parser

template <class GridType> std::shared_ptr<GridType> Dune::Dorie::build_grid_cube(const Dune::ParameterTree &inifile, const Dune::MPIHelper &helper)¶

Read parameter file specifications, create a grid, and return a shared pointer to the grid.

Refine the grid globally and call default load balancing for multi-processing support.

Return

Shared pointer to grid

Template Parameters

GridType: Type of the grid. Must be supported by Dune::StructuredGridFactory

Parameters

inifile: Parameter file parser

Simulation ¶

The Simulation template receives the assembled Grid Function Space and initializes the PDELab Backends to compute the solution.

template <typename Traits> class Simulation¶

Basic Simulation class providing objects and functions for computing the solution.

Subclassed by Dune::Dorie::KnoFuInterface< Traits >

Richards Equation Classes ¶

To compute the solution, the various parts of the Richards Equation and the chosen parameterization have to be represented. This is achieved by building objects which provide appropriate queries for Boundary Conditions, Parameters and Sink/Source terms, depending on both time and space. There are three classes with standardized function queries for these tasks:

template <typename T, class Interpolator> class ParameterizationBase¶

Base class for all parameterizations.

Subclassed by Dune::Dorie::MualemVanGenuchten< T, Interpolator >

Public Functions

void read_parameters()¶: Iterate through parameter array and read each one from the H5 file.

RF headRefToMiller(const RF &h, const Domain &pos) const¶

Transformation: Head: Reference to Miller similar (if activated)

\[ h^* = h \exp\left( - (v - \text{Var}(v))\right) \]

with the raw random field value \(v\) and its variance.

Return

Transformed value

Parameters

h: Matric head
pos: Position in global coordinates

RF headMillerToRef(const RF &h, const Domain &pos) const¶

Transformation: Head: Miller similar to Reference (if activated)

See: headRefToMiller

RF condRefToMiller(const RF &K, const Domain &pos) const¶

Transformation: Conductivity: Reference to Miller similar (if activated)

\[ K^* = K \exp\left( 2 (v - \text{Var}(v))\right) \]

with the raw random field value \(v\) and its variance.

Return

Transformed value

Parameters

h: Matric head
pos: Position in global coordinates

RF condMillerToRef(const RF &k, const Domain &pos) const¶

Transformation: Conductivity: Miller similar to Reference (if activated)

See: condRefToMiller

const Vector &gravity() const¶: Return gravity vector.

template <typename T, class Interpolator> class MualemVanGenuchten : public Dune::Dorie::ParameterizationBase<T, Interpolator>¶

Mualem-van Genuchten parameterization class.

Public Functions

RF saturation(const RF &h, const Domain &pos) const¶

Return water saturation for given matric head at current position.

\[ \Theta(h) = \left( 1 + \left( \alpha h \right)^n \right)^{-m} \]

with \(m = 1 - 1/n\)

Return

Soil water saturation

Parameters

h: Matric head

RF waterContent(const RF &h, const Domain &pos) const¶

Return volumetric water content for given matric head at current position.

\[ \theta(h) = \Theta(h) (\theta_s - \theta_r) + \theta_r \]

Return

Soil Water Content

Parameters

h: Matric head

RF K(const Domain &pos) const¶

Return saturated conductivity at current position.

\[ K = K_0 \]

Return: Saturated conductivity

RF condFactor(const RF &sat, const Domain &pos) const¶

Return conductivity factor for given water saturation at current position.

Saturated conductivity K is multiplied with factor condFactor < 1 to receive saturation-dependent conductivity.

\[ K_\text{fac} = \Theta^\tau \left( \Theta - \left( 1 - \Theta^{1/m} \right)^m \right)^2 \]

with \(m = 1 - 1/n\)

Return

Conductivity Factor

Parameters

sat: Soil water saturation

void read_parameters(): Iterate through parameter array and read each one from the H5 file.

RF headRefToMiller(const RF &h, const Domain &pos) const

Transformation: Head: Reference to Miller similar (if activated)

\[ h^* = h \exp\left( - (v - \text{Var}(v))\right) \]

with the raw random field value \(v\) and its variance.

Return

Transformed value

Parameters

h: Matric head
pos: Position in global coordinates

RF headMillerToRef(const RF &h, const Domain &pos) const

Transformation: Head: Miller similar to Reference (if activated)

See: headRefToMiller

RF condRefToMiller(const RF &K, const Domain &pos) const

Transformation: Conductivity: Reference to Miller similar (if activated)

\[ K^* = K \exp\left( 2 (v - \text{Var}(v))\right) \]

with the raw random field value \(v\) and its variance.

Return

Transformed value

Parameters

h: Matric head
pos: Position in global coordinates

RF condMillerToRef(const RF &k, const Domain &pos) const

Transformation: Conductivity: Miller similar to Reference (if activated)

See: condRefToMiller

const Vector &gravity() const: Return gravity vector.

template <typename Traits> class FlowBoundary¶

Boundary type and condition value queries.

This class containts functions that return the type of boundary conditions and the values of the boundary condition parameters. Both types of queries can be time dependent.

Public Functions

FlowBoundary(const Dune::ParameterTree &config)¶

Read BC file type and create data handling object.

See

Dune::Dorie::BCReadoutInterface

Exceptions

Dune::IOError: BC File Type not supported

RF getNextTimeStamp(const RF &time) const¶

Return next time at which any boundary condition changes. Return -1.0 by default.

Parameters

time: Current time

BoundaryCondition::Type bc(const Intersection &is, const ID &x, const RF &time) const¶

Return Boundary Condition Type at certain position and time.

Return

Boundary condition type enumerator

Parameters

is: Intersection enitity
x: Position in local entity coordinates
time: Time value
xGlobal: Global coordinate

RF g(const Intersection &is, const ID &x, const RF &time) const¶

Dirichlet boundary condition at certain position and time.

Return

Value of matric head

Parameters

is: Intersection enitity
x: Position in local entity coordinates
time: Time value
xglobal: Global coordinate

RF j(const Intersection &is, const ID &x, const RF &time) const¶

Neumann boundary condition at certain position and time.

Return

Value of flux

Parameters

is: Intersection enitity
x: Position in local entity coordinates
time: Time value
xglobal: Global coordinate

template <typename Traits, typename Parameter> class FlowSource¶

Sink and source terms for flow equation.

This class implements sinks (q<0) and sources (q>0) inside of the grid.

Public Functions

RF q(const Element &elem, const Domain &x, RF time) const¶

Return value of the sink/source at certain position and time.

Return

Sink/source term value

Parameters

elem: Element enitity
x: Position in local entity coordinates
time: Time value
xglobal: Global coordinates

Richards Equation Local Operators ¶

The local operators bind to the PDELab backend and compute the local residual along quadrature points on the grid. They represent the numeric implementation of the PDE’s strong formulation. DORiE uses the mass conserving Discontinuous Galerkin (DG) discretization scheme.

template <typename Traits, typename Parameter, typename Boundary, typename SourceTerm, typename FEM, bool adjoint> class RichardsDGSpatialOperator : public Dune::PDELab::NumericalJacobianVolume<RichardsDGSpatialOperator<Traits, Parameter, Boundary, SourceTerm, FEM, adjoint>>, public Dune::PDELab::NumericalJacobianSkeleton<RichardsDGSpatialOperator<Traits, Parameter, Boundary, SourceTerm, FEM, adjoint>>, public Dune::PDELab::NumericalJacobianBoundary<RichardsDGSpatialOperator<Traits, Parameter, Boundary, SourceTerm, FEM, adjoint>>, public FullSkeletonPattern, public FullVolumePattern, public LocalOperatorDefaultFlags, public Dune::PDELab::InstationaryLocalOperatorDefaultMethods<Traits::RangeField>¶

A local operator for solving the Richards equation using upwinding.

Written by Ole Klein, IWR, Heidelberg and modified by the DORiE developers, IUP, Heidelberg

In theory it is possible that one and the same local operator is called first with a finite element of one type and later with a finite element of another type. Since finite elements of different type will usually produce different results for the same local coordinate they cannot share a cache. Here we use a vector of caches to allow for different orders of the shape functions, which should be enough to support p-adaptivity. (Another likely candidate would be differing geometry types, i.e. hybrid meshes.)

Public Functions

RichardsDGSpatialOperator(const Dune::ParameterTree &config, const Parameter &param_, const GridView view_, Boundary &boundary_, const SourceTerm &sourceTerm_, RichardsDGMethod::Type method_ = RichardsDGMethod::SIPG, RichardsDGUpwinding::Type upwinding_ = RichardsDGUpwinding::none, RichardsDGWeights::Type weights_ = RichardsDGWeights::weightsOn, int intorderadd_ = 2, int quadrature_factor_ = 2)¶

Create operator.

See

RichardsDGMethod::Type

See

RichardsDGUpwinding::Type

See

RichardsDGWeights::Type

Parameters

param_: Object providing equation parameters
view_: Grid view
boundary_: Object prividing boundary conditions
sourceTerm_: Object providing source term information
method_: DG method type

Parameters

upwinding_: Upwinding type

Parameters

weights_: DG weigths switch

Parameters

intorderadd_: Addend to integration order
quadrature_factor_: Factor to FEM integration order

template <typename EG, typename LFSU, typename X, typename LFSV, typename R> void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const¶

Volume integral depending on test and ansatz functions.

Richards Equation - Strong formulation of volume integral

\[ \int_{\Omega_e} K (\Theta (h_m^*)) \cdot (\nabla u - \mathbf{\hat{g}}) \cdot \nabla v \; dx \]

Miller Similarity Transformations:

\[\begin{split}\begin{eqnarray*} h_m \; &\xrightarrow{\text{MS}^{-1}} \; & h_m^* \\ K^*(\Theta(h_m^*)) \; &\xrightarrow{\text{MS}} \; & K(\Theta(h_m^*)) \end{eqnarray*}\end{split}\]

template <typename IG, typename LFSU, typename X, typename LFSV, typename R> void alpha_skeleton(const IG &ig, const LFSU &lfsu_s, const X &x_s, const LFSV &lfsv_s, const LFSU &lfsu_n, const X &x_n, const LFSV &lfsv_n, R &r_s, R &r_n) const¶

Skeleton integral depending on test and ansatz functions.

Each face is only visited ONCE!

Weights for symmetry term:

\[ \omega_x = \frac{ K_{0,y} }{ K_{0,x} + K_{0,y}} \; , \quad x,y \in \lbrace s,n \rbrace , \; x \neq y \]

Solution jump over face:

\[ \Delta u = u_s - u_n \]

Penalty factor:

\[\begin{split}\begin{eqnarray*} p &= \mathsf{p} \cdot k (k+d-1) / h_f \\ h_f &= \frac{\min \lbrace \Omega_s,\Omega_n \rbrace }{\Omega_i}\\ d \; &: \text{spatial dimension}\\ k \; &: \text{polynomial order}\\ \Omega_s,\Omega_n \in \mathbb{R}^d \; &: \text{element volumes}\\ \Omega_i \in \mathbb{R}^{d-1} \; &: \text{interface surface} \end{eqnarray*}\end{split}\]

template <typename IG, typename LFSU, typename X, typename LFSV, typename R> void alpha_boundary(const IG &ig, const LFSU &lfsu_s, const X &x_s, const LFSV &lfsv_s, R &r_s) const¶

Boundary integral depending on test and ansatz functions.

template <typename EG, typename LFSV, typename R> void lambda_volume(const EG &eg, const LFSV &lfsv, R &r) const¶: Volume integral depending only on test functions.

void setTime(RF t)¶: set operator time

void preStep(RF t, RF dt, int stages)¶: set operator time and erase BC cache map

template <typename Traits, typename Parameter, typename FEM, bool adjoint> class RichardsDGTemporalOperator : public Dune::PDELab::NumericalJacobianVolume<RichardsDGTemporalOperator<Traits, Parameter, FEM, adjoint>>, public Dune::PDELab::NumericalJacobianApplyVolume<RichardsDGTemporalOperator<Traits, Parameter, FEM, adjoint>>, public FullSkeletonPattern, public FullVolumePattern, public LocalOperatorDefaultFlags, public Dune::PDELab::InstationaryLocalOperatorDefaultMethods<Traits::RangeField>¶

DG local operator for the temporal derivative of the Richards equation.

Public Functions

RichardsDGTemporalOperator(const Dune::ParameterTree &config, const Parameter &param_, const typename Traits::GridView view, int intorderadd_ = 0, int quadrature_factor_ = 2)¶: Constructor.

template <typename EG, typename LFSU, typename X, typename LFSV, typename R> void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const¶: Volume integral depending on test and ansatz functions.

void setTime(RF t)¶: set operator time

Time Step Controller ¶

The time step interval is controlled by an own class. It is adjusted according to the performance of the non-linear Newton Solver and possible changes in the boundary conditions.

template <typename R, typename TSC> class CalculationController¶

Class to control time and maximum number of iterations for the Newton solver.

Public Functions

CalculationController(const ParameterTree &config, const TSC &tsc_, const Dune::MPIHelper &helper)¶

Read relevant time values, check for consistency and adapt time step to BC change, if necessary.

Parameters

tsc_: Class with public function getNextTimeStamp(time)

R getTime() const¶: Return current time.

R getDT() const¶: Return timestep for next calculation.

int getIterations() const¶: Return maximum number of Newton iterations allowed for next calculation.

bool doStep() const¶: Return boolean whether new time step shall be computed.

void set_t_end(const R _t_end)¶: Set a new finish time for the computation.

bool validate(const bool exc)¶

Validation of solution and time step calculation.

If the Newton solver throws an exception (because the solution did not converge in the allowed steps or the timestep is too large), exc=true is passed to this function and the timestep will be reduced. Returning false causes the Problem solver to recalculate the timestep with new parameters provided by the controller. The next time step is also adjusted to the next time the BC changes

Return

Boolean if solution is valid

Parameters

exc: Boolean if Newton Class threw an exception

Exceptions

Dune::Exception: Newton solver does not converge at dt>=dtmin and it<=itmax