In this tutorial we solve the problem
$$\begin{cases} -\Delta u = f, & \text{in } \Omega,\\ u = g, & \text{on } \partial\Omega, \end{cases}$$
where $\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 \in V(\Omega_1), u_2 \in V(\Omega_2), \eta \in 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 = 0, \qquad \forall v_1 \in V(\Omega_1), v_2 \in V(\Omega_2) $$ and $$ \int_{\Gamma} \eta (u_1 - u_2) ds = 0, \qquad \forall \eta \in E(\Gamma) $$ where 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
r = 3
mesh_size = 1. / 4.
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.000388364s, CPU 0.000801s) Info : Meshing 2D... Info : [ 0%] Meshing surface 1 (Plane, Frontal-Delaunay) Info : [ 50%] Meshing surface 2 (Plane, Frontal-Delaunay) Info : Done meshing 2D (Wall 0.0113973s, CPU 0.011614s) 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)
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)
dS = ufl.Measure("dS")(subdomain_data=boundaries_and_interfaces)
dS = dS(2) # restrict to the interface, which has facet ID equal to 2
viskex.dolfinx.plot_mesh(mesh)
viskex.dolfinx.plot_mesh_tags(mesh, subdomains, "subdomains")
viskex.dolfinx.plot_mesh_tags(mesh, boundaries_and_interfaces, "boundaries and interfaces")
# 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 0x7f0d10466890>
# 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 0x7f0d10466070>
viskex.dolfinx.plot_scalar_field(u1, "u1")
viskex.dolfinx.plot_scalar_field(u2, "u2")
viskex.dolfinx.plot_scalar_field(l, "l")
# 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")
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