DuMu^x applications

Example application

Simulation goal

Gas injection / immiscible two phase flow

Mass balance equations for two fluid phases: \[ \begin{aligned} \frac{\partial \left(\phi \varrho_\alpha S_\alpha \right)}{\partial t} + \nabla \cdot \left(\varrho_\alpha \boldsymbol{v}_\alpha \right) - q_\alpha = 0, \quad \alpha \in \lbrace w, n \rbrace. \end{aligned} \]

Momentum balance equations (multiphase-phase Darcy’s law): \[ \begin{aligned} \boldsymbol{v}_\alpha = -\frac{k_{r\alpha}}{\mu_\alpha} \mathbf{K} \left(\nabla\, p_\alpha - \varrho_{\alpha} \mathbf{g} \right), \quad \alpha \in \lbrace w, n \rbrace. \end{aligned} \]

Gas injection / immiscible two phase flow

\[ \begin{aligned} \frac{\partial \left(\phi \varrho_\alpha S_\alpha \right)}{\partial t} - \nabla \cdot \left( \varrho_\alpha \frac{k_{r\alpha}}{\mu_\alpha} \mathbf{K} \left(\nabla\, p_\alpha - \varrho_{\alpha} \mathbf{g} \right) \right) - q_\alpha = 0 \end{aligned} \]

\(p_w\), \(p_n\): wetting and non-wetting fluid phase pressure
\(\varrho_\alpha\), \(\mu_\alpha\): fluid phase density and dynamic viscosity
\(\phi\), \(\mathbf{K}\), \(k_{r\alpha}\): porosity, intrinsic permeability, relative permeability
\(S_w\), \(S_n\): wetting and non-wetting saturation
\(\mathbf{g}\), \(q_\alpha\): gravitational potential, source term

Gas injection / immiscible two phase flow

Constitutive relations: \(p_n := p_w + p_c\), \(p_c := p_c(S_w)\), \(k_{r\alpha}\) = \(k_{r\alpha}(S_w)\)
Physical constraint: \(S_w + S_n = 1\)
Primary variables: \(p_w\), \(S_n\) (wetting phase pressure, non-wetting phase saturation)

A DuMu^x application

The source code for a typical DuMu^x application consists of:

Property specializations (e.g. properties.hh)
Problem class definition (e.g. problem.hh) and often a spatial parameter class definition (e.g. spatialparams.hh)
A parameter configuration file (e.g. params.input)
The main program (e.g. main.cc)
CMakeLists.txt (CMake build system configuration)

Property specializations

Properties

“Compile-time configuration”

“Tell DuMu^x what type of classes to use”

Details on properties in Part II: Property system

*properties*.hh

Create tag, choose model and discretization scheme:

namespace Dumux::Properties::TTag {

// the application type tag (choose any name)
struct Injection2pCC { using InheritsFrom = std::tuple<TwoP, CCTpfaModel>; };

} // end namespace Dumux::Properties::TTag

Specialize the Grid property for tag TTag::Injection2pCC to set the grid type:

namespace Dumux::Properties {

template<class TypeTag>
struct Grid<TypeTag, TTag::Injection2pCC>
{ using type = Dune::YaspGrid<2>; }; // Yet Another Structured Parallel Grid

} // end namespace Dumux::Properties

Specialize the Problem property to set the problem type:

namespace Dumux::Properties {

template<class TypeTag>
struct Problem<TypeTag, TTag::Injection2pCC>
{ using type = InjectionProblem2P<TypeTag>; };

} // end namespace Dumux::Properties

The problem class InjectionProblem2P is discussed in one of the following sections.

Specialize the SpatialParams property to set the spatial parameter type:

namespace Dumux::Properties {

template<class TypeTag>
struct SpatialParams<TypeTag, TTag::Injection2pCC>
{
    using GridGeometry
        = GetPropType<TypeTag, Properties::GridGeometry>;
    using Scalar
        = GetPropType<TypeTag, Properties::Scalar>;
    using type = InjectionSpatialParams<GridGeometry, Scalar>;
};

} // end namespace Dumux::Properties

Specialize the FluidSystem property to set the fluid system type:

namespace Dumux::Properties {

template<class TypeTag>
struct FluidSystem<TypeTag, TTag::Injection2pCC>
{
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
    using Policy = FluidSystems::H2ON2DefaultPolicy<
        /*fastButSimplifiedRelations=*/true>;
    using type = FluidSystems::H2ON2<Scalar, Policy>;
};

} // end namespace Dumux::Properties

Problem definition

Problem

A “problem” class in DuMu^x implements a specific scenario:

*problem*.hh

template<class TypeTag>
class InjectionProblem2P : public PorousMediumFlowProblem<TypeTag>
{
    // - boundary conditions
    // - initial conditions
    // - source terms
    // - scenario name (for output)
}

Inherit from base class PorousMediumFlowProblem, only override scenario-specific functions (static polymorphism).

The boundary condition types:

BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
{
    BoundaryTypes bcTypes;
    // Use Dirichlet boundary conditions on the left
    if (globalPos[0] < eps_)
        bcTypes.setAllDirichlet();
    // Use Neumann boundary conditions on the rest of the boundary
    else
        bcTypes.setAllNeumann();
    return bcTypes;
}

Dirichlet boundary condition values (only called on Dirichlet boundaries)

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

PrimaryVariables is an array of primary variables (here, the size of the array is 2).

More general Dirichlet boundary functions:

PrimaryVariables dirichlet(const Element &element, 
                           const SubControlVolumeFace &scvf) const
{
    ...
}

PrimaryVariables dirichlet(const Element &element, 
                           const SubControlVolume &scv) const
{
    ...
}

Neumann boundary condition values / boundary fluxes (only called on Neumann boundaries)

NumEqVector neumannAtPos(const GlobalPosition& globalPos) const
{
    NumEqVector values(0.0);
    if (injectionActive() && onInjectionBoundary(globalPos))
    {
        values[Indices::conti0EqIdx + FluidSystem::N2Idx] = -1e-4;
        values[Indices::conti0EqIdx + FluidSystem::H2OIdx] = 0.0;
    }

    return values;
}

NumEqVector is an array of equations (here, the size of the array is 2).

More general Neumann boundary function:

template<class ElementVolumeVariables, class ElementFluxVariablesCache>
NumEqVector neumann(const Element& element,
                    const FVElementGeometry& fvGeometry,
                    const ElementVolumeVariables& elemVolVars,
                    const ElementFluxVariablesCache& elemFluxVarsCache,
                    const SubControlVolumeFace& scvf) const
{
    ...
}

Initial conditions:

PrimaryVariables initialAtPos(const GlobalPosition& globalPos) const
{
    PrimaryVariables values(0.0);
    const Scalar densityW
        = FluidSystem::H2O::liquidDensity(temperature(), 1.0e5);
    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;
}

Source/sink terms:

NumEqVector sourceAtPos(const GlobalPosition &globalPos) const
{ return NumEqVector(0.0); }

More general source/sink function:

template<class ElementVolumeVariables>
NumEqVector source(const Element &element,
                   const FVElementGeometry& fvGeometry,
                   const ElementVolumeVariables& elemVolVars,
                   const SubControlVolume &scv) const
{
    ...
}

Spatial Parameters definition

Spatial parameters

*spatialparams*.hh

template<class FVGridGeometry, class Scalar>
class InjectionSpatialParams
: public FVPorousMediumFlowSpatialParamsMP<
    FVGridGeometry, Scalar,
    InjectionSpatialParams<FVGridGeometry, Scalar>
> {
// parameters that are typically position-dependent
// e.g. porosity, permeability
};

Inherit from FVPorousMediumFlowSpatialParamsMP where FV: finite volumes, MP: multi-phase flow.

A function returning the instrinsic permeability \(K\):

auto permeabilityAtPos(const GlobalPosition& globalPos) const
{
    if (isInAquitard_(globalPos))
        return aquitardK_;
    return aquiferK_;
}

A function returning the porosity:

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

A function returning a parameterized constitutive law describing fluid-matrix interactions (\(p_c(S_w)\), \(k_r(S_w)\)):

auto fluidMatrixInteractionAtPos(const GlobalPosition& globalPos) const
{
    if (isInAquitard_(globalPos))
        return makeFluidMatrixInteraction(aquitardPcKrSwCurve_);
    return makeFluidMatrixInteraction(aquiferPcKrSwCurve_);
}

Set the (constant) temperature field in the domain:

Scalar temperatureAtPos(const GlobalPosition& globalPos) const
{
    return 273.15 + 20.0;  // 20°C
}

Runtime parameters

Dune INI syntax ([Group] and Key = Value pairs)

params.input

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

[Problem]
Name = test

See Part I: Runtime parameters for details.

Main program and main function

Main program

Each problem has one main file (test_name.cc or main.cc)
Here: 2pmain.cc and 2pnimain.cc
Contains the main function (mandatory in C/C++ programs)

Main function

Calling Dumux::initialize is mandatory at the beginning of each DuMu^x program to make sure the environment is correctly set up.

int main(int argc, char** argv)
{
    // initialize MPI+X backends (mandatory)
    Dumux::initialize(argc, argv);

Parse runtime parameters:

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

See Part I: Runtime parameters for details.

Define an alias for the test problem type tag

using namespace Dumux;
using TypeTag = Properties::TTag::Injection2pCC;

The tag (tag alias) is used to extract types and values via the property system (details on properties in Part II: Property system).

Create the computational grid:

// create a grid (take parameters from input file)
GridManager<GetPropType<TypeTag, Properties::Grid>> gridManager;
gridManager.init();

// obtain the grid and a grid view
const auto& grid = gridManager.grid();
const auto& leafGridView = grid.leafGridView();

More details on the grid in Part I: Grid.

GridGeometry and Problem instantiation:

// create the finite volume grid geometry
// (represents the chosen discretization scheme)
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);

Initial solution and secondary variables:

// 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);

Setting up VTK output:

// initialize the vtk output module
VtkOutputModule<GridVariables, SolutionVector>
    vtkWriter(*gridVariables, x, problem->name());

// optionally add a velocity post-processor
using VelocityOutput
    = GetPropType<TypeTag, Properties::VelocityOutput>;
vtkWriter.addVelocityOutput(
    std::make_shared<VelocityOutput>(*gridVariables));

// add input/output defaults for the chosen model
using IOFields = GetPropType<TypeTag, Properties::IOFields>;
IOFields::initOutputModule(vtkWriter);

Some typical main function structures:

Steady-state linear problems
Steady-state non-linear problems
Time-dependent non-linear problems

Assembling the system for a steady-state linear problem:

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

// 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);

and solve it using an iterative method:

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

(*r) *= -1.0; // solve Ax = -r to save update and copy
linearSolver->solve(*A, x, *r);

Simplified code for steady-state linear problem using Dumux::LinearPDESolver

auto assembler = ...;  // construct as before
auto linearSolver = ...;  // construct as before

using Solver = LinearPDESolver<Assembler, LinearSolver>;
Solver solver(assembler, linearSolver);

// assemble & solve
solver.solve(x);

For steady-state non-linear problems, use Dumux::NewtonSolver:

auto assembler = ...;  // construct as before
auto linearSolver = ...;  // construct as before

using Solver = NewtonSolver<Assembler, LinearSolver>;
Solver solver(assembler, linearSolver);

// linearize & solve
solver.solve(x);

For time-dependent problems, pass a TimeLoop instance and a solution vector carrying the solution on the previous/initial time level to the assembler:

SolutionVector xOld = x;

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, xOld // <-- new arguments
);

// assemble linear/non-linear problem as before
// ...

The skeleton of a time loop:

timeLoop->start(); do {

    // Time loop body

} while (!timeLoop->finished());

Solving a time step:

// time loop
timeLoop->start(); do
{
    // linearize & solve
    nonLinearSolver.solve(x, *timeLoop);
    // make the new solution the old solution
    xOld = x;
    gridVariables->advanceTimeStep();

Advancing the time loop to the next step:

    // 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());
// print final report
timeLoop->finalize(leafGridView.comm());

Build system (CMake)

Typical C++ build steps

configure
compile
install (copy programs to some location)

Step 3 is usually skipped by DuMu^x developers. But you can of course “install” Dune modules and DuMu^x.

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 portable and supports cross-platform compilation
Build file generator: Support various backends called generators (default: GNU Makefiles)
Control files written in the CMake language

Build system - CMakeLists.txt

Every directory in a project contains a file called CMakeLists.txt, which is written in the CMake language
Upon configure, the top-level CMakeLists.txt is executed
Subfolders are included with the add_subdirectory(...) function

Build system - configuring

Configure build-time compiler parameters / linking information
Create “targets” that can be build to create executables
Find and link external libraries / software dependencies

Build system - dunecontrol

For Dune, CMake is triggered via the dunecontrol script

./dune-common/bin/dunecontrol [options] <command>

useful commands:
- configure: configure libraries, variables, paths, files
- make all: build libraries
- all: includes configure and make all

Build system - dunecontrol

For Dune, CMake is triggered via the dunecontrol script

./dune-common/bin/dunecontrol [options] <command>

useful options:
- --opts=<optionfile.opts> specify e.g. compiler flags; DuMu^x provides dumux/cmake.opts.
- --build-dir=<build directory> specify path for out-of-source build; Default: every module contains its own build directory called build-cmake/
- --only <module> only apply command to a specific module/s (comma-separated list)

Build system - dunecontrol

You may have to reconfigure (rerun dunecontrol) whenever a dependency changes or a Dune library is updated. CMake caches. To reconfigure, you may need to remove the cache first:

./dune-common/bin/dunecontrol bexec rm -r CMakeFiles CMakeCache.txt

removes the CMake cache files from all Dune modules. (Use --only <module> to only remove caches of specific modules.)

Build system - basic CMake functions

Use add_subdirectory for recursively adding subdirectories
Subdirectory must contain a CMakeLists.txt (can be empty)
Executables are added via add_executable(<name> source1 [source2 ...])
Tests are added via dumux_add_test which also add a test executable to the test suite
Symlinks can be added via dune_symlink_to_source_files

Build system

CMakeLists.txt for an example application: create a test target exercise_basic_2p by defining a name and source file:

dumux_add_test(NAME exercise_basic_2p
               SOURCES 2pmain.cc)

Build system

Another example: add extra compile definitions and commands (typical for DuMu^x regression tests):

dumux_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 definitions for an executable:

dune_symlink_to_source_files(FILES "params.input")
# or: add_input_file_links()

# 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

Useful resources:

CMake online documentation
Dune CMake online documentation for CMake functions available with Dune.

Parallelism

Distributed memory parallelism with MPI

MPI stands for Message Passing Interface
MPI manages communication between processes that run in parallel

Distributed memory parallelism with MPI

Main idea is the concept of domain decomposition
The grid manages the grid partitioning
Each local subdomain is assembled on an individual process (rank)
Solvers have to handle decomposed vectors/matrices
Most solvers in DuMu^x (iterative Krylov subspace methods based on dune-istl) are capable of MPI parallel solving

Distributed memory parallelism with MPI

Run with:

   mpirun -np [n_cores] [executable_name]

Distributed memory parallelism with MPI

Partitioned VTK output files:

Each rank creates its own *.vtu/*.vtp file
Combined into *.pvtu/*.pvtp files for each time step
*.pvd groups all *.pvtu/*.pvtp files into a time series

Shared-memory parallelism

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

Shared-memory parallelism

Restrict the number of threads

Number of threads can be restricted via the environment variable DUMUX_NUM_THREADS

    DUMUX_NUM_THREADS=2 ./executable

Shared-memory parallelism

Multi-threaded assembly

Control use of multi-threading during assembly in input file:

[Assembly]
Multithreading = false

Defaults to true.

Exercises

Basic application setup:

Go to Exercise basic

The main function:

Go to Exercise main file

DuMux applications

Example application

Simulation goal

Gas injection / immiscible two phase flow

Gas injection / immiscible two phase flow

Gas injection / immiscible two phase flow

A DuMux application

A DuMux application

Property specializations

Properties

Problem definition

Problem

Spatial Parameters definition

Spatial parameters

Runtime parameters

Runtime parameters

Main program and main function

Main program

Main function

Build system (CMake)

Typical C++ build steps

Build system - what is CMake?

Build system - CMakeLists.txt

Build system - configuring

Build system - dunecontrol

Build system - dunecontrol

Build system - dunecontrol

Build system - basic CMake functions

Build system

Build system

Build system

Build system

Parallelism

Distributed memory parallelism with MPI

Distributed memory parallelism with MPI

Distributed memory parallelism with MPI

Distributed memory parallelism with MPI

Shared-memory parallelism

Shared-memory parallelism

Shared-memory parallelism

Exercises

Exercises

DuMu^x applications

A DuMu^x application

A DuMu^x application