Mixed-dimensional problem

A strength of Dune-MMesh is its application to mixed-dimensional problems. So, let us consider a mixed-dimensional Burger’s equation on a domain with a lower-dimensional interface \(\Gamma\) that represents a physical domain with aperture \(d > 0\).

Find \(u, u_\Gamma\) s.t \begin{align} \renewcommand{\jump}[1]{[\mskip-5mu[ #1 ]\mskip-5mu]} u_t + \operatorname{div} \mathbf{F}(u) &= 0, &&\quad \text{in } \Omega, \\ (d u_\Gamma)_t + \operatorname{div} c \mathbf{F}(u_\Gamma) &= \jump{\mathbf{F}(u) \cdot \mathbf{n}}, &&\quad \text{in } \Gamma, \end{align} where \(\mathbf{F}(u) = \frac{1}{2} u^2 \mathbf{a}\) with \(\mathbf{a} = [1,0]^T\) and \(c > 0\) is a speed scaling for the interface problem.

We impose the interior boundary condition for the bulk problem \begin{align} F(u) &= F(u_\Gamma), \quad \text{on } \Gamma, \\ \end{align} and consider inflow from left on the T-junction grid.

[1]:
from dune.grid import reader
from dune.mmesh import mmesh

file = "grids/tjunction.msh"

gridView = mmesh((reader.gmsh, file), 2)
igridView = gridView.hierarchicalGrid.interfaceGrid
[2]:
from ufl import *
from dune.ufl import Constant

d = Constant(0.01, name="d")
a = as_vector([1, 0])
c = Constant(100, name="c")

def F(u):
    return 0.5 * u**2 * a

def upwind(u_l, u_r, n):
    return inner( conditional( inner(a, n) > 0, F(u_l), F(u_r) ), n )

dt = 0.25
tau = Constant(dt, name="tau")

Bulk problem

We use a finite volume space and define the UFL form of the bulk problem. Hereby, we used the interface indicator I to distinguish between interface and non-interface facets.

[3]:
from dune.fem.space import finiteVolume
space = finiteVolume(gridView)

u  = TrialFunction(space)
uu = TestFunction(space)

x = SpatialCoordinate(space)
n = FacetNormal(space)

left = conditional(x[0] < 1e-6, 1, 0)
uD = 1

uh = space.interpolate(0, name="uh")
uhOld = uh.copy()

from dune.mmesh import interfaceIndicator
I = interfaceIndicator(igridView)

A  = (u - uhOld) * uu * dx
A += tau * upwind(u('+'), u('-'), n('+')) * jump(uu) * (1-I)*dS
A += tau * upwind(u, uD, n) * uu * left * ds

Interface problem

Similarly, we define the UFL form of the interface problem.

[4]:
ispace = finiteVolume(igridView)

u_g  = TrialFunction(ispace)
uu_g = TestFunction(ispace)

n_g = FacetNormal(ispace)

uh_g = ispace.interpolate(0, name="uh_g")
uhOld_g = uh_g.copy()

A_g  = d * (u_g - uhOld_g) * uu_g * dx
A_g += tau * c * upwind(u_g('+'), u_g('-'), n_g('+')) * jump(uu_g) * dS

Coupling condition

The coupling condition can be incorporated using the trace and skeleton functionality.

[5]:
from dune.mmesh import skeleton, trace, normals

u_gamma = avg(skeleton(uh_g))

A += tau * upwind(u('+'), u_gamma, n('+')) * uu('+') * I*dS
A += tau * upwind(u('-'), u_gamma, n('-')) * uu('-') * I*dS


traceu = trace(uh)
inormal = normals(igridView)

A_g -= tau * upwind(traceu('+'), u_g, -inormal) * uu_g * dx
A_g -= tau * upwind(traceu('-'), u_g, -inormal) * uu_g * dx

Monolithic solution strategy

Dune-MMesh provides coupled solution strategies to solve bulk and interface schemes together, either iteratively or monolithically. Here, we choose the monolithicSolve method that implements a mixed-dimensional Newton method.

[6]:
from dune.fem.scheme import galerkin
scheme   = galerkin([A   == 0])
scheme_g = galerkin([A_g == 0])

from dune.mmesh import monolithicSolve
def solve():
    uhOld.assign(uh)
    uhOld_g.assign(uh_g)
    monolithicSolve(schemes=(scheme, scheme_g), targets=(uh, uh_g), verbose=False)
[7]:
import matplotlib.pyplot as plt
from dune.fem.plotting import plotPointData as plot

fig, axs = plt.subplots(1, 4, figsize=(20,5))

for i in range(4):
    solve()

    plot(uh, figure=(fig, axs[i]), gridLines=None, clim=[0,1])
    plot(uh_g, figure=(fig, axs[i]), linewidth=0.04, clim=[0,1], colorbar=None)
../_images/examples_mixeddim_12_0.png