problem.md

title: DuMu^x^ applications

Example application

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 \boldsymbol{v}_\alpha

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 = \varrho_\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

$\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}$

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

A DuMu^x^ application

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 FVGridGeometry
        = GetPropType<TypeTag, Properties::FVGridGeometry>;
    using Scalar
        = GetPropType<TypeTag, Properties::Scalar>;
    using type = InjectionSpatialParams<FVGridGeometry, 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).

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).

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

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

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)
• 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: