Solving two Poisson problems on subdomains of a parent mesh yielding non-deterministic success

General
1

I’m trying to solve a pretty straightforward problem:

Find u_0 \in V_0(\Omega_0) and u_1 \in V_1(\Omega_1) such that

\begin{align} (\nabla u_0, \nabla v_0)_{\Omega_0} &= (1, v_0)_{\Omega_0} \\ (\nabla u_1, \nabla v_1)_{\Omega_1} &= (1, v_1)_{\Omega_1} \end{align}

for all v_0 \in V_{0,E}(\Omega_0) and v_1 \in V_{1,E}(\Omega_1). Here V_{\cdot, E}(\cdot) corresponds to the space modified to accommodate homogeneous Dirichlet boundary data.

The following code is my attempt to approximate the solution of the above problem.

from mpi4py import MPI
import dolfinx
import dolfinx.fem.petsc
import numpy as np
import ufl

# Parent mesh
mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 8, 8)

# Extract submeshes, mesh0 left, mesh1 right
cells = np.arange(mesh.topology.index_map(mesh.topology.dim).size_local, dtype=np.int32)
sub_cells0 = cells[dolfinx.mesh.compute_midpoints(mesh, mesh.topology.dim, cells)[:, 0] <= 0.5]
sub_cells1 = cells[dolfinx.mesh.compute_midpoints(mesh, mesh.topology.dim, cells)[:, 0] >= 0.5]
sub_mesh0, mapping_entity0, mapping_vertex0, mapping_geom0 = dolfinx.mesh.create_submesh(mesh, mesh.topology.dim, sub_cells0)
sub_mesh1, mapping_entity1, mapping_vertex1, mapping_geom1 = dolfinx.mesh.create_submesh(mesh, mesh.topology.dim, sub_cells1)

# Create domain markers
sub_mt0 = dolfinx.mesh.meshtags(mesh, mesh.topology.dim, sub_cells0, np.ones_like(sub_cells0))
sub_mt1 = dolfinx.mesh.meshtags(mesh, mesh.topology.dim, sub_cells1, np.ones_like(sub_cells1))

# Create spaces on each subdomain and mixed element
V0 = dolfinx.fem.functionspace(sub_mesh0, ("CG", 1))
V1 = dolfinx.fem.functionspace(sub_mesh1, ("CG", 1))
W = ufl.MixedFunctionSpace(V0, V1)

# Integration measures on each subdomain
dx0 = ufl.Measure("dx", domain=mesh, subdomain_data=sub_mt0)
dx1 = ufl.Measure("dx", domain=mesh, subdomain_data=sub_mt1)

# FE formulation
u0, u1 = ufl.TrialFunctions(W)
v0, v1 = ufl.TestFunctions(W)

a = ufl.inner(ufl.grad(u0), ufl.grad(v0)) * dx0
a += ufl.inner(ufl.grad(u1), ufl.grad(v1)) * dx1

L = ufl.inner(dolfinx.fem.Constant(sub_mesh0, 1.0), v0) * dx0
L += ufl.inner(dolfinx.fem.Constant(sub_mesh1, 1.0), v1) * dx1

# Extract block system
a_blocked = ufl.extract_blocks(a)
L_blocked = ufl.extract_blocks(L)

# Create BCs
sub_mesh0.topology.create_connectivity(sub_mesh0.topology.dim, sub_mesh0.topology.dim - 1)
sub_mesh0.topology.create_connectivity(sub_mesh0.topology.dim - 1, sub_mesh0.topology.dim)
u0_bdry = dolfinx.fem.Function(V0)
u0_bc = dolfinx.fem.dirichletbc(
    u0_bdry, dolfinx.fem.locate_dofs_topological(
        V0, sub_mesh0.topology.dim - 1, dolfinx.mesh.exterior_facet_indices(sub_mesh0.topology)
))

sub_mesh1.topology.create_connectivity(sub_mesh1.topology.dim, sub_mesh1.topology.dim - 1)
sub_mesh1.topology.create_connectivity(sub_mesh1.topology.dim - 1, sub_mesh1.topology.dim)
u1_bdry = dolfinx.fem.Function(V1)
u1_bc = dolfinx.fem.dirichletbc(
    u1_bdry, dolfinx.fem.locate_dofs_topological(
        V1, sub_mesh1.topology.dim - 1, dolfinx.mesh.exterior_facet_indices(sub_mesh1.topology)
))

# Solve problem
bcs = [u0_bc, u1_bc]
u0h = dolfinx.fem.Function(V0)
u1h = dolfinx.fem.Function(V1)
problem = dolfinx.fem.petsc.LinearProblem(
    a_blocked,
    L_blocked,
    u=[u0h, u1h],
    bcs=bcs,
    entity_maps=[mapping_entity0, mapping_entity1],
    petsc_options={
        "ksp_type": "preonly",
        "pc_type": "lu",
        "pc_factor_mat_solver_type": "mumps",
        "ksp_error_if_not_converged": True,
        "ksp_monitor": None,
    },
    petsc_options_prefix="monolithic_",
)
problem.solve()

I notice that I’m getting non-deterministic success. Sometimes this solves giving me the expected result, and sometimes it does not, yielding the error:

Traceback (most recent call last):
  File "$SCRIPT_PATH/poisson_subdomain_coupled.py", line 65, in <module>
    problem = dolfinx.fem.petsc.LinearProblem(
        a_blocked,
    ...<11 lines>...
        petsc_options_prefix="monolithic_",
    )
  File "$LIB_PATH/python3.14/site-packages/dolfinx/fem/petsc.py", line 804, in __init__
    self._A = create_matrix(self._a, kind=kind)
              ~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^
  File "$LIB_PATH/python3.14/site-packages/dolfinx/fem/petsc.py", line 193, in create_matrix
    return _cpp.fem.petsc.create_matrix_block(_a, kind)
           ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^^^^^^^^^^
RuntimeError: Cannot insert rows that do not exist in the IndexMap.
Exception ignored while calling deallocator <function LinearProblem.__del__ at 0x12aee5380>:
Traceback (most recent call last):
  File "$LIB_PATH/python3.14/site-packages/dolfinx/fem/petsc.py", line 882, in __del__
    lambda obj: obj is not None, (self._solver, self._A, self._b, self._x, self._P_mat)
AttributeError: 'LinearProblem' object has no attribute '_solver'

Am I doing something silly in my original script?

>>> import dolfinx
>>> dolfinx.git_commit_hash
'338f32f692a1500571c910a29e69a3ff1d5c6549'
>>> dolfinx.__version__
'0.11.0.dev0'
2

I think I may have poorly defined the subdomain integration measures. Making the following change

# Integration measures on each subdomain
dx0 = ufl.Measure("dx", domain=mesh, subdomain_data=sub_mt0)
dx1 = ufl.Measure("dx", domain=mesh, subdomain_data=sub_mt1)

# Integration measures on each subdomain
dx0 = ufl.Measure("dx", domain=mesh, subdomain_data=sub_mt0, subdomain_id=1)
dx1 = ufl.Measure("dx", domain=mesh, subdomain_data=sub_mt1, subdomain_id=1)

Seems to give me consistently good results