This tutorial is part of multiphenicsx.

Tutorial 03: weak imposition of Dirichlet BCs by a Lagrange multiplier (interface problem)

In this tutorial we solve the problem

$$\begin{cases} -\Delta u = f, & \text{in } \Omega,\\ u = g, & \text{on } \partial\Omega, \end{cases}$$

where $f=1$, $g=0$ and $\Omega$ is a ball in 2D, using a domain decomposition approach for $\Omega = \Omega_1 \cup \Omega_2$, and introducing a lagrange multiplier to handle the continuity of the solution across the interface $\Gamma$ between $\Omega_1$ and $\Omega_2$.

The resulting weak formulation is: $$ \text{find } u_1, u_2, \lambda \in V(\Omega_1) \times V(\Omega_2) \times E(\Gamma) $$ s.t. $$ \int_{\Omega_1} \nabla u_1 \cdot \nabla v_1 dx + \int_{\Omega_2} \nabla u_2 \cdot \nabla v_2 dx + \int_{\Gamma} \lambda (v_1 - v_2) ds + \int_{\Gamma} \eta (u_1 - u_2) ds = \int_{\Omega_1} f v_1 dx + \int_{\Omega_2} f v_2 dx, \qquad \forall (v_1, v_2, \eta) \in V(\Omega_1) \times V(\Omega_2) \times E(\Gamma). $$

Equivalenty this equation can be written as a system of equations $$ \begin{align*} &\int_{\Omega_1} \nabla u_1 \cdot \nabla v_1 dx & & & &+ \int_\Gamma \lambda v_1 ds & &= \int_{\Omega_1} f v_1 ds & &\forall v_1 \in V(\Omega_1) \\ & & &\int_{\Omega_2} \nabla u_2 \cdot \nabla v_2 dx & &- \int_\Gamma \lambda v_2 ds & &= \int_{\Omega_2} f v_2 ds & &\forall v_2 \in V(\Omega_2) \\ &\int_\Gamma \eta u_1 ds & &- \int_\Gamma \eta u_2 ds & & & &= 0 & &\forall \eta \in E(\Gamma). \end{align*} $$

Boundary conditions on $\partial\Omega$ are embedded in $V(\Omega_i) \subset H^1(\Omega_i)$, $i = 1, 2$, and $E(\Gamma) \subset L^2(\Gamma)$.

This example is a prototypical case of problems containing interface restricted variables (the Lagrange multiplier, in this case).

import dolfinx.fem
import dolfinx.fem.petsc
import dolfinx.io
import gmsh
import mpi4py.MPI
import numpy as np
import petsc4py.PETSc
import ufl
import viskex

import multiphenicsx.fem
import multiphenicsx.fem.petsc

Geometrical parameters

r = 3
mesh_size = 1. / 4.

Mesh

gmsh.initialize()
gmsh.model.add("mesh")

p0 = gmsh.model.geo.addPoint(0.0, 0.0, 0.0, mesh_size)
p1 = gmsh.model.geo.addPoint(0.0, +r, 0.0, mesh_size)
p2 = gmsh.model.geo.addPoint(0.0, -r, 0.0, mesh_size)
c0 = gmsh.model.geo.addCircleArc(p1, p0, p2)
c1 = gmsh.model.geo.addCircleArc(p2, p0, p1)
l0 = gmsh.model.geo.addLine(p2, p1)
line_loop_left = gmsh.model.geo.addCurveLoop([c0, l0])
line_loop_right = gmsh.model.geo.addCurveLoop([c1, -l0])
semicircle_left = gmsh.model.geo.addPlaneSurface([line_loop_left])
semicircle_right = gmsh.model.geo.addPlaneSurface([line_loop_right])

gmsh.model.geo.synchronize()
gmsh.model.addPhysicalGroup(1, [c0, c1], 1)
gmsh.model.addPhysicalGroup(1, [l0], 2)
gmsh.model.addPhysicalGroup(2, [semicircle_left], 1)
gmsh.model.addPhysicalGroup(2, [semicircle_right], 2)
gmsh.model.mesh.generate(2)

Info    : Meshing 1D...
Info    : [  0%] Meshing curve 1 (Circle)
Info    : [ 40%] Meshing curve 2 (Circle)
Info    : [ 70%] Meshing curve 3 (Line)
Info    : Done meshing 1D (Wall 0.000332925s, CPU 0.00061s)
Info    : Meshing 2D...
Info    : [  0%] Meshing surface 1 (Plane, Frontal-Delaunay)
Info    : [ 60%] Meshing surface 2 (Plane, Frontal-Delaunay)
Info    : Done meshing 2D (Wall 0.0115016s, CPU 0.011678s)
Info    : 589 nodes 1201 elements

partitioner = dolfinx.mesh.create_cell_partitioner(dolfinx.mesh.GhostMode.shared_facet)
mesh, subdomains, boundaries_and_interfaces = dolfinx.io.gmshio.model_to_mesh(
    gmsh.model, comm=mpi4py.MPI.COMM_WORLD, rank=0, gdim=2, partitioner=partitioner)
gmsh.finalize()

# Create connectivities required by the rest of the code
mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)
mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim - 1)
mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim)

cells_Omega1 = subdomains.indices[subdomains.values == 1]
cells_Omega2 = subdomains.indices[subdomains.values == 2]
facets_partial_Omega = boundaries_and_interfaces.indices[boundaries_and_interfaces.values == 1]
facets_Gamma = boundaries_and_interfaces.indices[boundaries_and_interfaces.values == 2]

# Define associated measures
dx = ufl.Measure("dx", subdomain_data=subdomains)

viskex.dolfinx.plot_mesh(mesh)

error: XDG_RUNTIME_DIR is invalid or not set in the environment.
MESA: error: ZINK: failed to choose pdev
glx: failed to create drisw screen

viskex.dolfinx.plot_mesh_tags(mesh, subdomains, "subdomains")

error: XDG_RUNTIME_DIR is invalid or not set in the environment.

MESA: error: ZINK: failed to choose pdev
glx: failed to create drisw screen

viskex.dolfinx.plot_mesh_tags(mesh, boundaries_and_interfaces, "boundaries and interfaces")

error: XDG_RUNTIME_DIR is invalid or not set in the environment.
MESA: error: ZINK: failed to choose pdev
glx: failed to create drisw screen

Define dS measure over the interface

Interior facet integrals have no consistent ordering in dolfinx, resulting in the default assignment to "+" and "-" sides to be arbitrary. To restore a consistent ordering we define a custom dS measure.

integration_entities_on_Gamma = dolfinx.fem.compute_integration_domains(
    dolfinx.fem.IntegralType.interior_facet, mesh.topology, facets_Gamma, mesh.topology.dim - 1)
integration_entities_on_Gamma_reshaped = integration_entities_on_Gamma.reshape(-1, 4)
connected_cells_to_Gamma = integration_entities_on_Gamma_reshaped[:, [0, 2]]
subdomain_ordering = (
    subdomains.values[connected_cells_to_Gamma[:, 0]] < subdomains.values[connected_cells_to_Gamma[:, 1]])
if len(subdomain_ordering) > 0 and any(subdomain_ordering):
    integration_entities_on_Gamma_reshaped[subdomain_ordering] = integration_entities_on_Gamma_reshaped[
        subdomain_ordering][:, [2, 3, 0, 1]]
integration_entities_on_Gamma = integration_entities_on_Gamma_reshaped.flatten()

dS = ufl.Measure("dS", domain=mesh, subdomain_data=[(2, np.array(integration_entities_on_Gamma, dtype=np.int32))])
dS = dS(2)

Check correctness of the subdomain measure by integrating a piecewise function defined as (1, 1) on subdomain 1, and (2, 2) and subdomain 2.

DG = dolfinx.fem.functionspace(mesh, ("DG", 0, (mesh.topology.dim, )))

dg_function = dolfinx.fem.Function(DG)
dg_function.interpolate(lambda x: np.full((mesh.topology.dim, x.shape[1]), 1.0, dtype=np.float64), cells0=cells_Omega1)
dg_function.interpolate(lambda x: np.full((mesh.topology.dim, x.shape[1]), 2.0, dtype=np.float64), cells0=cells_Omega2)

n = ufl.FacetNormal(mesh)

dg_function_check_from_left_subdomain = mesh.comm.allreduce(
    dolfinx.fem.assemble_scalar(dolfinx.fem.form(ufl.inner(dg_function("-"), n("-")) * dS)), op=mpi4py.MPI.SUM)
dg_function_check_from_right_subdomain = mesh.comm.allreduce(
    dolfinx.fem.assemble_scalar(dolfinx.fem.form(ufl.inner(dg_function("+"), n("+")) * dS)), op=mpi4py.MPI.SUM)
print("Check correctness when integrating from the left subdomain", dg_function_check_from_left_subdomain)
print("Check correctness when integrating from the right subdomain", dg_function_check_from_right_subdomain)
assert np.isclose(dg_function_check_from_left_subdomain, 6., atol=1.e-10)
assert np.isclose(dg_function_check_from_right_subdomain, -12., atol=1.e-10)

Check correctness when integrating from the left subdomain 6.000000000000001
Check correctness when integrating from the right subdomain -12.000000000000002

With domain decomposition

# Define function spaces
V = dolfinx.fem.functionspace(mesh, ("Lagrange", 2))
V1 = V.clone()
V2 = V.clone()
M = V.clone()

# Define restrictions
dofs_V1_Omega1 = dolfinx.fem.locate_dofs_topological(V1, subdomains.dim, cells_Omega1)
dofs_V2_Omega2 = dolfinx.fem.locate_dofs_topological(V2, subdomains.dim, cells_Omega2)
dofs_M_Gamma = dolfinx.fem.locate_dofs_topological(M, boundaries_and_interfaces.dim, facets_Gamma)
restriction_V1_Omega1 = multiphenicsx.fem.DofMapRestriction(V1.dofmap, dofs_V1_Omega1)
restriction_V2_Omega2 = multiphenicsx.fem.DofMapRestriction(V2.dofmap, dofs_V2_Omega2)
restriction_M_Gamma = multiphenicsx.fem.DofMapRestriction(M.dofmap, dofs_M_Gamma)
restriction = [restriction_V1_Omega1, restriction_V2_Omega2, restriction_M_Gamma]

# Define trial and test functions
(u1, u2, l) = (ufl.TrialFunction(V1), ufl.TrialFunction(V2), ufl.TrialFunction(M))
(v1, v2, m) = (ufl.TestFunction(V1), ufl.TestFunction(V2), ufl.TestFunction(M))

# Define problem block forms
zero = dolfinx.fem.Constant(mesh, petsc4py.PETSc.ScalarType(0))
a = [[ufl.inner(ufl.grad(u1), ufl.grad(v1)) * dx(1), None, ufl.inner(l("-"), v1("-")) * dS],
     [None, ufl.inner(ufl.grad(u2), ufl.grad(v2)) * dx(2), - ufl.inner(l("+"), v2("+")) * dS],
     [ufl.inner(u1("-"), m("-")) * dS, - ufl.inner(u2("+"), m("+")) * dS, None]]
f = [ufl.inner(1, v1) * dx(1), ufl.inner(1, v2) * dx(2), ufl.inner(zero, m("-")) * dS]
a_cpp = dolfinx.fem.form(a)
f_cpp = dolfinx.fem.form(f)

# Define boundary conditions
dofs_V1_partial_Omega = dolfinx.fem.locate_dofs_topological(
    V1, boundaries_and_interfaces.dim, facets_partial_Omega)
dofs_V2_partial_Omega = dolfinx.fem.locate_dofs_topological(
    V2, boundaries_and_interfaces.dim, facets_partial_Omega)
bc1 = dolfinx.fem.dirichletbc(zero, dofs_V1_partial_Omega, V1)
bc2 = dolfinx.fem.dirichletbc(zero, dofs_V2_partial_Omega, V2)
bcs = [bc1, bc2]

# Assemble the block linear system
A = multiphenicsx.fem.petsc.assemble_matrix_block(a_cpp, bcs=bcs, restriction=(restriction, restriction))
A.assemble()
F = multiphenicsx.fem.petsc.assemble_vector_block(f_cpp, a_cpp, bcs=bcs, restriction=restriction)

# Solve
u1u2l = multiphenicsx.fem.petsc.create_vector_block(f_cpp, restriction=restriction)
ksp = petsc4py.PETSc.KSP()
ksp.create(mesh.comm)
ksp.setOperators(A)
ksp.setType("preonly")
ksp.getPC().setType("lu")
ksp.getPC().setFactorSolverType("mumps")
ksp.getPC().setFactorSetUpSolverType()
ksp.getPC().getFactorMatrix().setMumpsIcntl(icntl=7, ival=4)
ksp.setFromOptions()
ksp.solve(F, u1u2l)
u1u2l.ghostUpdate(addv=petsc4py.PETSc.InsertMode.INSERT, mode=petsc4py.PETSc.ScatterMode.FORWARD)
ksp.destroy()

<petsc4py.PETSc.KSP at 0x7f09594501d0>

# Split the block solution in components
(u1, u2, l) = (dolfinx.fem.Function(V1), dolfinx.fem.Function(V2), dolfinx.fem.Function(M))
with multiphenicsx.fem.petsc.BlockVecSubVectorWrapper(
        u1u2l, [V1.dofmap, V2.dofmap, M.dofmap], restriction) as u1u2l_wrapper:
    for u1u2l_wrapper_local, component in zip(u1u2l_wrapper, (u1, u2, l)):
        with component.x.petsc_vec.localForm() as component_local:
            component_local[:] = u1u2l_wrapper_local
u1u2l.destroy()

<petsc4py.PETSc.Vec at 0x7f08d8b70d10>

viskex.dolfinx.plot_scalar_field(u1, "u1")

error: XDG_RUNTIME_DIR is invalid or not set in the environment.
MESA: error: ZINK: failed to choose pdev
glx: failed to create drisw screen

viskex.dolfinx.plot_scalar_field(u2, "u2")

error: XDG_RUNTIME_DIR is invalid or not set in the environment.
MESA: error: ZINK: failed to choose pdev
glx: failed to create drisw screen

viskex.dolfinx.plot_scalar_field(l, "l")

error: XDG_RUNTIME_DIR is invalid or not set in the environment.
MESA: error: ZINK: failed to choose pdev
glx: failed to create drisw screen

Without domain decomposition

# Define trial and test functions
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

# Define problem forms
a_ex = ufl.inner(ufl.grad(u), ufl.grad(v)) * dx
f_ex = ufl.inner(1, v) * dx

# Define Dirichlet BC object on Gamma
dofs_V_partial_Omega = dolfinx.fem.locate_dofs_topological(
    V, boundaries_and_interfaces.dim, facets_partial_Omega)
bc_ex = dolfinx.fem.dirichletbc(zero, dofs_V_partial_Omega, V)

# Solve
u_ex = dolfinx.fem.Function(V)
problem_ex = dolfinx.fem.petsc.LinearProblem(
    a_ex, f_ex, bcs=[bc_ex], u=u_ex,
    petsc_options={
        "ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps",
        "mat_mumps_icntl_7": 4
    })
problem_ex.solve()
u_ex.x.petsc_vec.ghostUpdate(addv=petsc4py.PETSc.InsertMode.INSERT, mode=petsc4py.PETSc.ScatterMode.FORWARD)

viskex.dolfinx.plot_scalar_field(u_ex, "u")

error: XDG_RUNTIME_DIR is invalid or not set in the environment.
MESA: error: ZINK: failed to choose pdev
glx: failed to create drisw screen

Comparison and error computation

u_ex1_norm = np.sqrt(mesh.comm.allreduce(
    dolfinx.fem.assemble_scalar(dolfinx.fem.form(ufl.inner(u_ex, u_ex) * dx(1))), op=mpi4py.MPI.SUM))
u_ex2_norm = np.sqrt(mesh.comm.allreduce(
    dolfinx.fem.assemble_scalar(dolfinx.fem.form(ufl.inner(u_ex, u_ex) * dx(2))), op=mpi4py.MPI.SUM))
err1_norm = np.sqrt(mesh.comm.allreduce(
    dolfinx.fem.assemble_scalar(dolfinx.fem.form(ufl.inner(u_ex - u1, u_ex - u1) * dx(1))), op=mpi4py.MPI.SUM))
err2_norm = np.sqrt(mesh.comm.allreduce(
    dolfinx.fem.assemble_scalar(dolfinx.fem.form(ufl.inner(u_ex - u2, u_ex - u2) * dx(2))), op=mpi4py.MPI.SUM))
print("Relative error on subdomain 1", err1_norm / u_ex1_norm)
print("Relative error on subdomain 2", err2_norm / u_ex2_norm)
assert np.isclose(err1_norm / u_ex1_norm, 0., atol=1.e-10)
assert np.isclose(err2_norm / u_ex2_norm, 0., atol=1.e-10)

Relative error on subdomain 1 1.9093314700379406e-15
Relative error on subdomain 2 1.791646777841337e-15