problem.md

---
title: Setting up a test problem / application and using the build system (CMake)
---

# Test Problems / Applications

## Test Problems / Applications
A test problem / application consists of:

* A property file (properties.hh)
* a problem file (often problem.hh)
* a spatial parameters file (often spatialparams.hh)
* an input file (often params.input)
* a main file (often main.cc)
* a build system file (CMakeLists.txt)

## Example for an immiscible two phase injection test case
Mass balance:

$\begin{equation}
\phi \frac{\partial \varrho_\alpha S_\alpha}{\partial t}
 -
 \text{div} \left\{ \boldsymbol{v}_\alpha \right\}
 -
 q_\alpha = 0
\end{equation}$

Momentum balance:

$\begin{equation}
\boldsymbol{v}_\alpha = \varrho_\alpha \frac{k_{r\alpha}}{\mu_\alpha} \mathbf{K} \left(\textbf{grad}\, p_\alpha - \varrho_{\alpha} \mathbf{g} \right)
\end{equation}$

# The properties file (properties.hh)

## The properties file
Lists all the properties and dependencies of the current problem. The injection test case inherits from the 2p model:

```cpp
namespace Dumux::Properties {
// define the TypeTag for this problem with a cell-centered two-point flux approximation spatial discretization.
// Create new type tags
namespace TTag {
struct Injection2p { using InheritsFrom = std::tuple<TwoP>; };
} // end namespace TTag
// Include all necessary properties within the namespace Dumux::Properties
} // end namespace Dumux::Properties
```

## The properties file
Often specifies the discretization method:

```cpp
namespace Dumux::Properties {
// define the TypeTag for this problem with a cell-centered two-point flux approximation spatial discretization.
// Create new type tags
namespace TTag {
struct Injection2p { using InheritsFrom = std::tuple<TwoP>; };
struct Injection2pCC { using InheritsFrom = std::tuple<Injection2p, CCTpfaModel>; };
} // end namespace TTag
} // end namespace Dumux::Properties
```

## The properties file
The grid type:

```cpp
// Set the grid type
template<class TypeTag>
struct Grid<TypeTag, TTag::Injection2p> { using type = Dune::YaspGrid<2>; };
```

## The properties file
The problem type:

```cpp
// Set the problem property
template<class TypeTag>
struct Problem<TypeTag, TTag::Injection2p> { using type = InjectionProblem2P<TypeTag>; };
```

## The properties file
The spatial parameters:

```cpp
// Set the spatial parameters
template<class TypeTag>
struct SpatialParams<TypeTag, TTag::Injection2p> {
private:
    using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
public:
    using type = InjectionSpatialParams<FVGridGeometry, Scalar>; };
```

## The properties file
The fluidsystem:

```cpp
// Set fluid configuration
template<class TypeTag>
struct FluidSystem<TypeTag, TTag::Injection2p> { using type = FluidSystems::H2ON2<GetPropType
                    <TypeTag, Properties::Scalar>, FluidSystems::H2ON2DefaultPolicy</*fastButSimplifiedRelations=*/true>>; };
```

## The properties file
The property file can also incorporate many more properties depending on the utilized model and test case.

# The problem file (problem.hh)

## The problem file
A problem in DuMu$^\mathsf{x}$ implements a specific model scenario:

```cpp
template<class TypeTag>
class InjectionProblem2P : public PorousMediumFlowProblem<TypeTag>
{
// Details of the model scenario (BoundaryConditions, InitialConditions, etc.)
}
```

## The problem file
Defines boundary conditions:

```cpp
BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
{
    BoundaryTypes bcTypes;
    // set the left of the domain (with the global position in "0 = x" direction as a Dirichlet boundary
    if (globalPos[0] < eps_)
        bcTypes.setAllDirichlet();
    // set all other as Neumann boundaries
    else
        bcTypes.setAllNeumann();
    return bcTypes;
}
```


## The problem file
Evaluates boundary conditions:

```cpp
PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
{
    return initialAtPos(globalPos);
}
```

## The problem file
```cpp
PrimaryVariables neumannAtPos(const GlobalPosition &globalPos) const
{
    // initialize values to zero, i.e. no-flow Neumann boundary conditions
    PrimaryVariables values(0.0);

    // using < boundary + eps_ or > boundary - eps_ is safer for floating point comparisons
    // than using <= or >= as it is robust with regard to imprecision introduced by rounding errors.
    if (time_ < injectionDuration_ && globalPos[1] < 15 + eps_ && globalPos[1] > 7 - eps_
        && globalPos[0] > 0.9*this->fvGridGeometry().bBoxMax()[0])
    {
        // inject nitrogen. Negative values mean injection
        // units kg/(s*m^2)
        values[Indices::conti0EqIdx + FluidSystem::N2Idx] = -1e-4;
        values[Indices::conti0EqIdx + FluidSystem::H2OIdx] = 0.0;
    }
    return values;
}
```

## The problem file
Defines initial conditions
```cpp
PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
{
    PrimaryVariables values(0.0);
    // get the water density at atmospheric conditions
    const Scalar densityW = FluidSystem::H2O::liquidDensity(temperature(), 1.0e5);

    // assume an initially hydrostatic liquid pressure profile
    const Scalar pw = 1.0e5 - densityW*this->gravity()[dimWorld-1]*(aquiferDepth_ - globalPos[dimWorld-1]);
    values[Indices::pressureIdx] = pw;
    values[Indices::saturationIdx] = 0.0;
    return values;
}
```


## The problem file
Defines source/sink terms

```cpp
NumEqVector sourceAtPos(const GlobalPosition &globalPos) const
{
    NumEqVector source(0.0)

    // The units must be according to either using mole or mass fractions.
    // (mole/(m^3*s) or kg/(m^3*s))
    // extract nitrogen at specific point. Positive values mean extraction.
    if (globalPos[0] < 5 + eps_ && globalPos[0] > 4 - eps_ && globalPos[1] < 5 + eps_ && globalPos[1] > 4 - eps_)
        source[Indices::conti0EqIdx + FluidSystem::N2Idx] = 1e-6;

    return source;
}
```

# The spatial parameters (spatialparams.hh)

## The spatial parameters
The spatialparams define spatial parameters of the porous material.

## The spatial parameters
Permeability:

```cpp
PermeabilityType permeabilityAtPos(const GlobalPosition& globalPos) const
{
    // here, either aquitard or aquifer permeability are returned, depending on the global position
    if (isInAquitard_(globalPos))
        return aquitardK_;
    return aquiferK_;
}
private:
// provides a convenient way distinguishing whether a given location is inside the aquitard
bool isInAquitard_(const GlobalPosition &globalPos) const
{
    // globalPos[dimWorld-1] is the y direction for 2D grids or the z direction for 3D grids
    return globalPos[dimWorld-1] > aquiferHeightFromBottom_ + eps_;
}
```

## The spatial parameters
Porosity:

```cpp
Scalar porosityAtPos(const GlobalPosition& globalPos) const
{
    // here, either aquitard or aquifer porosity are returned
    if (isInAquitard_(globalPos))
        return aquitardPorosity_;

    return aquiferPorosity_;
}
```

## The spatial parameters
Capillary pressure - saturation relationship:

More information in a later lecture on the materialsystem!

```cpp
const auto fluidMatrixInteractionAtPos(const GlobalPosition& globalPos) const
{
    if (isInAquitard_(globalPos))
        return makeFluidMatrixInteraction(aquitardPcKrSwCurve_);

    return makeFluidMatrixInteraction(aquiferPcKrSwCurve_);
}
```

## The spatial parameters
Temperature:

```cpp
Scalar temperatureAtPos(const GlobalPosition& globalPos) const
{
    // constant temperature of 20°C
    return 273.15 + 20.0;
}
```

# The input file (*.input)

## The input file
DUNE INI syntax:

```cpp
[Grid]
LowerLeft = 0 0
UpperRight = 60 40
Cells = 24 16

[Problem]
Name = test
```
Input files are specified as arguments to the executable

`./myexecutable params.input`

## The input file
If no input file is given it defaults to the file `params.input` or `myexecutablename.input`

Parameters can be overwritten through the command line like via:

`./executable –Problem.Name myNewName`

# The main source file (*.cc)

## The main source file
* Each problem has a specific main file (`test_name.cc` or `main.cc`) which sets up the program structure and calls assembler and solvers to assemble and solve the PDEs.
* Depending on the complexity of the problem the main file can be either set up to solve a linear problem, a non-linear problem and stationary as well as instationary problems.
* The main file usually includes the problem, the solvers, the assembler, the VTK output module and the gridmanager.

## The main source file
Common structure for most main files:

```cpp
// define the type tag for this problem
using TypeTag = Properties::TTag::OnePIncompressible;

// maybe initialize MPI and/or multithreading backend
const auto& mpiHelper = Dune::MPIHelper::instance();

// print dumux start message
if (mpiHelper.rank() == 0)
    DumuxMessage::print(/*firstCall=*/true);

// parse command line arguments and input file
Parameters::init(argc, argv);
```

## The main source file
```cpp
// try to create a grid (from the given grid file or the input file)
GridManager<GetPropType<TypeTag, Properties::Grid>> gridManager;
gridManager.init();

// we compute on the leaf grid view
const auto& leafGridView = gridManager.grid().leafGridView();

// create the finite volume grid geometry
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
auto gridGeometry = std::make_shared<GridGeometry>(leafGridView);

// the problem (initial and boundary conditions)
using Problem = GetPropType<TypeTag, Properties::Problem>;
auto problem = std::make_shared<Problem>(gridGeometry);
```

## The main source file
```cpp
// the solution vector
using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
SolutionVector x(gridGeometry->numDofs());

// the grid variables
using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
auto gridVariables = std::make_shared<GridVariables>(problem, gridGeometry);
gridVariables->init(x);

// initialize the vtk output module
VtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables, x, problem->name());
using VelocityOutput = GetPropType<TypeTag, Properties::VelocityOutput>;
vtkWriter.addVelocityOutput(std::make_shared<VelocityOutput>(*gridVariables));
using IOFields = GetPropType<TypeTag, Properties::IOFields>;
IOFields::initOutputModule(vtkWriter); //!< Add model specific output fields
```

## The main source file

Differences for various problem cases:

* Stationary linear problem
* Stationary non-linear problem
* Instationary non-linear problem

## The main source file

A stationary linear problem:

```cpp
// the assembler for stationary problems
using Assembler = FVAssembler<TypeTag, DiffMethod::numeric>;
auto assembler = std::make_shared<Assembler>(problem, gridGeometry, gridVariables);

// the linear solver
using LinearSolver = ILUBiCGSTABIstlSolver<LinearSolverTraits<GridGeometry>,
                                           LinearAlgebraTraitsFromAssembler<Assembler>>;
auto linearSolver = std::make_shared<LinearSolver>(gridGeometry->gridView(), gridGeometry->dofMapper());

// the discretization matrices for stationary linear problems
using JacobianMatrix = GetPropType<TypeTag, Properties::JacobianMatrix>;
auto A = std::make_shared<JacobianMatrix>();
auto r = std::make_shared<SolutionVector>();
assembler->setLinearSystem(A, r);
assembler->assembleJacobianAndResidual(x);
```

## The main source file
```cpp
// we solve Ax = -r to save update and copy
(*r) *= -1.0;
linearSolver->solve(*A, x, *r);
// the grid variables need to be up to date for subsequent output
gridVariables->update(x);
```

## The main source file
A stationary non-linear problem:

```cpp
using Assembler = FVAssembler<TypeTag, DiffMethod::numeric>;
auto assembler = std::make_shared<Assembler>(problem, gridGeometry, gridVariables);
// the linear solver
using LinearSolver = ILUBiCGSTABIstlSolver<LinearSolverTraits<GridGeometry>,
                                           LinearAlgebraTraitsFromAssembler<Assembler>>;
auto linearSolver = std::make_shared<LinearSolver>(gridGeometry->gridView(), gridGeometry->dofMapper());
// the non-linear solver
using NewtonSolver = Dumux::NewtonSolver<Assembler, LinearSolver>;
NewtonSolver nonLinearSolver(assembler, linearSolver);
// linearize & solve
nonLinearSolver.solve(x);
```

## The main source file
An instationary non-linear problem:

```cpp
// instantiate time loop
auto timeLoop = std::make_shared<TimeLoop<Scalar>(0.0, dt, tEnd);
timeLoop->setMaxTimeStepSize(maxDt);
using Assembler = FVAssembler<TypeTag, DiffMethod::numeric>;
auto assembler = std::make_shared<Assembler>(problem, gridGeometry, gridVariables, timeLoop);

// the linear solver
using LinearSolver = ILUBiCGSTABIstlSolver<LinearSolverTraits<GridGeometry>,
                                           LinearAlgebraTraitsFromAssembler<Assembler>>;
auto linearSolver = std::make_shared<LinearSolver>(gridGeometry->gridView(), gridGeometry->dofMapper());

// the non-linear solver
using NewtonSolver = Dumux::NewtonSolver<Assembler, LinearSolver>;
NewtonSolver nonLinearSolver(assembler, linearSolver);
```

## The main source file
```cpp
// time loop
timeLoop->start(); do
{
    // Calculate solution within each time step
} while (!timeLoop->finished());
```

## The main source file
The timeloop for the instationary non-linear problem:

```cpp
// time loop
timeLoop->start(); do
{
    // set previous solution for storage evaluations
    assembler->setPreviousSolution(xOld);

    // linearize & solve
    nonLinearSolver.solve(x, *timeLoop);

    // make the new solution the old solution
    xOld = x;
    gridVariables->advanceTimeStep();
```

## The main source file
```cpp
    // advance to the time loop to the next step
    timeLoop->advanceTimeStep();
    // report statistics of this time step
    timeLoop->reportTimeStep();
    // set new dt as suggested by the newton solver
    timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()));
} while (!timeLoop->finished());
timeLoop->finalize(leafGridView.comm());
```

# Build system (CMakeLists.txt)

## Build system - what is CMake?

* Open source build system tool developed by KITware.
* Offers a one-tool-solution to all building tasks, like configuring, building, linking, testing and packaging.
* Is a build system generator: It supports a set of backends called generators.
* Is portable and supports cross-platform compilation.
* Is controlled by ONE rather simple language.
* Every directory in a project contains a file called `CMakeLists.txt`, which is written in the CMake language. You can think of these as a distributed configure script. Upon configure, the top-level `CMakeLists.txt` is executed.

## Build system - configuring

* Configure build time compiler parameters / linking information.
* Create „targets“ that can be build to create executables.

## Build system - configuring

* Build with the script `dune-common/bin/dunecontrol <options>` which takes care of all dependencies and modular dune structure.
* Option `all`: build all libraries and executables.
* Option `--opts=<optionfile.opts>` specify e.g. compiler flags,  for DuMu$^\mathsf{x}$ the `optionfile` is `dumux/cmake.opts`.
* Option `--build-dir=<build directory>` specify path for out-of-source build.
* Default: every module contains its own build directory `build-cmake/`.
* You have to reconfigure (possibly deleting all build directories first) whenever a dependency changes or a Dune library is updated.

`./dune-common/bin/dunecontrol --opts=dumux/cmake.opts all`

## Build system - important basic commands
* Use `add_subdirectory` for recursively adding subdirectories.
* The subdirectory has to contain a `CMakeLists.txt` file (can be empty).
* Executables are added via `add_executable(<name> source1 [source2 ...])`.
* Tests are added via `dune_add_test(...)` which also add a test executable to the test suite.
* Symlinks can be added via `dune_symlink_to_source_files(FILES file1 [file2 ...])`.

## Build system
Simplest incorporation of a test by defining name, source file and command line arguments:
```cmake
dune_add_test(NAME test_2p_incompressible_box
SOURCES test_2p_fv.cc
CMD_ARGS test_2p.input)
```

## Build system
Add extra compile definitions and commands:
```cmake
dune_add_test(NAME test_2p_incompressible_box
SOURCES test_2p_fv.cc
COMPILE_DEFINITIONS TYPETAG=TwoPIncompressibleBox
COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
CMD_ARGS --script fuzzy
         --files ${CMAKE_SOURCE_DIR}/test/references/lensbox-reference.vtu
                 ${CMAKE_CURRENT_BINARY_DIR}/2p_box-00007.vtu
         --command "${CMAKE_CURRENT_BINARY_DIR}/test_2p_incompressible_box
                   test_2p.input -Problem.Name 2p_box")
```

## Build system
Extra: Create linked parameter file in `build-cmake` folder, separately add executable and set compile defintions for an executable:
```cmake
dune_symlink_to_source_files(FILES "params.input")
# using tpfa
add_executable(test_2p_incompressible_tpfa EXCLUDE_FROM_ALL main.cc)
target_compile_definitions(test_2p_incompressible_tpfa PUBLIC TYPETAG=TwoPIncompressibleTpfa)
```

## Build system
Instead of the box discretization use the tpfa discretization:
```cmake
dumux_add_test(NAME test_2p_incompressible_tpfa
              TARGET test_2p_incompressible_tpfa
              LABELS porousmediumflow 2p
              COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
              CMD_ARGS --script fuzzy
                       --files ${CMAKE_SOURCE_DIR}/test/references/test_2p_incompressible_cc-reference.vtu
                               ${CMAKE_CURRENT_BINARY_DIR}/test_2p_incompressible_tpfa-00007.vtu
                       --command "${CMAKE_CURRENT_BINARY_DIR}/test_2p_incompressible_tpfa params.input -Problem.Name test_2p_incompressible_tpfa")
```

## Build system
Important basic commands:

* See also Dune build system documentation on [https://www.dune-project.org/sphinx/core/](https://www.dune-project.org/sphinx/core/) for a comprehensive CMake online documentation.

# Parallelism

## Distributed memory parallelism with MPI

* MPI stands for Message Passing Interface
* Main idea is the concept of domain decomposition
* Each local subdomain is solved on an individual process(rank)
* MPI manages the communication between the ranks
* Most solvers in DuMu$^\mathsf{x}$ are capable of parallel solving

## Distributed memory parallelism with MPI

Run with:
```cpp
   mpirun -np [n_cores] [executable_name]
```
Handling results:

* Each rank creates its own `*.vtu`/`*.vtp` file.
* These are combined into `*.pvtu`/`*.pvtp` files for each time step.
* A normal `*.pvd` file is created from the `*.pvtu`/`*.pvtp` files.

## Shared-memory parallelism and multi-threaded applications

* Dumux can exploit parallelism with the shared memory model
* Used in the Dumux::FVAssembler by default to assemble the residual and stiffness matrix
* Is enabled when a multi-threading backend is found
* Backend is selected by `CMAKE` during configuration and stored in `DUMUX_MULTITHREADING_BACKEND`
* Possible examples are `OpenMP`, `TBB`, C++ parallel algorithms and `Kokkos`
* Can be turned of in `params.input` with:
```cpp
[Assembly]
Multithreading = false
```
* Important for working on clusters: Number of threads can also be restricted via manipulating the environment variable `DUMUX_NUM_THREADS=2 ./executable`

# Exercises

## Exercises:
Exercise about setting boundary conditions, the problem file etc:

* Go to [https://git.iws.uni-stuttgart.de/dumux-repositories/dumux-course/tree/master/exercises/exercise-basic](https://git.iws.uni-stuttgart.de/dumux-repositories/dumux-course/tree/master/exercises/exercise-basic) and check out the README

Exercise for the main-files:

* Go to [https://git.iws.uni-stuttgart.de/dumux-repositories/dumux-course/tree/master/exercises/exercise-mainfile](https://git.iws.uni-stuttgart.de/dumux-repositories/dumux-course/tree/master/exercises/exercise-basic) and check out the README