multidomain.md

title: The DuMux Multidomain Framework

Coupling Multiple Domains

Coupling Multiple Domains

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

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)

/*!*
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];
    ...
}

Coupling Residual / Data Transfer (2/3)

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

Coupling Residual / Data Transfer (3/3)

template<class ElementVolumeVariables>
void pointSource(PointSource& source...) const
{
    ...
    const Scalar sourceValue = rho*(pressure3D - pressure1D)
                            /meanDistance*matrixPerm/mu;
    source = sourceValue*source.quadratureWeight()
            *source.integrationElement();
    ...
}

Coupling Residual (1/2)

 /*! 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!");
    ...
}

Coupling Residual (2/2)

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

Coupling Stencil

Coupling Stencil (1/2)

/*!
 * 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!");
    ...
}

Coupling Stencil (2/2)

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

Examples

Embedded Fracture Model

Root Water Uptake