multidomain.md

---
title: The DuMux Multidomain Framework
---

# Coupling Multiple Domains

## Coupling Multiple Domains
![](./img/multidomain.svg){width=900px}

## Design Goals

* Reuse existing DuMux models in the coupled setting
* Extensible to more than two subdomains
* Common structure for different coupling modes

# General Structure

## General Structure
![](./img/mdstructure.png)

## The Sub Problems

* any DuMux problem with
    - initial conditions
    - boundary conditions
    - associated DuMux model

* get a pointer to the coupling manager

## The Coupling Manager
* Transfer data from one subproblem to another
    - e.g. give me the soil pressure (from the well domain)
    - e.g. give me the radius of the embedded well (from the soil domain)

* Compute the coupling stencil
* Compute the coupling residual (numerical differentiation)

# Example: Embedded Fracture Model

## Coupling Residual

$$
\begin{aligned}
    \frac{\partial \rho \phi}{\partial t} - \nabla\cdot \left[ \frac{\rho}{\mu}\mathbf{K_m} (\nabla p_m - \rho \mathbf{g}) \right] - q \delta_\Gamma = 0 \quad \text{in} \, \Omega \\
    \frac{\partial \rho \phi}{\partial t} - \nabla_\tau\cdot \left[ \frac{\rho}{\mu}\mathbf{K_f} (\nabla p_f - \rho \mathbf{g}) \right] + q = 0 \quad \text{in} \, \Lambda \\
    q = \zeta (p_f - p_m) \quad \text{and} \quad \Gamma = \Omega \cap \Lambda
\end{aligned}
$$

## Coupling Residual / Data Transfer (1/3)

```c++
/*!*
Evaluates the point sources (added by addPointSources)
* for all phases within a given sub-control volume.!*/
template<class ElementVolumeVariables>
void pointSource(PointSource& source...) const
{
    // compute source at every integration point
    const Scalar pressure3D = this->couplingManager()
                .bulkPriVars(source.id())[Indices::pressureIdx];
    const Scalar pressure1D = this->couplingManager()
                .lowDimPriVars(source.id())[Indices::pressureIdx];
    ...
}
```
<span style="font-size: 0.4em; position: relative; top: -38px; color: gray;">File: `dumux/test/multidomain/embedded/2d3d/1p_1p/problem_matrix.hh`</span>

## Coupling Residual / Data Transfer (2/3)

```c++
template<class ElementVolumeVariables>
void pointSource(PointSource& source...) const
{
...
    // calculate the source
    const Scalar meanDistance = this->couplingManager()
                                .averageDistance(source.id());
    static const Scalar matrixPerm = getParamFromGroup<Scalar>("...");
    static const Scalar rho = getParam<Scalar>("...");
    static const Scalar mu = getParam<Scalar>("...")*rho;
...
}
```
<span style="font-size: 0.4em; position: relative; top: -38px; color: gray;">File: `dumux/test/multidomain/embedded/2d3d/1p_1p/problem_matrix.hh`</span>

## Coupling Residual / Data Transfer (3/3)

```c++
template<class ElementVolumeVariables>
void pointSource(PointSource& source...) const
{
    ...
    const Scalar sourceValue = rho*(pressure3D - pressure1D)
                            /meanDistance*matrixPerm/mu;
    source = sourceValue*source.quadratureWeight()
            *source.integrationElement();
    ...
}
```
<span style="font-size: 0.4em; position: relative; top: -38px; color: gray;">File: `dumux/test/multidomain/embedded/2d3d/1p_1p/problem_matrix.hh`</span>


## Coupling Residual (1/2)

```c++
 /*! evaluates the element residual of a coupled element of
    domain i which depends on the variables at the degree of
    freedom with index dofIdxGlobalJ of domain j */
template<std::size_t i, std::size_t j, class LocalAssemblerI>
decltype(auto) evalCouplingResidual(Dune::index_constant<i> domainI,
                                    const LocalAssemblerI& la,
                                    Dune::index_constant<j> domainJ,
                                    std::size_t dofIdxGlobalJ)
{
    static_assert(i != j, "A domain cannot be coupled to itself!");
    ...
}
```
<span style="font-size: 0.4em; position: relative; top: -38px; color: gray;">File: `dumux/dumux/multidomain/embedded/couplingmanagerbase.hh`</span>

## Coupling Residual (2/2)

```c++
decltype(auto) evalCouplingResidual(...)
{
    ...
    for (const auto& scv : scvs(fvGeometry))
    {
        auto couplingSource = this->problem(domainI)
                                .scvPointSources(...);
        couplingSource += this->problem(domainI).source(...);
        couplingSource *= -GridGeometry<i>::Extrusion::volume(...)
                          *curElemVolVars[scv].extrusionFactor();
        residual[scv.indexInElement()] = couplingSource;
    }
    return residual;
}
```
<span style="font-size: 0.4em; position: relative; top: -38px; color: gray;">File: `dumux/dumux/multidomain/embedded/couplingmanagerbase.hh`</span>

## Coupling Stencil

![](./img/disc.png)

## Coupling Stencil (1/2)

```c++
/*!
 * returns an iterable container of all indices of degrees of
 *freedom of domain j that couple with / influence the element
 *residual of the given element of domain i */
template<std::size_t i, std::size_t j>
const CouplingStencil<j>& couplingStencil(domainI,
                                          const Element<i>& element,
                                          domainJ) const
{
    static_assert(i != j, "A domain cannot be coupled to itself!");
    ...
}
```
<span style="font-size: 0.4em; position: relative; top: -38px; color: gray;">File: `dumux/dumux/multidomain/embedded/couplingmanagerbase.hh`</span>

## Coupling Stencil (2/2)

```c++
template<std::size_t i, std::size_t j>
const CouplingStencil<j>& couplingStencil(domainI,
                                          const Element<i>& element,
                                          domainJ) const
{
    ...
    if (couplingStencils(domainI).count(eIdx))
        return couplingStencils(domainI).at(eIdx);
    else
        return emptyStencil(domainI);
}
```
<span style="font-size: 0.4em; position: relative; top: -38px; color: gray;">File: `dumux/dumux/multidomain/embedded/couplingmanagerbase.hh`</span>

# Examples
## Embedded Fracture Model
![](./img/fracture.gif)

## Root Water Uptake
![](./img/root.gif)