DuMu$^\text{x}$ 3 -- an open-source simulator for solving flow and transport problems in porous media with a focus on model coupling

We present version 3 of the open-source simulator for flow and transport processes in porous media DuMu$^\text{x}$. DuMu$^\text{x}$ is based on the modular C++ framework Dune (Distributed and Unified Numerics Environment) and is developed as a research code with a focus on modularity and reusability. We describe recent efforts in improving the transparency and efficiency of the development process and community-building, as well as efforts towards quality assurance and reproducible research. In addition to a major redesign of many simulation components in order to facilitate setting up complex simulations in DuMu$^\text{x}$, version 3 introduces a more consistent abstraction of finite volume schemes. Finally, the new framework for multi-domain simulations is described, and three numerical examples demonstrate its flexibility.


Introduction
DuMu x , Dune for multi-{phase, component, scale, physics, domain . . .} flow and transport in porous media, is a free and open-source simulator for flow and transport processes in porous media [1] (dumux.org).Its main intention is to provide a sustainable, consistent and modular framework for the implementation and application of porous media model concepts and constitutive relations.It has been successfully applied to greenhouse gas and CO 2 storage [2,3,4,5], radioactive waste disposal [6], environmental remediation problems [7], transport of therapeutic agents through biological tissue [8,9,10,11], fractured porous media [12,13,14,15], and subsurface-atmosphere coupling [16,17].For a more complete list of publications that have been achieved with the help of DuMu x , see dumux.org/publications.
DuMu x is based on the Distributed Unified Numerics Environment (Dune) [18,19,20], an open-source scientific numerical software framework for solving partial differential equations, and is thus part of a larger community that goes beyond the simulation of fluid-mechanical processes in porous media.Dune and DuMu x are written in C++, using modern C++ programming techniques and C++ template meta programming for efficiency and generic interfaces.The Dune core modules provide, among other things, multiple grid managers implementing a versatile common grid interface, linear algebra abstraction and an iterative solver back-end [21], as well as abstractions facilitating parallel computing.DuMu x is designed and developed as a Dune module depending on the Dune core modules and optionally interacts with a number of other Dune-based extension modules.The key features of DuMu x are its comprehensive library of multi-phase and multi-component flow and transport models, the flexible and modular fluid and material framework for constitutive relations [1], abstractions for finite volume discretization schemes (see Section 3), and its focus on model coupling (see Section 4).
This paper can be seen as a continuation of [1] which describes the code and the project at the beginning of the 2.X release series.In the remainder of this introductory section, we provide a brief chronology of the development of DuMu x and some information on activities and measures that go beyond the development of the mainline code base.Section 2 introduces the general structure of the code and its design principles.Section 3 goes into more details by explaining abstractions and software concepts for general finite volume schemes encountered in DuMu x .Section 4 is devoted to simulations coupling two or more computational domains and models.In Section 5, we briefly address new features in the new version DuMu x 3. Numerical examples that particularly highlight the strengths of the new multi-domain framework are presented in Section 6.Finally, Section 7 discusses current limitations and perspectives.

History
The following section briefly outlines the DuMu x project history and its development.More details can be found in [34].
Initial development and first release.The development of DuMu x started in January 2007 at the Department of Hydromechanics and Modelling of Hydrosystems (LH2) at the University of Stuttgart.After an evaluation of the PDE software frameworks available at that time, the decision fell to build the code on top of the C++ toolbox Dune [19,18].In March 2007, a Subversion repository was set up for DuMu x to host and control the code development.
From that point on, every new doctoral student and postdoc at the LH2 performed his/her modeling tasks by using and enhancing the program code in that repository.Up until today, most developers of DuMu x are associated to the LH2 working group.This naturally results in a continuous major goal of the project that is to provide a tool enabling all developers to perform their research and possibly also teaching tasks.In July 2009, DuMu x 1.0 was released under the GNU General Public License Version 2 (or later) [35].
The 2.X release series.DuMu x 1.0 consisted of a subset of the code stored in the Subversion repository because not everything in there was adequate for public release.Since selecting a subset of a private code base for public release proved to be rather impractical, the repository was split into a stable part dumux and a development part dumux-devel that was dependent on the stable part.Public read access to the repository of the stable part was granted.Since the code design still exhibited several shortcomings concerning dependencies, generality and modularity, a major refactoring of the code base was performed in the following 1.5 years.As a result, DuMu x 2.0 was released in February 2011 which also lead to the main reference [1].After this release, interfaces were kept stable for at least one release cycle in order to achieve more sustainability and security for the growing number of users and developers.
Following the tendency from centralized to distributed version control, the DuMu x Subversion repository was converted into Git repositories in September 2015, following the transition of the Dune framework to Git.GitLab was employed as a version control management system, and the Git repositories are publicly accessible at git.iws.uni-stuttgart.de.The release process has been streamlined and since release 2.4 in October 2013, DuMu x has been released semiannually in spring and autumn every year, the last release of the 2.X series being 2.12 in late 2017.Since 2.7, every release tarball is uploaded to Zenodo, thereby receiving a DOI, as for example [36].
Transition to DuMu x 3.During the 2.X release series, several new features were added to the code base.Many of these additions fell in line with the main intention of DuMu x to be a framework for the implementation of porous-media model concepts and constitutive relations by actually providing such implementations.However, some more central additions had to be rather forced into the code base and proved to be inefficient, inconsistent with the original design ideas or increasingly difficult to maintain.This lead to a new major release cycle that was initiated by branching off the main development line in November 2016.Due to the large amount of changes and their extent into all parts of the code base, the requirement of backward compatibility was dropped.Right after the release of DuMu x 2.12 in December 2017, the development branch was integrated back into the main line and an alpha release was published.It took another year before DuMu x 3.0 finally was released in December 2018.Starting with the release of DuMu x 3.1 scheduled for October 2019, the aforementioned semiannual release cycle is resumed.

DuMu x as a framework
In the following, we address aspects going beyond the mainline code features and development.After discussing quality assurance and reproducibility, we introduce the two modules dumux-course and dumux-lecture.This is followed by a brief sketch of the Open Porous Media (OPM) initiative and looking at community building.Again, more details are presented in [34].
Quality assurance and reproducibility.To assure the quality of the developed software, DuMu x is currently accompanied by about 400 unit and system tests automatically built and run by a BuildBot CI server at git.iws.uni-stuttgart.de/buildbot after each commit to the master branch.We regularly assess that the tests cover a large number of lines of code and work towards increasing the coverage1 .One guideline for developers is that every newly added feature has to be accompanied by a corresponding test and a comprehensive documentation.In order to achieve reproducibility of the computational results, the project DuMu x -Pub has been initiated in 2015.It provides a set of tools for the researcher to extract the code that has been used in a publication into a separate Dune module, together with an install script to download and compile the code, including all necessary dependencies.Since 2015, every journal publication at the LH2 as well as every bachelor, master and doctoral thesis has to be accompanied by such a module.All resulting modules are published at git.iws.uni-stuttgart.de/dumux-pub.First efforts have been undertaken to provide also a complete runtime environment in form of a Docker container, see for example git.iws.uni-stuttgart.de/dumux-pub/Koch2017a.These efforts are continued as part of the ongoing DFG project "Sustainable infrastructure for the improved usability and archivability of research software on the example of the porous-media-simulator DuMu x " (SusI), which also aims to provide browser frontends to operate corresponding DuMu x Docker containers., where computer exercises using DuMu x are an essential part.The module dumux-lecture contains all the example applications employed in these lectures, with most of the applications belonging to the lecture "Multiphase Modeling".Each application has its own folder containing problem setup and spatial parameter specifications, as well as an input file where typically those runtime parameters can be specified that are of educational value for the problem.Each example is accompanied by explanations and tasks.
The Open Porous Media (OPM) initiative.In 2009, the OPM initiative was born with the principal objective "to develop a simulation suite that is capable of modeling industrially and scientifically relevant flow and transport processes in porous media and bridge the gap between the different application areas of porous media modeling" [37] (opm-project.org).Currently, the main focus of OPM is on reservoir engineering with a particular emphasis on being able to compete with proprietary industry-standard tools.As such, the black-oil simulator Flow is a main product of OPM.Flow is built upon the OPM modules opm-material and ewoms, both of which originated from a fork of DuMu x 2.2 in March 2012.In turn, DuMu x has a suggested dependency on the OPM module opm-grid which enables DuMu x users to use corner-point grids, the de facto standard in the petroleum industry.

Community building.
By going open source, external users of DuMu x have been welcome ever since the initial release in 2009.In order to get in contact with the users, a first DuMu x user meeting was held in Stuttgart in June 2015.26 participants were counted, 10 of them were external users.For attracting new users, a first DuMu x course was given in October 2017, followed by a second one in July 2018, and a short course at the InterPore conference 2019.In May 2019, the kick-off workshop for the aforementioned project SusI took place in Stuttgart, where the project investigators and 12 external users focused on identifying current shortcomings and possible remedies associated with the usability of DuMu x .With GitLab, all technical possibilities are available for users to upload contributions to the code base in form of merge requests as well as for developers to review, discuss, improve and possibly integrate them.
However, contributions from outside the LH2 are very scarce until now.Current notable exceptions are members from the neighboring working group VEGAS and the Institute of Applied Analysis and Numerical Simulation (IANS) at the University of Stuttgart as well as researchers from the Federal Waterways Engineering and Research Institute (BAW) in Karlsruhe.

Structure and design principles
DuMu x is designed as a research code framework.Foremost, this means that it is designed with an emphasis on modularity.DuMu x user code is usually developed in a separate Dune module that lists DuMu x as its dependency.To maintain flexibility for a user, DuMu x follows the principle that all components of a simulation should be easily replaceable with a new implementation, without modifying code in the DuMu x module itself.For example, it is simple to change from one of the implemented capillary-pressure relationships to a custom function, to modify the numerical flux computation, and to change from a cell-centered to a vertex-centered finite volume discretization.DuMu x profits from the modular design of Dune.For example, an unstructured grid implementation can be exchanged for an efficient structured grid implementation by changing a single line in the user code.
To maintain computational efficiency, much of the modularity is realized through the use of C++ templates and generic programming techniques such as traits and policies.Moreover, the developers of DuMu x try to follow wellknown object-oriented design principles and the separation of algorithms and data structures, as for example in the design of the DuMu x material framework [38].DuMu x 3, which requires a C++14-compatible compiler, increasingly makes use of modern C++ features such as lambdas, type deduction, or smart pointers, to the benefit of usability, modularity, and efficiency.
In the DuMu x environment, the term model is used to describe a system of coupled partial differential equations (PDEs) including constitutive equations needed for closure.Many models, in particular PDEs describing general nonisothermal multi-component multi-phase flow processes in porous media, are already implemented in DuMu x .Furthermore, a user can choose from a variety of constitutive laws to describe closure relations, and multiple fluid and solid systems and components.More importantly, the code design also facilitates using custom implementations.
The main components of a DuMu x simulation are represented in the code by corresponding C++ classes.When using one of the implemented models outof-the-box, a user usually implements at least two such classes in addition to the program's main function.The Problem class defines boundary conditions, initial conditions (if necessary), and volumetric source terms (if any), by implementing a defined class interface.The SpatialParams class defines the spatial distribution of material parameters, such as porosity, permeability, or for example, parameters for the van Genuchten water retention model (if applicable).The Problem class stores a pointer to an instance of the SpatialParams class.Moreover, each simulation currently creates at least one tag (a C++ struct) which is used to specify several compile time options (properties) of the simulation by means of partial template specialization.Code Example 1 shows how to set properties for a newly defined tag.Properties can be extracted as traits of such a tag in other parts of the code.Depending on the chosen model, more of such properties have to be defined and examples are provided for all existing models in DuMu x .User defined models can provide their own properties.A rather significant improvement in DuMu x 3 concerns the way the program's main function is written.In DuMu x 2, the main function called a predefined start function that determined the program flow.Modifications to the program flow were made possible by numerous hook functions.For example, the Problem base class implemented a hook function postTimeStep that could be overloaded by the user's Problem implementation, injecting the notion of time and program flow into the Problem class; a hook function addVtkOutputFields injected a notion of file I/O.This design concept introduced many dependencies between classes and (as shown here for the Problem class) lead to the violation of the well-known single responsibility principle, hindering modular design.Furthermore, due to the concept of hook functions, it was very difficult to write simulations outside the predefined program flow in start.For example, all simulations were considered transient and non-linear.Consequently, solving stationary equations was often realized by solving a single time step and making sure that the assembled storage term is zero.Solving linear equations was realized by executing a single iteration of a Newton method.On the positive side, a user had the possibility to realize most changes within a single header file.
In contrast, the program's main function in DuMu x 3 contains all of the main steps of the simulation, see Code Example 2. Therefore, to solve a stationary PDE, the TimeLoop instance (l.41) and the time loop (l.53ff) can simply be removed from the program flow.Similarly, the NewtonSolver class can be replaced by a LinearPDESolver to solve a PDE.Finally, the I/O logic and other modifications to the main program flow can be directly written inside the main function.As a result, the main program flow becomes arguably more transparent, which has been confirmed by feedback from the DuMu x user community.For more information on the steps of a DuMu x simulation, we refer the interested reader to the DuMu x handbook [39].

Abstractions and concepts for general finite volume schemes
One of the most important abstractions in DuMu x is the grid geometry.The concept partly existed in DuMu x 2 but was significantly redesigned and has been casted into an object-oriented representation in DuMu x 3. A grid geometry is a wrapper around a grid view on a Dune grid instance.Dune grids are general hierarchical grids and grid views provide read-only access to certain parts of the grid.In particular, a leaf grid view is a view on all grid entities without descendants in the hierarchy (which are not refined), thus covering the whole domain (in sequential simulations) or the part of the grid assigned to a single processor (in parallel simulations), while a level grid view is a view on all entities of a given level of the refinement hierarchy.The grid geometry constructs, from such a grid view, all the geometrical and topological data necessary to evaluate the discrete equations resulting from a given finite volume scheme.This abstraction allows us to implement many different finite volume schemes in a unified way.In the following sections, we will present and motivate the mathematical abstractions behind the grid geometry concept, describe their realizations in the code in form of C++ classes and provide an exemplary code snippet that illustrates how to use them.

Finite volume discretization
Let us consider a domain Ω ⊂ R n , 1 ≤ n ≤ 3, with boundary ∂Ω, which is further decomposed into two subsets on which Dirichlet and Neumann boundaries are specified, that is ∂Ω = ∂Ω D ∪ ∂Ω N .Furthermore, let us consider a general conservation equation, ) where m = m(u) is the conserved quantity depending on the primary variable u, ψ = ψ(u) is a flux term and q = q(u) a source term.All terms may non-linearly depend on u.Equations (3.1a) to (3.1d) constitute a non-linear problem in u with initial and boundary conditions.We introduce the primary grid M with elements E ∈ M, such that Ω h = E∈M E is a discrete approximation of Ω.Furthermore, we introduce a tessellation T of Ω, such that Ω h = K∈T K, where each K ∈ T is a control volume with measure |K| > 0. The control volumes K do not necessarily need to coincide with the elements E. However, for cell-centered finite volume schemes, usually T ≡ M. Each control volume is partitioned into sub-control-volumes κ, where κ = K is one admissible partition.Moreover, the boundary of the control volume ∂K consists of a finite number of faces f ⊂ ∂K which are either inner faces f I = ∂K ∩ ∂L, where L ∈ T is an adjacent control volume, or boundary faces f B = ∂K ∩ ∂Ω h .Depending on the discretization method, it may be useful to partition the faces f into sub-(control-volume)-faces σ.We denote by E K the entire set of such sub-faces on ∂K, forming a disjoint partition such that ∂K = σ∈E K σ.Accordingly, T K denotes the set of sub-control volumes embedded in a control volume K with K = κ∈T K κ.
Integration of Eq. (3.1a) over a control volume and application of the divergence theorem yields where n is a unit outer normal vector on ∂K.Using the concept of sub-controlvolumes and sub-control-volume-faces, allows to reformulate Eq. (3.2) as Let us now replace the exact flux over sub-face σ by the approximation F K,σ ≈ σ ψ • n dΓ and approximate the volume integrals in Eq. (3.3) to arrive at the discrete, control-volume-local formulation of Eq. (3.1a), where the approximate integrals m κ and q κ depend on the specific finite volume method.The time derivative is approximated by a backward or forward finite difference 2 , determining the time level t or t + ∆t at which the flux and source terms F K,σ and q are evaluated.The expressions for the discrete fluxes F K,σ depend on the actual underlying finite volume scheme.Fore more detailed information, we refer to [40,41] for the box scheme, and [42] for two-point and multi-point-flux-approximation finite volume schemes (tpfa and mpfa).
If a backward Euler time discretization is chosen, the nonlinear system corresponding to Eq. (3.4) on each control volume K is solved by Newton's method.In each Newton step n, we assemble the discrete PDE system in residual form where u now is the vector of primary variables for each degree of freedom, and A = ∂r ∂u is the Jacobian of the residual r, i.e. the discrete equation evaluated at u n .The residual for a control volume K ∈ T is given by the left-hand-side of Eq. (3.4).In the default configuration, the derivatives in the Jacobian matrix are approximated by numerical differentiation in DuMu x .

Element-wise assembly
Dune grids and the connectivity information provided by the grid are rather element-centered (elements are grid entities of co-dimension 0).For this reason, it is natural in Dune to choose an element-wise assembly of the residual.The choice of T and hence the sub-control-volumes κ and sub-control-volume-faces σ is motivated by the goal of a convenient implementation of an element-wise assembly for a given control volume scheme.Figure 1 depicts the configurations of control volumes and sub-faces as chosen in DuMu x , exemplarily for the tpfa, mpfa and the box scheme.A corresponding illustration for the mac scheme can be found, for example, in [43].As illustrated in Fig. 1, the sub-controlvolumes for the box scheme are chosen such that the residual associated with a degree of freedom on a vertex can be split into contributions from all neighboring elements.For the mpfa-o scheme, the sub-faces are defined such that each subface can be associated to a grid vertex.This is motivated by the fact that the expressions for the discrete fluxes across the sub-faces depend on the degrees of freedom within so-called interaction regions, constructed around the grid vertices [44].

Representation in the software implementation
Sub-control-volumes κ and sub-control-volume-faces σ, as introduced in Section 3.1, are represented in the code by objects of the classes SubControlVolume and SubControlVolumeFace. Depending on the caching model, which can be selected by specifying a template argument, these objects are either stored in the GridGeometry instance, or constructed on-the-fly on a per-element basis during the assembly.Furthermore, DuMu x has the notion of a local view on the grid geometry.The local view is an element-centered view that gives access to all κ and σ on an element, as well as all connectivity and geometry information to assemble storage, flux, and source terms.For example, for a tpfa scheme, the neighboring control volumes and associated primary and secondary variables have to be accessed to calculate the fluxes over the control-volume interfaces.Code Example 3 shows how to assemble the total mass contained in a domain for a given grid M. The presented code is identical for all introduced finite volume schemes.Similar to the range-based-for syntax to iterate over all elements of a grid view provided by Dune, the local view on the grid geometry in DuMu x makes it possible to iterate over all sub-control-volumes associated with an element in a concise and readable way (l.17).Analogously, a range-generator for iterating over all sub-control-volume-faces associated with an element is provided.Using the presented abstraction, the following discretization schemes have been implemented in DuMu x : tpfa, mpfa-o, mpfa-l, mac (staggered grid), finite volumes with non-linear two-point and multi-point flux approximation (nl-tpfa, nl-mpfa), and mimetic finite differences (with additional degrees of freedom on the grid faces) [45,46], from which tpfa, box, mpfa-o, and the mac scheme are available in the latest DuMu x version.

Multi-domain simulations
Model coupling is one of the key features in DuMu x [47,9] and has been a driving force in the development of DuMu x 3, which features a new, consistent framework for implementing coupled models, that is, models composed of multiple coupled sub-models.Figure 2 shows different types of model coupling that can be realized in this framework.In this context, the term multi-domain sim- ulation covers a broad class of simulations featuring model coupling.The idea is to solve a coupled PDE system, of which a subset (a subdomain model) can be formulated on a different domain, potentially with different dimensionality, and may be discretized using different meshes or different spatial discretization schemes.The framework is designed to minimize software coupling, such that two existing DuMu x models can be coupled without modifications in the core components.To be minimally invasive, only the Problem class is required to be slightly modified to become a sub-problem.Sub-problems define their own initial and boundary conditions, and source term.However, boundary and source terms may depend on quantities of other domains.For the data transfer, we introduce the concept of a coupling manager which is represented in the code by classes that implement the CouplingManager interface.A coupling manager has to provide the information which degrees of freedom in one domain are coupled to which degrees of freedom in another domain (the coupling stencil ).To this end, every coupling manager has to implement the function shown in Code Example 4. Code Example 4: The member function couplingStencil has to be implemented by all deriving coupling manager implementations.The template arguments i and j are the indices of one pair of coupled subdomains.They are deduced from the two index objects domainI and domainJ, passed as arguments to the function.In case two subdomain are not coupled, the function is required to return a reference to an empty stencil vector.For domain i, an instance of the element is passed, for which the residual in the element-wise assembly is to be computed.The function returns a vector of all indices of degrees of freedom coupled to one of the degrees of freedom associated with element elementI.
Moreover, the coupling manager has to transfer data between the sub-problems.This general concept can be used to implement a wide class of coupling schemes.Finding the coupling stencils often involves intersecting two grids.For this purpose, DuMu x provides efficient intersection algorithms based on axis-aligned bounding box volume hierarchy data structures [48,49] computed for Dune grids.Code Example 5 shows how to intersect two given grid geometries and use the connectivity information to construct a coupling stencil, that can be used, for example, in a non-conforming embedded fracture model 3 .Code Example 5: Sample code which intersects two grid geometries (see Section 3.1) and uses the resulting connectivity information to construct coupling stencils for a cell-centered finite volume scheme.A bulk grid (bulk) is intersected with an overlapping grid with lower dimension (lowDim) which may discretize a three-dimensional rock domain and a network of two-dimensional fractures.(The code snippet is identical for 3D-2D and 2D-1D.)An intersection is Γ = E b ∩ E l , where E b (targetEntity) is an element of the bulk grid and E l (domainEntity) an element of the lower-dimensional embedded grid.If there exist multiple such intersections with identical geometry (e.g. a fracture element coincides with the intersection of two bulk elements), they are represented in the code by a single intersection object with access to all elements ("entities") associated with the intersection ("neighbors").The interface of the glue object is similar to that implemented in the Dune module dune-grid-glue [50] which implements an advancing front algorithm instead of an algorithm based on spatial data structures used here.
Furthermore, DuMu x provides an assembler class for multi-domain models, which assembles the discrete PDE system in residual form where A i is the Jacobian of the discrete PDE system for subdomain i, and C ij is the coupling Jacobian with derivatives of residuals of domain i with respect to degrees of freedom of domain j, C ij = ∂ri ∂uj .The assembler class and the matrix class are generic, so that the sub-vectors u i and sub-matrices A i can themselves have a block structure, and support an arbitrary number of subdomains.The block structure can be exploited, for example, for constructing preconditioners for a monolithic solver, or to algebraically implement schemes where the subdomain systems are solved successively in an iterative algorithm.
The presented multi-domain concept has been successfully used to implement models with coupled flow and transport processes in vascularized brain tissue [10,51], root-soil interaction models in the vadose zone [9], flow and transport models for fractured rock systems [15,52], coupled porous medium flow and atmosphere flow (Darcy-Navier-Stokes) at the soil surface [43], and a model that couples a pore-network model with a Navier-Stokes model [53].

New features in DuMu x 3
Numerous new features are added in DuMu x 3 in comparison with the 2.X series.We briefly mention the most important changes.A complete redesign of many high-level class abstractions such as assembler, linear and non-linear solvers, grid readers and grid geometry, and file I/O leads to more readable and flexible main functions, see also Section 2. The numerous models are improved in terms of code reuse and modularity, so that code duplication is minimized and readability improved.In addition to the fluid system concept, a solid system concept is introduced, which facilitated the implementation of new models including mineralization or precipitation of substances that potentially modify the porous matrix structure [54,55].Porous material properties such as intrinsic permeability and porosity can be implemented to depend (linearly or non-linearly) on the primary variables.We generalized thermal and chemical non-equilibrium models to be combinable with any porous medium model.Multi-component diffusion can now be modeled by Maxwell-Stefan diffusion.A versatile implementation of cell-centered mpfa-o scheme is now usable with all models [15].The Navier-Stokes models are redesigned to use a mac scheme on a staggered grid [56,43,53], including a Reynolds-averaged Navier-Stokes model [57,58] with a variety of turbulence models (e.g.k-, k-ω), as well as a second-order upwind scheme.We can now solve problems based on the twodimensional shallow water equations.Finally, the multi-domain module adds the functionality of versatile model coupling as described in Section 4.

Numerical examples
In the following section, three numerical examples demonstrate the flexibility of the multi-domain framework in DuMu x 3.In Section 6.1, free flow over a porous medium is modeled by coupling a Navier-Stokes model to a pore-network model.The domains are coupled at a common interface.Section 6.2 shows a simulation of two-phase flow in a fractured rock matrix.The fracture flow is computed on lower-dimensional domains conforming with the control-volume faces of the three-dimensional rock matrix domain discretization, which allows to model highly conductive fractures as well as impermeable fractures.The example shows the differences between a tpfa and an mpfa-o finite volume scheme for a rock matrix with anisotropic permeability.Finally, an example of root water uptake and tracer transport is given.The roots are represented by a network of tubes embedded into the soil matrix.The mass exchange between the two non-conforming domains is realized with adequate source terms.The source code for all examples as well as instructions for the reproduction of the results can be found in the DuMu x -pub module to this publication at git.iws.uni-stuttgart.de/dumux-pub/dumux2019.

Coupling a free flow model with a pore-network model
This example is adapted from [53], where a coupled model of free channel flow adjacent to a porous medium is presented with a detailed model description.We model transient multi-component flow over a random porous structure in a two-dimensional model domain.The channel flow is governed by the Navier-Stokes equations, neglecting gravity and dilatation [59], with the fluid density ρ, the fluid velocity v, the fluid pressure p, and the dynamic viscosity µ.A molar balance equation for each component κ, models advective and diffusive transport of the fluid components, where x κ is the component mole fraction, ρ m the fluid molar density, and D the binary diffusion coefficient.Diffusive fluxes are described by Fick's law.
For modeling flow and transport in the porous medium, a pore-network model is used [60,61].The complex pore-scale geometry of the porous medium is transferred to a simplified geometry, the pore bodies and pore throats.For each component κ, Eq. (6.2) is formulated discretely on each pore body where the primary variables are located.The advective and diffusive fluxes are evaluated on the pore throats.We refer to [53] for further details.
Coupling conditions are formulated at the interface between the two models to ensure thermodynamic consistency.At the locations where no throat intersects with the interface, a no-flow/no-slip condition for the Navier-Stokes model is enforced.At the actual intersections, we prescribe the continuity of normal forces [62] resulting in a Neumann boundary condition for the free flow domain, where the superscripts ff and pnm mark quantities of the free-flow and porenetwork model, respectively.We use the tangential part of the average velocity within the throat at the boundary [v] pnm as a slip condition for the free flow, where [t] ff is a unit tangential vector to the coupling interface.The velocity within the throat is given by where Q t is the volume flow within the throat, A t is its respective cross-sectional area and n t is a unit vector parallel to the pore throat, pointing towards the coupling interface.Finally, we require the conservation of mass and enforce the continuity of mole fractions at the interface, The Navier-Stokes equations are discretized with a mac scheme on a staggered grid.The pore-network model is also implemented in DuMu x and will become part of stable code basis in an upcoming release.The pore-network is described as a one-dimensional network embedded in a two-dimensional domain, using the Dune grid implementation dune-foamgrid [63].We used the DuMu x multi-domain framework in order to achieve a fully monolithic coupling between the two sub-models.Since only elements on the domain boundary are coupled in this example, we employed a simplified intersection algorithm for creating the coupling stencils in comparison with Code Example 5. Instead of intersecting the entire grid, only the end points (boundary pores) of the pore-network grid geometry are intersected with the channel domain grid geometry to compute the coupling stencils.Furthermore, for the staggered grid discretization, we compute separate coupling stencils for the degrees of freedom located at cell centers and those located on cell faces.Newton's method is used to solve the nonlinear system of equations in combination with SuiteSparse's UMFPack [64] as a direct linear solver.Implementation details can be found in the folder dumux/multidomain/boundary located in the DuMu x repository and in the DuMu x -pub module accompanying this paper (see above).
The given example is discussed in detail in [53].Fluid density and viscosity are assumed constant, with ρ = 1 • 10 3 kg/m 3 and µ = 1 • 10 −3 Pa s.A tracer injected at the bottom of the pore network is transported upwards until it reaches the free-flow channel through which it leaves the system at its left or right sides where fixed pressures are set (see Fig. 3).All other sides of the channel are closed and no-flow/no-slip conditions hold.In the pore network, Dirichlet conditions for p and x κ are set at the bottom while Neumann no-flow boundaries are assigned to the lateral sides.Varying the pressure gradients in the free-flow channel yields three different scenarios.Figure 3 shows the resulting velocity fields where distinct preferential flow paths in the network and the influence of the inclined throats at the interface become visible.Figure 4 shows the temporal development of the concentration fields.At t 1 , an average mole fraction of 5 • 10 −4 is reached while t 2 corresponds to a value of 9 • 10 −4 .While the first two cases reach these points after 14 and 50 seconds, the higher pressure gradient in the channel for the case with Re = 55.24 repels the tracer, keeping it longer in the network.As there is no imposed flow for the first case, the tracer spreads equally to the left and right side of the channel while in the other two cases, the pressure gradient drives the tracer towards the right outlet.The formation of a boundary layer at the interface can be observed which becomes thinner for the higher Re case.

Two-phase flow through fractured porous media
The example shown in this section is inspired by an exercise of the DuMu xcourse 4 and considers the buoyancy-driven upwards migration of gaseous nitrogen in an initially fully water-saturated fractured porous medium.In this model, the fractures are assumed to represent very thin heterogeneities with substantially differing material parameters in comparison to the surrounding porous medium.The thin nature of these inclusions favors a dimension-reduced description of the fractures as entities of co-dimension 1, on which cross-section averaged PDEs are solved and appropriate coupling conditions describe the interaction with the surrounding porous medium.Such approaches have been widely reported in the literature for both single-phase flow (see e.g.[65,66,67]) and two-phase flow (see e.g.[68,69,70]) in fractured porous media.In the approach presented in this example, the two subdomains are not discretized independently, but it requires the facets of the higher-dimensional grid (bulk grid) to be conforming with the elements of the lower-dimensional grid for the fractures.Therefore, we currently rely on grid file formats that allow the extraction of both grids from a single file together with the connectivity information between them.The implementation in this example uses mesh files generated by Gmsh [71].This facilitates the determination of the coupling stencils, making it obsolete to intersect the grids, in contrast to the procedure presented in Code Example 5 for non-conforming methods.All classes and functions related to this approach can be found in the dumux/multidomain/facet folder of the DuMu x repository.For further details on the numerical scheme, we refer to [15,52].
The fracture network geometry used in this example is taken from [67] and is shown together with the boundary conditions in Fig. 5.The computational f .Making use of this notation, we state the governing system of equations as follows: ) for i ∈ {m, f} and β ∈ {H 2 O, N 2 }.Furthermore, ρ β and µ β are the density and viscosity of a fluid phase β, while p β,i , S β,i and k rβ,i denote the fluid phase pressure, saturation and relative permeability as prevailing in subdomain i. Correspondingly, K i and φ i denote the intrinsic permeability and the porosity in a subdomain, and v β,i is the respective Darcy velocity.Finally, K ⊥ f denotes the normal part of the permeability in the fracture domain and a is the aperture of the fracture.Equation (6.8a) states that Darcy's law is valid within all subdomains, while Eqs.(6.8b) and (6.8c) are the mass balance equations for the fluid phases in the matrix and fracture subdomain, respectively.Note that the term ρ β v β,m • n describes the jump in the normal flux of a phase β across the fracture and acts as an additional source or sink term, caused by the interaction with the surrounding matrix.On blocking fractures, flux and pressure continuity are enforced by means of Eq. (6.8d), while on conductive fractures, the pressure continuity condition Eq. (6.8e) is used, which assumes the pressure jump across the fracture to be negligible.The system of equations is completed by the following set of boundary conditions: Here, the subscript hs refers to hydrostatic conditions and the boundary saturation S D N2 is defined as depicted in Fig. 5.In both domains, the water pressure p H2O and the nitrogen saturation S N2 are used as primary variables and the system is closed via the constitutive relationships S H2O,i +S N2,i = 1 and p N2,i = p H2O,i + p c,i , where p c,i = p c,i (S H2O,i ) denotes the capillary pressure as a function of the water saturation.For this example, we use the van Genuchten-Mualem model for capillary pressure and relative permeability [72,73,74].An overview over the material parameter distributions used in the subdomains is given in Table 1.VG: van Genuchten model [72].
Note that the permeability of a subdomain is defined as where R is the two-dimensional rotation matrix.In addition to that, the aperture is set to a = 0.05 m for all fractures.The relationships of the water density and viscosity on temperature and pressure are taken from [75], while the nitrogen gas viscosity is modeled after [76].The density of gaseous nitrogen is described with the ideal gas law.The problem has been solved using the tpfa, mpfa-o and the box scheme.Switching between the different schemes is realized by changing at most two lines in the code.The grid used in this example consists of 20 910 2-dimensional and 421 1-dimensional elements.Figure 6 depicts the results for the box and the mpfa-o scheme, as well as the difference in the solution between tpfa and mpfa-o at the final simulation time of 75 000 s.
It can be seen that as nitrogen migrates upwards driven by buoyancy, it avoids entering the blocking fractures, which act as both capillary and hydraulic barriers, and preferentially flows around them.On the other hand, after entering conductive fractures, the nitrogen rapidly flows vertically and accumulates at the fracture tips due to capillary forces and the limited transport capacity of the surrounding medium.Furthermore, the nitrogen gas accumulates below the upper no-flow boundary.The differences between tpfa and mpfa-o highlight the importance of using consistent schemes on unstructured grids and in presence of full tensor permeabilities.This example shows that the model is able to qualitatively capture the hydraulic effects of fractures acting as both hydraulic and capillary barriers as well as conduits.For a more quantitative assessment of the performance of this model we refer to [15,52].

Root-soil interaction
In the last example, we simulate root-water uptake including tracer transport.The employed model concepts and numerical methods are described in detail in [9], where a similar example is discussed.A comparable but slightly simpler example is included in the DuMu x test suite 5 .
Root-water uptake is modeled in a cylindrical soil domain Ω with a radius of 5 cm and a height of 10 cm containing a 2-week-old white lupin root system (with center-line skeleton Λ), reconstructed from MRI data [77].The soil domain Ω is discretized with an unstructured tetrahedral grid refined around the root system, using dune-uggrid, while Λ is represented by an independent grid of line segments forming a root network, using dune-foamgrid [63].At the root collar a transpiration rate of r T = min{2.15• 10 −8 kg s −1 , r T,c } is pre-scribed, where r T,c is the transpiration rate for which the root collar pressure is p r,c = −1.4MPa (wilting point pressure).The soil domain contains a tracer of initially uniformly distributed concentration (x κ w = 3 • 10 −7 ).The tracer (binary diffusion coefficient in water D κ w = 2.3 • 10 −9 m s −2 ) does not enter the roots.As a consequence, it is expected to accumulate at the locations of the largest root water uptake.The soil flow is modeled with the Richards equation, and the flow inside the root xylem is modeled using a bundle-of-tubes approach, resulting in a Darcy-type flow model, cf.[9].The nonlinear PDE system where with the average soil pressure ps = 1 2π 2π 0 p s | R dθ, and an analogous definition for xκ w,s [9,78], is solved for water pressure p w and tracer mole fraction x κ w for a period of 3 d with a maximum time step size of 1 h.The symbols and parameter values are given in Table 2.The subscripts r and s denote root and soil quantities wherever they need to be distinguished.Equations (6.11a) and (6.11b) are spatially discretized using the box scheme, while Eq.(6.11c) is discretized with the tpfa scheme.The corresponding coupling manager is generic, in the sense that it can compute the coupling stencils between the root and the soil domain if the discretization schemes in both domains differ.The box scheme is superior to the tpfa scheme in the soil domain for the presented case, because it usually uses less degrees of freedom on tetrahedral grids, and it is well-known that the tpfa scheme is not consistent on non-K-orthogonal grids [45].In comparison with Code Example 5, the procedure to compute the coupling stencil is slightly more complex for this example, since the integration of the source term q w requires the evaluation of quantities on the root segment surface, resulting in non-local stencils.For details on the discretization method for this mixeddimension embedded problem, we refer to [9].The required classes and functions can be found in the DuMu x repository's dumux/multidomain/embedded folder.
The resulting spatial distribution of water and the tracer is shown in Fig. 7.A tracer accumulation in close vicinity to the roots is evident.This type of simulation may present a valuable tool to interpret experimental tracer data.Due to the non-linear relation of the tracer concentration distribution to the water uptake rate and its dependence on water distribution, it is otherwise difficult to obtain quantitative result on local root water uptake rates from tracer concentration measurements.

Current limitations and perspectives
As a research code under active development, DuMu x currently has certain limitations, from which many may be resolved in future versions.We briefly mention some of the limitations.
Multi-domain simulations in DuMu x currently do not run in parallel when involving multiple grids.In Dune, the Dune grid implementations are responsible for managing MPI-based distributed memory parallelism.When having two or more grid instances, data has to be communicated between grids.First steps towards such a feature within the Dune framework have been undertaken with the dune-grid-glue module [50].Efficient load balancing in complicated multi-domain, or mixed-dimensional setups is challenging and is an active field of research targeted by codes such as the MOOSE framework [32], or preCICE [79].
DuMu x currently supports forward and backward Euler time discretizations only.Implementing other time discretization schemes requires some non-trivial refactoring of the assembly process.
As discussed in Section 2, DuMu x uses a C++ programming technique based traits of tags to define what we call properties of a model.Due to the dependency on a tag, many of such properties (often corresponding to C++ types) are grouped together.Some classes have such a tag as a template argument (see FVAssembler in Code Example 2, l. 44).While this significantly reduces the number of template arguments of a class, it leads to an nontransparent way of injecting dependencies.Furthermore, FVAssembler<Tag1> and FVAssembler<Tag2> are different types, even if all type traits (properties) extracted from the tags are identical.Modularity and re-usability of such classes is impeded.Thus, the use of tags as template arguments has been significantly reduced in DuMu x 3 in comparison with previous versions.A medium-term goal is to replace all such classes with classes with explicit dependencies.The property technique is also the main reason, why it is currently difficult to add useful Python wrappers (a feature desired by many DuMu x users) for DuMu x .
The focus of DuMu x is on providing many usable and extensible models, and less on flexible linear algebra.While efficient and flexible data structures and easy-to-use and extensible, preconditioned linear solvers are available through dune-istl [21], DuMu x currently does not facilitate algebraic manipulations of linear equation systems, such as different orderings of the Jacobian matrix, quasi-Newton schemes, or complex linear solver strategies based on matrix decomposition.
The future development of DuMu x will continue to reflect the most recent advances in the field of porous medium research within the DuMu x community.For example, the goal of a current project is the improvement of solvers for free-flow and shallow water models coupled with porous medium flow models.Pore-network models for single-and multi-phase flow will be included into the DuMu x framework in an upcoming release, including static and dynamic approaches.In collaboration with the development team of preCICE [79], the model coupling capabilities of DuMu x with internal and external modules will be further improved.The usability and archivability of DuMu x will be further improved within the aforementioned project SusI which constitutes a step towards reproducible research.
DuMu x course.The dumux-course module has been created for a DuMu x course offered 2018 in Stuttgart.It is continuously improved and enhanced along with the mainline development.While the module is used in conventional DuMu x courses with participants attending personally, it is also designed to be used independently by everyone who would like to learn how to work with DuMu x .All course exercises are documented and contain task descriptions and solutions.Their contents range from a very basic first building and running experience to fairly advanced applications based on coupled model equations.The module also contains the slides from the course providing broader context and background information.DuMu x lecture.The Department of Hydromechanics and Modelling of Hydrosystems at the University of Stuttgart (LH2) offers Master level courses for different study programs (Environmental and Civil Engineering, Water Resources Engineering and Management, Simulation Technology, or Computational Mechanics of Materials and Structures) using InheritsFrom = std :: tuple < OneP , CCTpfaModel >; }; } // end namespace TTag // set the Dune grid type template < class Tag > struct Grid < Tag , TTag :: Example > { using type = Dune :: UGGrid <3 >; }; // set the problem type template < class Tag > struct Problem < Tag , TTag :: Example > { using type = UserProblem ; }; // set the spatial parameters type template < class Tag > struct SpatialParams < Tag , TTag :: Example > { using type = U s e r S p a t i a l P a r a m s ; }; // set the fluid system type template < class Tag > struct FluidSystem < Tag , TTag :: Example > { using type = FluidSystems :: OnePLiquid < double , Components :: H2O < double > >; }; } // end namespace Properties } // end namespace Dumux Code Example 1: An example of a property setting for DuMu x 3. The simulation is configured to use a one-phase porous medium model (OneP), a cell-centered two-point-flux-approximation discretization scheme (CCTpfaModel), and liquid water as a fluid.The UserProblem and UserSpatialParams types are implementation-and user-defined.Header includes are omitted for brevity.The types attached to the Example tag, can be extracted as traits, see GetPropType in Code Example 2.

Figure 1 :
Figure1: Illustration of the configurations of control volumes and sub-faces on a computational grid for different discretization schemes.The tpfa and mpfa scheme use the grid elements as control volumes (E ≡ K ≡ κ), while in the box scheme the control volume is constructed around the grid vertices by connecting the barycenters of all adjacent geometrical entities.The partition of the boundary ∂K by the sub-faces σ is illustrated by depicting the corners of each sub-face.
const double porosity = 0.4; const double density = 1000.0;double total WaterMas s = 0.0; // iterate over all elements of the leaf grid view of the grid for ( const auto & element : elements ( leafGridView ) ) { // construct a local view on the grid geometry auto fvGeometry = localView (* f vGridGe ometry ) ; // the view has to be bound to the current element // this is a no -op if the geometry is cached fvGeometry .bind ( element ) ; // iterate over all sub -control -volumes in the local view // and accumulate the mass for ( const auto & scv : scvs ( fvGeometry ) ) total WaterMa ss += saturation [ scv .dofIndex () ] * porosity * density * scv .volume () ; } std :: cout << " The total water mass is " << tot alWater Mass << " kg / s " << std :: endl ; Code Example 3: The code to compute the total water mass in the domain Ω h given a discrete water saturation field with values for each degree of freedom, stored in the vector saturation.Most variables have been introduced in Code Example 2.

Figure 2 :
Figure 2: Different types of model coupling in DuMu x multi-domain simulations.(a) multiphysics models (on the same grid), (b) multiple non-overlapping domains with sharp conforming or non-conforming interface, (c) multiple overlapping domains with different discretizations, (d) conforming and non-conforming (embedded) mixed-dimensional domains (1D-2D, 1D-3D, 2D-3D).The different coupling modes can also be combined.Typical mixeddimensional simulations also solve multi-physics problems, or use different discretization schemes in the subdomains.The number of subdomains is not limited to two.

/
/ required member function of every Co up l in gM an a ge r class // return the coupling stencil for element elementI ( domain i ) // with respect to degrees of freedom in domain j template < std :: size_t i , std :: size_t j > const std :: vector < std :: size_t >& c ou pl i ng St en c il ( Dune :: index_constant <i > domainI , const Element <i >& elementI , Dune :: index_constant <j > domainJ ) const ;

/
/ which degrees of freedom of the other domain does the element residual depend on ?using Stencils = std :: size_t < std :: set < std :: size_t > >; Stencils lowDimCouplingStencils , b u l k C o u p l i n g S t e n c i l s ; // resize to the number of grid elements l o w D i m C o u p l i n g S t e n c i l s .resize ( l o w D i m G r i d G e o m e t r y .gridView () .size ( /* codim = */ 0) ) ; b u l k C o u p l i n g S t e n c i l s .resize ( b u l k G r i d G e o m e t r y .gridView () .size ( /* codim = */ 0) ) ; // intersect two grid geometries // we choose domain ( first argument ) : lowDim , target ( second argument ) : bulk // see < dumux / multidomain / glue .hh > const auto glue = makeGlue ( lowDimGridGeometry , b u l k G r i d G e o m e t r y ) ; // interate over all intersections for ( const auto & is : intersections ( glue ) ) { // the element index of the lowDim element of this intersection const auto domainIdx = l o w D i m G r i d G e o m e t r y .elementMapper () .index ( is .domainEntity (0) ) ; // there might be multiple bulk elements associated with this intersection for ( unsigned int i = 0; i < is .n u m T a r g e t N e i g h b o r s () ; ++ i ) { // the element index of the bulk element of this intersection const auto targetIdx = b u l k G r i d G e o m e t r y .elementMapper () .index ( is .targetEntity ( i ) ) ; // insert target -domain index pair into the respective stencil l o w D i m C o u p l i n g S t e n c i l s [ domainIdx ]. insert ( targetIdx ) ; b u l k C o u p l i n g S t e n c i l s [ targetIdx ]. insert ( domainIdx ) ; } }

Figure 3 :
Figure 3: Velocity fields for the three scenarios.Re is based on the averaged velocity within the channel.Note the different color scale for the network in the third scenario (bottom).pout = 0. Figure adapted from [53] (license: CC BY 4.0).

pin = 0 Figure 4 :
Figure 4: Distribution of the mole fraction x t for the three different scenarios at different times.Figure adapted from [53] (license: CC BY 4.0).

Figure 5 :
Figure 5: Domain and Dirichlet boundary conditions for the two-phase flow example through fractured porous media.The subscript hs refers to hydrostatic pressure conditions.

Figure 6 :
Figure 6: Nitrogen saturation distribution at the final simulation time t = 75 000 s obtained with the mpfa-o scheme (upper left), the box scheme (upper right) and the tpfa scheme (lower left) for the example application of two-phase flow through a fractured porous medium.The lower right image shows the difference in the saturations obtained with the mpfa-o and the tpfa scheme.

Figure 7 :
Figure 7: Simulation of root-water uptake of a white lupin and simultaneous tracer transport in the soil, shown at t = 3 d.Three horizontal cuts through the soil domain are shown.The tracer, with mole fraction x κw , is not taken up by the roots and accumulates, particularly where the root water uptake rate is highest.On the bottom slice, the water saturation Sw is shown.The saturation slightly decreases close to the roots.Its spatial gradient depends on the flow resistances in soil and root, the current water distribution, and the prescribed transpiration rate r T .

Table 1 :
Parameter choices for the two-phase flow through fractured porous media example.

Table 2 :
Symbols and parameter values for the root-soil interaction example.