MPC backsubstitution in parallel

Dear all,

I am wondering whether mpc.backsubstitution also works in parallel. I solve the Stokes equations with periodic boundary conditions in 2D and the solution with a single rank is true. However, errors (only) occur at the slave nodes in parallel. These are the lines after solving the problems:

ksp.solve(b, func_sol.vector)
func_sol.x.scatter_forward()
mpc.backsubstitution(func_sol.vector)

ve, pr = func_sol.sub(0).collapse(), func_sol.sub(1).collapse()
ve.x.scatter_forward()
pr.x.scatter_forward()

Or am I missing something else after the solution of the problem? I appreciate any help! Thank you.

Here is the MWE:

import gmsh
import dolfinx
import dolfinx_mpc
from dolfinx.fem import FunctionSpace, VectorFunctionSpace, TensorFunctionSpace, Function, Constant, Expression, assemble_scalar, form
from dolfinx.io import XDMFFile
from dolfinx.mesh import locate_entities_boundary, meshtags
from dolfinx.io.gmshio import model_to_mesh

from ufl import grad, div, inner, Measure, TrialFunctions, TestFunctions, FiniteElement, VectorElement
import ufl
from petsc4py import PETSc
import numpy as np
from mpi4py import MPI

# MPI initialization
comm = MPI.COMM_WORLD
rank = MPI.COMM_WORLD.rank

#############################################
#              Mesh generation              #
#############################################
gmsh.initialize()
L = 1                       # length of rectangle
f = 0.3                     # volume fraction of inclusion
r = np.sqrt(f*L**2/np.pi)   # radius of inclusion

gdim = 2         # Geometric dimension of the mesh
fdim = gdim - 1  # facet dimension

if rank == 0:
    inclusion = gmsh.model.occ.addDisk(L / 2, L / 2, 0, r, r)
    matrix = gmsh.model.occ.addRectangle(0,0,0, L, L)
    whole_domain = gmsh.model.occ.fragment([(2, matrix)], [(2, inclusion)])
    gmsh.model.occ.synchronize()

    matrix_physical = gmsh.model.addPhysicalGroup(gdim, [matrix], tag = 1)
    inclusion_physical = gmsh.model.addPhysicalGroup(gdim, [inclusion], tag = 2)

    # Get the interface elements between rectangle and circle and tag
    interface_elements = gmsh.model.getBoundary([(gdim, inclusion_physical)])
    interface_physical = gmsh.model.addPhysicalGroup(fdim, (5, 5), tag = 9)

    background_surfaces = []
    inclusion_surfaces = []
    for domain in whole_domain[0]:
        com = gmsh.model.occ.getCenterOfMass(domain[0], domain[1])
        mass = gmsh.model.occ.getMass(domain[0], domain[1])
        if np.isclose(mass, np.pi * (r ** 2)): # identify inclusion by its mass
            inclusion_surfaces.append(domain)
        else:
            background_surfaces.append(domain)

    gmsh.model.mesh.field.add("Distance", 1)
    edges = gmsh.model.getBoundary(inclusion_surfaces, oriented=False)
    gmsh.model.mesh.field.setNumbers(1, "EdgesList", [e[1] for e in edges])

    gmsh.model.mesh.field.add("Threshold", 2)
    gmsh.model.mesh.field.setNumber(2, "IField", 1)
    gmsh.model.mesh.field.setNumber(2, "LcMin", 0.01)
    gmsh.model.mesh.field.setNumber(2, "LcMax", 0.03)
    gmsh.model.mesh.field.setNumber(2, "DistMin", 0.01)
    gmsh.model.mesh.field.setNumber(2, "DistMax", 0.1)
    gmsh.model.mesh.field.setAsBackgroundMesh(2)

    gmsh.option.setNumber("Mesh.Algorithm", 6)
    gmsh.model.mesh.generate(gdim)

domain, ct, ft = model_to_mesh(gmsh.model, comm, 0, gdim=gdim)
gmsh.finalize()

#############################################
#             Mesh generation finished      #
#############################################

dx = Measure('dx', domain=domain, subdomain_data=ct)
num_cells = domain.topology.index_map(domain.topology.dim).size_global

# elements and funcionspace
###########################

# Taylor-Hood elements and Mixed Functionspace
P1 = FiniteElement("CG", domain.ufl_cell(), 1)
P2 = VectorElement("CG", domain.ufl_cell(), 2)
TH = ufl.MixedElement([P2, P1])
V = FunctionSpace(domain, TH)

# define function, trial function and test function
(ve, pr) = TrialFunctions(V)
(dve, dpr) = TestFunctions(V)

# define material parameter
S = FunctionSpace(domain, ("DG", 0))
vis = Function(S)
vis = Function(S)
fluid_visc = ct.find(1)
vis.x.array[fluid_visc] = np.full_like(fluid_visc, 1, dtype=PETSc.ScalarType)
inclusion_visc = ct.find(2)
vis.x.array[inclusion_visc]  = np.full_like(inclusion_visc, 1e+5, dtype=PETSc.ScalarType)
vis.x.scatter_forward()

# boundary conditions
#####################
## set all fluid velocities at the interface to zero (no slip boundary condition)
# create collapsed functionspace
V_coll0, _ = V.sub(0).collapse()
v_bc = Function(V_coll0)
v_bc.x.set(0)
v_bc.x.scatter_forward()
# dofs of interface and dirichlet bcs
dofs_v_interf = dolfinx.fem.locate_dofs_topological((V.sub(0), V_coll0)  , ft.dim, ft.find(9))  # interface tag: 9
bcs_v_interf = dolfinx.fem.dirichletbc(v_bc, dofs_v_interf, V.sub(0))

## set pressures at the corners to zero
def boundary_node(x):
    return np.logical_or(
        np.logical_and(
            np.isclose(x[0], 0.),
            np.logical_or(
                np.isclose(x[1], .0),
                np.isclose(x[1], L),
            )
        ),
        np.logical_and(
            np.isclose(x[0], L),
            np.logical_or(
                np.isclose(x[1], 0.),
                np.isclose(x[1], L),
            )
        ),
    )

boundary_facets = locate_entities_boundary(domain, 0, boundary_node) # it is zero dimension for a point
dofs_corners = dolfinx.fem.locate_dofs_topological(V.sub(1), 0, boundary_facets)  # it is zero dimension for a point
bc_corners = dolfinx.fem.dirichletbc(PETSc.ScalarType(0.0), dofs_corners, V.sub(1))

bcs = [bcs_v_interf, bc_corners]

## periodic boundary conditions
pbc_directions = []
pbc_slave_tags = []
pbc_is_slave = []
pbc_is_master = []
pbc_meshtags = []
pbc_slave_to_master_maps = []

mpc = dolfinx_mpc.MultiPointConstraint(V)

def generate_pbc_slave_to_master_map(i):
    def pbc_slave_to_master_map(x):
        out_x = x.copy()
        out_x[i] = x[i] - L
        return out_x
    return pbc_slave_to_master_map

def generate_pbc_is_slave(i):
    return lambda x: np.isclose(x[i], L)

def generate_pbc_is_master(i):
    return lambda x: np.isclose(x[i], 0.0)

for i in range(gdim):

    pbc_directions.append(i)
    pbc_slave_tags.append(i + 2)
    pbc_is_slave.append(generate_pbc_is_slave(i))
    pbc_is_master.append(generate_pbc_is_master(i))
    pbc_slave_to_master_maps.append(generate_pbc_slave_to_master_map(i))

    facets = locate_entities_boundary(domain, fdim, pbc_is_slave[-1])
    arg_sort = np.argsort(facets)
    pbc_meshtags.append(meshtags(domain,
                                 fdim,
                                 facets[arg_sort],
                                 np.full(len(facets), pbc_slave_tags[-1], dtype=np.int32)))

N_pbc = len(pbc_directions)
for i in range(N_pbc):
    if N_pbc > 1:
        def pbc_slave_to_master_map(x):
            out_x = pbc_slave_to_master_maps[i](x)
            idx = pbc_is_slave[(i + 1) % N_pbc](x)
            out_x[pbc_directions[i]][idx] = np.nan
            return out_x
    else:
        pbc_slave_to_master_map = pbc_slave_to_master_maps[i]

    for ij in range(gdim):
        mpc.create_periodic_constraint_topological(V.sub(0).sub(ij), pbc_meshtags[i],
                                                   pbc_slave_tags[i],
                                                   pbc_slave_to_master_map,
                                                   bcs)

    mpc.create_periodic_constraint_topological(V.sub(1), pbc_meshtags[i],
                                                   pbc_slave_tags[i],
                                                   pbc_slave_to_master_map,
                                                   bcs)

if len(pbc_directions) > 1:
    # Map intersection(slaves_x, slaves_y) to intersection(masters_x, masters_y),
    # i.e. map the slave dof at (1, 1)/(1,1,1) to the master dof at (0, 0)/(0,0,0)
    def pbc_slave_to_master_map(x):
        out_x = x.copy()
        out_x[0] = x[0] - domain.geometry.x[:, 0].max()
        out_x[1] = x[1] - domain.geometry.x[:, 1].max()
        idx = np.logical_and(pbc_is_slave[0](x), pbc_is_slave[1](x))
        out_x[0][~idx] = np.nan
        out_x[1][~idx] = np.nan
        return out_x

    for ij in range(gdim):
        mpc.create_periodic_constraint_topological(V.sub(0).sub(ij), pbc_meshtags[1],
                                                   pbc_slave_tags[1],
                                                   pbc_slave_to_master_map,
                                                   bcs)

    mpc.create_periodic_constraint_topological(V.sub(1), pbc_meshtags[1],
                                               pbc_slave_tags[1],
                                               pbc_slave_to_master_map,
                                               bcs)

mpc.finalize()

# Variational problem
#####################
f_constant = Constant(domain, PETSc.ScalarType((1,0)))

## weak formulation
a = vis * inner(grad(ve), grad(dve)) * dx + pr * div(dve) * dx + div(ve) * dpr * dx
L = inner(f_constant, dve) * dx
a_f = form(a)
L_f = form(L)

# solution function
func_sol = Function(V)

# assemble matrices
A = dolfinx_mpc.assemble_matrix(a_f, mpc, bcs=bcs)
A.assemble()
b = dolfinx_mpc.assemble_vector(L_f, mpc)
dolfinx_mpc.apply_lifting(b, [a_f], [bcs], mpc)
b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
dolfinx.fem.petsc.set_bc(b, bcs)
b.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)

# Preconditioner matrix form
Pf = dolfinx.fem.form(ufl.inner(ufl.grad(ve), ufl.grad(dve)) * dx + pr * dpr * dx)
P = dolfinx_mpc.assemble_matrix(Pf, mpc, bcs=bcs)
P.assemble()

# solver settings
ksp = PETSc.KSP().create(domain.comm)
ksp.setOperators(A, P)
ksp.setType(PETSc.KSP.Type.MINRES)
pc = ksp.getPC()
pc.setType("hypre")
pc.setHYPREType("boomeramg")
pc.setUp()

# solving and backsubstitution to obtain the values at the slave degrees of freedom
ksp.solve(b, func_sol.vector)
func_sol.x.scatter_forward()
mpc.backsubstitution(func_sol.vector)

# collapse mixed functionspace
ve, pr = func_sol.sub(0).collapse(), func_sol.sub(1).collapse()
ve.x.scatter_forward()
pr.x.scatter_forward()

# average velocity to see error immediately
v_avg = domain.comm.allreduce(assemble_scalar(form(ve[0] * dx)), op=MPI.SUM)
if rank == 0:
    print("Average velocity in x direction:", v_avg)

file_ve = XDMFFile(domain.comm, "results_2D_iterative/velocity.xdmf", "w")
file_ve.write_mesh(domain)
file_ve.write_function(ve)
file_ve.close()

Backsubstitution calls scatter forward:

Depending on what version of dolfinx_mpc you use, i would advice using:
mpc.backsubstitution(func_sol)
as petsc sometimes does funny things with communication if one uses non petsc functions to do the scattering, such as func_sol.x.scatter_forward()

Thank you for your fast reply!
Also also tried to use mpc.backsubstitution(func_sol), but then I get this error message: AttributeError: 'Function' object has no attribute 'localForm'. I am using docker image dolfinx_mpc:v0.6.1.post1. Do you have any idea why this error occurs?

It is because I changed this back in march: API Change: Replace PETSc.Vec in python functions when possible (#53) · jorgensd/dolfinx_mpc@5572550 · GitHub

I’m currently preparing a release compatible with 0.7.0 of DOLFINx (where this would be included). I’ll ping you once it is published.

2 Likes

Okay perfect. Thank you very much!

Thank you for publishing the new dolfinx_mpc version. Unfortunately, this does not resolve the problem. The backsubstitution does not work for me in parallel.

After solving the problem, the commands are:

ksp.solve(b, func_sol.vector)
mpc.backsubstitution(func_sol)

# collapse mixed functionspace
ve = func_sol.sub(0).collapse()
pr = func_sol.sub(1).collapse()

# average velocity
v_avg = domain.comm.allreduce(assemble_scalar(form(ve[0] * dx)), op=MPI.SUM)
if rank == 0:
    print("Average velocity in x direction:", v_avg)

and this gives different results with one/two cpus:

mpirun -n 1 python3 MWE2D.py
Average velocity in x direction: 0.00972682241169121
mpirun -n 2 python3 FH2D_PerRB_MWE_v5.py
Average velocity in x direction: 0.009613288617357399

The wrong results only appear at the slave DOFS. Do you have any idea what is going wrong or what I am missing?

Here is the MWE (again, but for v0.7.0):

import gmsh
import dolfinx
import dolfinx_mpc
from dolfinx.fem import FunctionSpace, Function, Constant, assemble_scalar, form
from dolfinx.mesh import locate_entities_boundary, meshtags
from dolfinx.io.gmshio import model_to_mesh

from ufl import grad, div, inner, Measure, TrialFunctions, TestFunctions, FiniteElement, VectorElement
import ufl
from petsc4py import PETSc
import numpy as np
from mpi4py import MPI

# MPI initialization
comm = MPI.COMM_WORLD
rank = MPI.COMM_WORLD.rank

#############################################
#              Mesh generation              #
#############################################
gmsh.initialize()
L = 1                       # length of rectangle
f = 0.3                     # volume fraction of inclusion
r = np.sqrt(f*L**2/np.pi)   # radius of inclusion

gdim = 2         # Geometric dimension of the mesh
fdim = gdim - 1  # facet dimension

if rank == 0:
    inclusion = gmsh.model.occ.addDisk(L / 2, L / 2, 0, r, r)
    matrix = gmsh.model.occ.addRectangle(0,0,0, L, L)
    whole_domain = gmsh.model.occ.fragment([(2, matrix)], [(2, inclusion)])
    gmsh.model.occ.synchronize()

    matrix_physical = gmsh.model.addPhysicalGroup(gdim, [matrix], tag = 1)
    inclusion_physical = gmsh.model.addPhysicalGroup(gdim, [inclusion], tag = 2)

    # Get the interface elements between rectangle and circle and tag
    interface_elements = gmsh.model.getBoundary([(gdim, inclusion_physical)])
    interface_physical = gmsh.model.addPhysicalGroup(fdim, (5, 5), tag = 9)

    background_surfaces = []
    inclusion_surfaces = []
    for domain in whole_domain[0]:
        com = gmsh.model.occ.getCenterOfMass(domain[0], domain[1])
        mass = gmsh.model.occ.getMass(domain[0], domain[1])
        if np.isclose(mass, np.pi * (r ** 2)): # identify inclusion by its mass
            inclusion_surfaces.append(domain)
        else:
            background_surfaces.append(domain)

    gmsh.model.mesh.field.add("Distance", 1)
    edges = gmsh.model.getBoundary(inclusion_surfaces, oriented=False)
    gmsh.model.mesh.field.setNumbers(1, "EdgesList", [e[1] for e in edges])

    gmsh.model.mesh.field.add("Threshold", 2)
    gmsh.model.mesh.field.setNumber(2, "IField", 1)
    gmsh.model.mesh.field.setNumber(2, "LcMin", 0.01)
    gmsh.model.mesh.field.setNumber(2, "LcMax", 0.03)
    gmsh.model.mesh.field.setNumber(2, "DistMin", 0.01)
    gmsh.model.mesh.field.setNumber(2, "DistMax", 0.1)
    gmsh.model.mesh.field.setAsBackgroundMesh(2)

    gmsh.option.setNumber("Mesh.Algorithm", 6)
    gmsh.model.mesh.generate(gdim)

domain, ct, ft = model_to_mesh(gmsh.model, comm, 0, gdim=gdim)
gmsh.finalize()
print("Mesh is loaded")
#############################################
#             Mesh generation finished      #
#############################################

dx = Measure('dx', domain=domain, subdomain_data=ct)
num_cells = domain.topology.index_map(domain.topology.dim).size_global

# elements and funcionspace
###########################

# Taylor-Hood elements and Mixed Functionspace
P1 = FiniteElement("CG", domain.ufl_cell(), 1)
P2 = VectorElement("CG", domain.ufl_cell(), 2)
TH = ufl.MixedElement([P2, P1])
V = FunctionSpace(domain, TH)

# define function, trial function and test function
(ve, pr) = TrialFunctions(V)
(dve, dpr) = TestFunctions(V)

# define material parameter
S = FunctionSpace(domain, ("DG", 0))
vis = Function(S)
fluid_visc = ct.find(1)
vis.x.array[fluid_visc] = np.full_like(fluid_visc, 1, dtype=PETSc.ScalarType)
inclusion_visc = ct.find(2)
vis.x.array[inclusion_visc]  = np.full_like(inclusion_visc, 1e+5, dtype=PETSc.ScalarType)
vis.x.scatter_forward()

# boundary conditions
#####################
## set all fluid velocities at the interface to zero (no slip boundary condition)
# create collapsed functionspace
V_coll0, _ = V.sub(0).collapse()
v_bc = Function(V_coll0)
v_bc.x.array[:] = 0.0
v_bc.x.scatter_forward()
# dofs of interface and dirichlet bcs
dofs_v_interf = dolfinx.fem.locate_dofs_topological((V.sub(0), V_coll0)  , ft.dim, ft.find(9))  # interface tag: 9
bcs_v_interf = dolfinx.fem.dirichletbc(v_bc, dofs_v_interf, V.sub(0))

## set pressures at the corners to zero
def boundary_node(x):
    return np.logical_or(
        np.logical_and(
            np.isclose(x[0], 0.),
            np.logical_or(
                np.isclose(x[1], .0),
                np.isclose(x[1], L),
            )
        ),
        np.logical_and(
            np.isclose(x[0], L),
            np.logical_or(
                np.isclose(x[1], 0.),
                np.isclose(x[1], L),
            )
        ),
    )

boundary_facets = locate_entities_boundary(domain, 0, boundary_node) # it is zero dimension for a point
dofs_corners = dolfinx.fem.locate_dofs_topological(V.sub(1), 0, boundary_facets)  # it is zero dimension for a point
bc_corners = dolfinx.fem.dirichletbc(PETSc.ScalarType(0.0), dofs_corners, V.sub(1))

bcs = [bcs_v_interf, bc_corners]

print("Set periodic boundary conditions...")
## periodic boundary conditions
pbc_directions = []
pbc_slave_tags = []
pbc_is_slave = []
pbc_is_master = []
pbc_meshtags = []
pbc_slave_to_master_maps = []

mpc = dolfinx_mpc.MultiPointConstraint(V)

def generate_pbc_slave_to_master_map(i):
    def pbc_slave_to_master_map(x):
        out_x = x.copy()
        out_x[i] = x[i] - L
        return out_x
    return pbc_slave_to_master_map

def generate_pbc_is_slave(i):
    return lambda x: np.isclose(x[i], L)

def generate_pbc_is_master(i):
    return lambda x: np.isclose(x[i], 0.0)

for i in range(gdim):

    pbc_directions.append(i)
    pbc_slave_tags.append(i + 2)
    pbc_is_slave.append(generate_pbc_is_slave(i))
    pbc_is_master.append(generate_pbc_is_master(i))
    pbc_slave_to_master_maps.append(generate_pbc_slave_to_master_map(i))

    facets = locate_entities_boundary(domain, fdim, pbc_is_slave[-1])
    arg_sort = np.argsort(facets)
    pbc_meshtags.append(meshtags(domain,
                                 fdim,
                                 facets[arg_sort],
                                 np.full(len(facets), pbc_slave_tags[-1], dtype=np.int32)))

N_pbc = len(pbc_directions)
for i in range(N_pbc):
    if N_pbc > 1:
        def pbc_slave_to_master_map(x):
            out_x = pbc_slave_to_master_maps[i](x)
            idx = pbc_is_slave[(i + 1) % N_pbc](x)
            out_x[pbc_directions[i]][idx] = np.nan
            return out_x
    else:
        pbc_slave_to_master_map = pbc_slave_to_master_maps[i]

    for ij in range(gdim):
        mpc.create_periodic_constraint_topological(V.sub(0).sub(ij), pbc_meshtags[i],
                                                   pbc_slave_tags[i],
                                                   pbc_slave_to_master_map,
                                                   bcs)

    mpc.create_periodic_constraint_topological(V.sub(1), pbc_meshtags[i],
                                                   pbc_slave_tags[i],
                                                   pbc_slave_to_master_map,
                                                   bcs)

if len(pbc_directions) > 1:
    # Map intersection(slaves_x, slaves_y) to intersection(masters_x, masters_y),
    # i.e. map the slave dof at (1, 1)/(1,1,1) to the master dof at (0, 0)/(0,0,0)
    def pbc_slave_to_master_map(x):
        out_x = x.copy()
        out_x[0] = x[0] - domain.geometry.x[:, 0].max()
        out_x[1] = x[1] - domain.geometry.x[:, 1].max()
        idx = np.logical_and(pbc_is_slave[0](x), pbc_is_slave[1](x))
        out_x[0][~idx] = np.nan
        out_x[1][~idx] = np.nan
        return out_x

    for ij in range(gdim):
        mpc.create_periodic_constraint_topological(V.sub(0).sub(ij), pbc_meshtags[1],
                                                   pbc_slave_tags[1],
                                                   pbc_slave_to_master_map,
                                                   bcs)

    mpc.create_periodic_constraint_topological(V.sub(1), pbc_meshtags[1],
                                               pbc_slave_tags[1],
                                               pbc_slave_to_master_map,
                                               bcs)

mpc.finalize()

# Variational problem
#####################
f_constant = Constant(domain, PETSc.ScalarType((1,0)))

## weak formulation
a = vis * inner(grad(ve), grad(dve)) * dx + pr * div(dve) * dx + div(ve) * dpr * dx
L = inner(f_constant, dve) * dx
a_f = form(a)
L_f = form(L)

# solution function
func_sol = Function(V)

## iterative solver
A = dolfinx_mpc.assemble_matrix(a_f, mpc, bcs=bcs)
A.assemble()
b = dolfinx_mpc.assemble_vector(L_f, mpc)
dolfinx_mpc.apply_lifting(b, [a_f], [bcs], mpc)
b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
dolfinx.fem.petsc.set_bc(b, bcs)
b.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)

# Preconditioner matrix form
Pf = dolfinx.fem.form(ufl.inner(ufl.grad(ve), ufl.grad(dve)) * dx + pr * dpr * dx)
P = dolfinx_mpc.assemble_matrix(Pf, mpc, bcs=bcs)
P.assemble()

# solver settings
ksp = PETSc.KSP().create(domain.comm)
ksp.setOperators(A, P)
ksp.setType(PETSc.KSP.Type.MINRES)
pc = ksp.getPC()
pc.setType("hypre")
pc.setHYPREType("boomeramg")
pc.setUp()

# solving and backsubstitution to obtain the values at the slave degrees of freedom
ksp.solve(b, func_sol.vector)
mpc.backsubstitution(func_sol)

# collapse mixed functionspace
ve = func_sol.sub(0).collapse()
pr = func_sol.sub(1).collapse()

# average velocity
v_avg = domain.comm.allreduce(assemble_scalar(form(ve[0] * dx)), op=MPI.SUM)
if rank == 0:
    print("Average velocity in x direction:", v_avg)

# intialize files for the results
with dolfinx.io.VTXWriter(MPI.COMM_WORLD, "results_2D_iterative/pressure.bp", [pr], engine="BP4") as vtx:
    vtx.write(0.0)
with dolfinx.io.VTXWriter(MPI.COMM_WORLD, "results_2D_iterative/velocity.bp", [ve], engine="BP4") as vtx:
    vtx.write(0.0)

I observe the same behavior, however, your code is quite complex, and I don’t have time to dig into it right now.
One thing to note, is that you are not calling func_sol.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
after the ksp-sol, but it doesn’t change the code, so something else is going on as well.

If you could reduce the problem to something on a built in mesh, it would be easier for me to pinpoint what goes wrong.

Thank you for your fast respone!
I modified the code to a built-in mesh (but with additional lines for cell tags and facet tags). The results are as expected. However, running in parallel gives different results again.
Viewing the fields in Paraview, only the values at the slave surfaces are wrong. So I guess there might happen something with the backsubstitution.
Here is the code: (If you have any whishes on how to modify it so you can easily see what is wrong, please tell me.)

import dolfinx
import dolfinx_mpc
from dolfinx.fem import FunctionSpace, Function, Constant, assemble_scalar, form
from dolfinx.mesh import locate_entities_boundary, meshtags

from ufl import grad, div, inner, Measure, TrialFunctions, TestFunctions, FiniteElement, VectorElement
import ufl
from petsc4py import PETSc
import numpy as np
from mpi4py import MPI
# MPI initialization
comm = MPI.COMM_WORLD
rank = MPI.COMM_WORLD.rank

#############################################
#              Mesh generation              #
#############################################
domain = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 50, 50, dolfinx.mesh.CellType.quadrilateral)
gdim = 2
fdim = 1
##  cells for different subdomains
def cells_1(x):
    return np.logical_or(np.logical_or(x[0] >= 0.6, x[0] <= 0.4), np.logical_or(x[1] >= 0.6, x[1] <= 0.4))
def cells_2(x):
    return np.logical_and(np.logical_and(x[0] <= 0.6, x[1] <= 0.6), np.logical_and(x[0] >= 0.4, x[1] >= 0.4))

# locate entities
cells_1_dom = dolfinx.mesh.locate_entities(domain, domain.topology.dim, cells_1)
cells_2_dom = dolfinx.mesh.locate_entities(domain, domain.topology.dim, cells_2)

## create cell tags
marked_cells = np.hstack([cells_1_dom, cells_2_dom])
marked_values = np.hstack([np.full_like(cells_1_dom, 1), np.full_like(cells_2_dom, 2)])
sorted_cells = np.argsort(marked_cells)
ct = dolfinx.mesh.meshtags(domain, gdim, marked_cells[sorted_cells], marked_values[sorted_cells])

## create facet tag at interface
def create_interface(facet_marker: np.typing.NDArray[np.int32],
                     tag: int,
                     fa: np.typing.NDArray[np.int32],
                     fo: np.typing.NDArray[np.int32],
                     cell_marker: np.typing.NDArray[np.int32]):
    for f in range(len(facet_marker)):
        cells = fa[fo[f]:fo[f + 1]]
        c_val = cell_marker[cells]
        if len(np.unique(c_val)) == 2:
            facet_marker[f] = tag


# Cold start of function
create_interface(np.empty(0, dtype=np.int32), 0, np.empty(
    0, dtype=np.int32), np.empty(0, dtype=np.int32), np.empty(0, dtype=np.int32))

cell_tag_1 = np.unique(ct.values)[0]
cell_tag_2 = np.unique(ct.values)[1]

domain.topology.create_connectivity(domain.topology.dim - 1, domain.topology.dim)
f_to_c = domain.topology.connectivity(domain.topology.dim - 1, domain.topology.dim)
fa = f_to_c.array
fo = f_to_c.offsets

facet_map = domain.topology.index_map(domain.topology.dim - 1)
num_facets = facet_map.size_local + facet_map.num_ghosts
facet_marker = np.zeros(num_facets, dtype=np.int32)
cell_map = domain.topology.index_map(domain.topology.dim)
num_cells = cell_map.size_local + cell_map.num_ghosts
cell_marker = np.full(num_cells, cell_tag_2, dtype=np.int32)
cell_marker[ct.indices] = ct.values

interface_marker = 9
create_interface(facet_marker, interface_marker, fa, fo, cell_marker)
ft = meshtags(domain, domain.topology.dim - 1,
                     np.arange(num_facets), facet_marker)

dx = Measure('dx', domain=domain, subdomain_data=ct)
#############################################
#             Mesh generation finished      #
#############################################



# elements and funcionspace
###########################

# Taylor-Hood elements and Mixed Functionspace
P1 = FiniteElement("CG", domain.ufl_cell(), 1)
P2 = VectorElement("CG", domain.ufl_cell(), 2)
TH = ufl.MixedElement([P2, P1])
V = FunctionSpace(domain, TH)

# define function, trial function and test function
(ve, pr) = TrialFunctions(V)
(dve, dpr) = TestFunctions(V)

# define material parameter
S = FunctionSpace(domain, ("DG", 0))
vis = Function(S)
fluid_visc = ct.find(1)
vis.x.array[fluid_visc] = np.full_like(fluid_visc, 1, dtype=PETSc.ScalarType)
inclusion_visc = ct.find(2)
vis.x.array[inclusion_visc]  = np.full_like(inclusion_visc, 1e+5, dtype=PETSc.ScalarType)
vis.x.scatter_forward()

# boundary conditions
#####################
## set all fluid velocities at the interface to zero (no slip boundary condition)
# create collapsed functionspace
V_coll0, _ = V.sub(0).collapse()
v_bc = Function(V_coll0)
v_bc.x.array[:] = 0.0
v_bc.x.scatter_forward()
# dofs of interface and dirichlet bcs
dofs_v_interf = dolfinx.fem.locate_dofs_topological((V.sub(0), V_coll0)  , ft.dim, ft.find(9))  # interface tag: 9
bcs_v_interf = dolfinx.fem.dirichletbc(v_bc, dofs_v_interf, V.sub(0))

## set pressures at the corners to zero
def boundary_node(x):
    return np.logical_or(
        np.logical_and(
            np.isclose(x[0], 0.),
            np.logical_or(
                np.isclose(x[1], .0),
                np.isclose(x[1], 1.0),
            )
        ),
        np.logical_and(
            np.isclose(x[0], 1.0),
            np.logical_or(
                np.isclose(x[1], 0.),
                np.isclose(x[1], 1.0),
            )
        ),
    )

boundary_facets = locate_entities_boundary(domain, 0, boundary_node) # it is zero dimension for a point
dofs_corners = dolfinx.fem.locate_dofs_topological(V.sub(1), 0, boundary_facets)  # it is zero dimension for a point
bc_corners = dolfinx.fem.dirichletbc(PETSc.ScalarType(0.0), dofs_corners, V.sub(1))

bcs = [bcs_v_interf, bc_corners]

## periodic boundary conditions
pbc_directions = []
pbc_slave_tags = []
pbc_is_slave = []
pbc_is_master = []
pbc_meshtags = []
pbc_slave_to_master_maps = []

mpc = dolfinx_mpc.MultiPointConstraint(V)

def generate_pbc_slave_to_master_map(i):
    def pbc_slave_to_master_map(x):
        out_x = x.copy()
        out_x[i] = x[i] - 1.0
        return out_x
    return pbc_slave_to_master_map

def generate_pbc_is_slave(i):
    return lambda x: np.isclose(x[i], 1.0)

def generate_pbc_is_master(i):
    return lambda x: np.isclose(x[i], 0.0)

for i in range(gdim):

    pbc_directions.append(i)
    pbc_slave_tags.append(i + 2)
    pbc_is_slave.append(generate_pbc_is_slave(i))
    pbc_is_master.append(generate_pbc_is_master(i))
    pbc_slave_to_master_maps.append(generate_pbc_slave_to_master_map(i))

    facets = locate_entities_boundary(domain, fdim, pbc_is_slave[-1])
    arg_sort = np.argsort(facets)
    pbc_meshtags.append(meshtags(domain,
                                 fdim,
                                 facets[arg_sort],
                                 np.full(len(facets), pbc_slave_tags[-1], dtype=np.int32)))

N_pbc = len(pbc_directions)
for i in range(N_pbc):
    if N_pbc > 1:
        def pbc_slave_to_master_map(x):
            out_x = pbc_slave_to_master_maps[i](x)
            idx = pbc_is_slave[(i + 1) % N_pbc](x)
            out_x[pbc_directions[i]][idx] = np.nan
            return out_x
    else:
        pbc_slave_to_master_map = pbc_slave_to_master_maps[i]

    for ij in range(gdim):
        mpc.create_periodic_constraint_topological(V.sub(0).sub(ij), pbc_meshtags[i],
                                                   pbc_slave_tags[i],
                                                   pbc_slave_to_master_map,
                                                   bcs)

    mpc.create_periodic_constraint_topological(V.sub(1), pbc_meshtags[i],
                                                   pbc_slave_tags[i],
                                                   pbc_slave_to_master_map,
                                                   bcs)

if len(pbc_directions) > 1:
    def pbc_slave_to_master_map(x):
        out_x = x.copy()
        out_x[0] = x[0] - domain.geometry.x[:, 0].max()
        out_x[1] = x[1] - domain.geometry.x[:, 1].max()
        idx = np.logical_and(pbc_is_slave[0](x), pbc_is_slave[1](x))
        out_x[0][~idx] = np.nan
        out_x[1][~idx] = np.nan
        return out_x

    for ij in range(gdim):
        mpc.create_periodic_constraint_topological(V.sub(0).sub(ij), pbc_meshtags[1],
                                                   pbc_slave_tags[1],
                                                   pbc_slave_to_master_map,
                                                   bcs)

    mpc.create_periodic_constraint_topological(V.sub(1), pbc_meshtags[1],
                                               pbc_slave_tags[1],
                                               pbc_slave_to_master_map,
                                               bcs)

mpc.finalize()

# Variational problem
#####################
f_constant = Constant(domain, PETSc.ScalarType((1,0)))

## weak formulation
a_f = form(vis * inner(grad(ve), grad(dve)) * dx + pr * div(dve) * dx + div(ve) * dpr * dx)
L_f = form(inner(f_constant, dve) * dx)

# solution function
func_sol = Function(V)

## iterative solver
A = dolfinx_mpc.assemble_matrix(a_f, mpc, bcs=bcs)
A.assemble()
b = dolfinx_mpc.assemble_vector(L_f, mpc)
dolfinx_mpc.apply_lifting(b, [a_f], [bcs], mpc)
b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
dolfinx.fem.petsc.set_bc(b, bcs)
b.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)

# Preconditioner matrix form
Pf = dolfinx.fem.form(ufl.inner(ufl.grad(ve), ufl.grad(dve)) * dx + pr * dpr * dx)
P = dolfinx_mpc.assemble_matrix(Pf, mpc, bcs=bcs)
P.assemble()

# solver settings
ksp = PETSc.KSP().create(domain.comm)
ksp.setOperators(A, P)
ksp.setType(PETSc.KSP.Type.MINRES)
pc = ksp.getPC()
pc.setType("hypre")
pc.setHYPREType("boomeramg")
pc.setUp()

# solving and backsubstitution to obtain the values at the slave degrees of freedom
ksp.solve(b, func_sol.vector)
func_sol.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
mpc.backsubstitution(func_sol)

# collapse mixed functionspace
ve = func_sol.sub(0).collapse()
pr = func_sol.sub(1).collapse()

# average velocity
v_avg = domain.comm.allreduce(assemble_scalar(form(ve[0] * dx)), op=MPI.SUM)
if rank == 0:
    print("Average velocity in x direction:", v_avg)

# write results to file
with dolfinx.io.VTXWriter(MPI.COMM_WORLD, "results_2D_square/pressure.bp", [pr], engine="BP4") as vtx:
    vtx.write(0.0)
with dolfinx.io.VTXWriter(MPI.COMM_WORLD, "results_2D_square/velocity.bp", [ve], engine="BP4") as vtx:
    vtx.write(0.0)

I still do not know how to handle the error of the backsubstitution that occurs running the code in parallel. I tested another physical problem that is easier to overview (no Preconditioning etc.). However, it still does not work. If we look at the figure, at some places, the back substitution works. No back substituion happens at the places indicated with the red rectangles - the values of the DOFs are still 0 (for mpirun -n 2).

Do you think that when mpc.backsubstitution is called, only the slave DOFs of one certain rank are obtained? Do you have an idea how to obtain the DOFs of all slave nodes?

Here is the MWE:

import gmsh
import dolfinx
import dolfinx_mpc
from dolfinx.fem import FunctionSpace, Function, Constant, assemble_scalar, form
from dolfinx.mesh import locate_entities_boundary, meshtags
from dolfinx.io.gmshio import model_to_mesh

from ufl import grad, Measure, TrialFunction, TestFunction, FiniteElement
import ufl
from petsc4py import PETSc
import numpy as np
from mpi4py import MPI

# MPI initialization
comm = MPI.COMM_WORLD
rank = MPI.COMM_WORLD.rank

#############################################
#              Mesh generation              #
#############################################
gmsh.initialize()
L = 1                       # length of rectangle
f = 0.3                     # volume fraction of inclusion
r = np.sqrt(f*L**2/np.pi)   # radius of inclusion

gdim = 2         # Geometric dimension of the mesh
fdim = gdim - 1  # facet dimension

if rank == 0:
    inclusion = gmsh.model.occ.addDisk(L / 2, L / 2, 0, r, r)
    matrix = gmsh.model.occ.addRectangle(0,0,0, L, L)
    whole_domain = gmsh.model.occ.fragment([(2, matrix)], [(2, inclusion)])
    gmsh.model.occ.synchronize()

    matrix_physical = gmsh.model.addPhysicalGroup(gdim, [matrix], tag = 1)
    inclusion_physical = gmsh.model.addPhysicalGroup(gdim, [inclusion], tag = 2)

    # Get the interface elements between rectangle and circle and tag
    interface_elements = gmsh.model.getBoundary([(gdim, inclusion_physical)])
    interface_physical = gmsh.model.addPhysicalGroup(fdim, (5, 5), tag = 9)

    background_surfaces = []
    inclusion_surfaces = []
    for domain in whole_domain[0]:
        com = gmsh.model.occ.getCenterOfMass(domain[0], domain[1])
        mass = gmsh.model.occ.getMass(domain[0], domain[1])
        if np.isclose(mass, np.pi * (r ** 2)): # identify inclusion by its mass
            inclusion_surfaces.append(domain)
        else:
            background_surfaces.append(domain)

    gmsh.model.mesh.field.add("Distance", 1)
    edges = gmsh.model.getBoundary(inclusion_surfaces, oriented=False)
    gmsh.model.mesh.field.setNumbers(1, "EdgesList", [e[1] for e in edges])

    gmsh.model.mesh.field.add("Threshold", 2)
    gmsh.model.mesh.field.setNumber(2, "IField", 1)
    gmsh.model.mesh.field.setNumber(2, "LcMin", 0.01)
    gmsh.model.mesh.field.setNumber(2, "LcMax", 0.03)
    gmsh.model.mesh.field.setNumber(2, "DistMin", 0.01)
    gmsh.model.mesh.field.setNumber(2, "DistMax", 0.1)
    gmsh.model.mesh.field.setAsBackgroundMesh(2)

    gmsh.option.setNumber("Mesh.Algorithm", 6)
    gmsh.model.mesh.generate(gdim)

domain, ct, ft = model_to_mesh(gmsh.model, comm, 0, gdim=gdim)
gmsh.finalize()
print("Mesh is loaded")
#############################################
#             Mesh generation finished      #
#############################################

dx = Measure('dx', domain=domain, subdomain_data=ct)
x = ufl.SpatialCoordinate(domain)

# elements and funcionspace
###########################

# Taylor-Hood elements and Mixed Functionspace
P1 = FiniteElement("CG", domain.ufl_cell(), 1)
V = FunctionSpace(domain, P1)

# define function, trial function and test function
T_fluc = TrialFunction(V)
dT_fluc = TestFunction(V)

# define material parameter
S = FunctionSpace(domain, ("DG", 0))
kappa = Function(S)
fluid_kappa = ct.find(1)
kappa.x.array[fluid_kappa] = np.full_like(fluid_kappa, 1, dtype=PETSc.ScalarType)
inclusion_kappa = ct.find(2)
kappa.x.array[inclusion_kappa]  = np.full_like(inclusion_kappa, 1e+5, dtype=PETSc.ScalarType)
kappa.x.scatter_forward()

# boundary conditions
#####################
bcs = []

print("Set periodic boundary conditions...")
## periodic boundary conditions
pbc_directions = []
pbc_slave_tags = []
pbc_is_slave = []
pbc_is_master = []
pbc_meshtags = []
pbc_slave_to_master_maps = []

mpc = dolfinx_mpc.MultiPointConstraint(V)

def generate_pbc_slave_to_master_map(i):
    def pbc_slave_to_master_map(x):
        out_x = x.copy()
        out_x[i] = x[i] - L
        return out_x
    return pbc_slave_to_master_map

def generate_pbc_is_slave(i):
    return lambda x: np.isclose(x[i], L)

def generate_pbc_is_master(i):
    return lambda x: np.isclose(x[i], 0.0)

for i in range(gdim):

    pbc_directions.append(i)
    pbc_slave_tags.append(i + 2)
    pbc_is_slave.append(generate_pbc_is_slave(i))
    pbc_is_master.append(generate_pbc_is_master(i))
    pbc_slave_to_master_maps.append(generate_pbc_slave_to_master_map(i))

    facets = locate_entities_boundary(domain, fdim, pbc_is_slave[-1])
    arg_sort = np.argsort(facets)
    pbc_meshtags.append(meshtags(domain,
                                 fdim,
                                 facets[arg_sort],
                                 np.full(len(facets), pbc_slave_tags[-1], dtype=np.int32)))

N_pbc = len(pbc_directions)
for i in range(N_pbc):
    if N_pbc > 1:
        def pbc_slave_to_master_map(x):
            out_x = pbc_slave_to_master_maps[i](x)
            idx = pbc_is_slave[(i + 1) % N_pbc](x)
            out_x[pbc_directions[i]][idx] = np.nan
            return out_x
    else:
        pbc_slave_to_master_map = pbc_slave_to_master_maps[i]

    mpc.create_periodic_constraint_topological(V, pbc_meshtags[i],
                                                   pbc_slave_tags[i],
                                                   pbc_slave_to_master_map,
                                                   bcs)

if len(pbc_directions) > 1:
    # Map intersection(slaves_x, slaves_y) to intersection(masters_x, masters_y),
    # i.e. map the slave dof at (1, 1)/(1,1,1) to the master dof at (0, 0)/(0,0,0)
    def pbc_slave_to_master_map(x):
        out_x = x.copy()
        out_x[0] = x[0] - domain.geometry.x[:, 0].max()
        out_x[1] = x[1] - domain.geometry.x[:, 1].max()
        idx = np.logical_and(pbc_is_slave[0](x), pbc_is_slave[1](x))
        out_x[0][~idx] = np.nan
        out_x[1][~idx] = np.nan
        return out_x

    mpc.create_periodic_constraint_topological(V, pbc_meshtags[1],
                                               pbc_slave_tags[1],
                                               pbc_slave_to_master_map,
                                               bcs)

mpc.finalize()

# Variational problem
#####################
T_gradient = Constant(domain, PETSc.ScalarType((1,0)))

## weak formulation
eqn = - kappa * ufl.dot(grad(ufl.dot(T_gradient , x) + T_fluc), grad(dT_fluc)) * dx
a_f = form(ufl.lhs(eqn))
L_f = form(ufl.rhs(eqn))

# solution function
func_sol = Function(V)

## iterative solver
A = dolfinx_mpc.assemble_matrix(a_f, mpc, bcs=bcs)
A.assemble()
b = dolfinx_mpc.assemble_vector(L_f, mpc)
dolfinx_mpc.apply_lifting(b, [a_f], [bcs], mpc)
b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
dolfinx.fem.petsc.set_bc(b, bcs)
b.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)

# solver settings
ksp = PETSc.KSP().create(domain.comm)
ksp.setOperators(A)
ksp.setType(PETSc.KSP.Type.PREONLY)
pc = ksp.getPC()
pc.setType("lu")
pc.setHYPREType("mumps")
pc.setUp()

# solving and backsubstitution to obtain the values at the slave degrees of freedom
ksp.solve(b, func_sol.vector)
func_sol.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
mpc.backsubstitution(func_sol)

# average quantity
T_fluc_avg = domain.comm.allreduce(assemble_scalar(form(func_sol * dx)), op=MPI.SUM)
if rank == 0:
    print(f"Average quantity:", T_fluc_avg)

# intialize files for the results
with dolfinx.io.VTXWriter(MPI.COMM_WORLD, f"results_2D_mwe/Quantity.bp", [func_sol], engine="BP4") as vtx:
    vtx.write(0.0)

I have quite alot of other commitments on my plate at the time.
As I’ve got several working implementations of periodic boundary conditions as demos, I’m quite suprised that this does not work.
This is on my list of priorities to debug in December.

Thank you very much for your efforts! I would be grateful if you would let me know if you find an error.

Hi, I just wanted to say that I independently stumbled over the exact same problem described here when trying to run a periodic 3D stokes flow case. For me as well, mpc.backsubstitution does not produce the expected outcome in parallel, whereas on a single processor it works just fine. I am very interested in any updates anyone might have on this topic, thank you so much!

Hi again, just a quick update on this, in case it might be helpful:

  1. I have verified that mpc.masters uses local indexing, as opposed to the documentation (dolfinx_mpc — DOLFINx-MPC: An extension to DOLFINx for multi point constraints) which states that these are global indices.

  2. The following piece of code works for me in parallel as a momentary replacement for mpc.backsubstitution(wh):

from bisect import bisect_right
mpc.homogenize(wh)
l2g = mpc.function_space.dofmap.index_map.local_to_global
n_local_dofs = mpc.function_space.dofmap.index_map.size_local
n_all_dofs = np.max(mpc.masters.array)+1 if mpc.masters.array.size > 0 else n_local_dofs
nonlocal_masters = l2g(np.arange(n_local_dofs, n_all_dofs, dtype=np.int32))
dof_offsets = np.cumsum(mesh.comm.allgather(n_local_dofs))
def global2local_dofs(idglobal):
    pr = bisect_right(dof_offsets, idglobal)
    if not pr:
        idlocal = idglobal
    else:
        idlocal = idglobal - dof_offsets[pr-1]
    return pr, idlocal
request_out = [[] for _ in range(mesh.comm.size)]
sorting = [[] for _ in range(mesh.comm.size)]
for i in range(nonlocal_masters.size):
    ownerproc, ownerdofid = global2local_dofs(nonlocal_masters[i])
    request_out[ownerproc].append(ownerdofid)
    sorting[ownerproc].append(i)
request_in = [mesh.comm.scatter(request_out, root=i) for i in range(mesh.comm.size)]
answer_out = [wh.x.array[request_in[i]] for i in range(mesh.comm.size)]
answer_in = [mesh.comm.scatter(answer_out, root=i) for i in range(mesh.comm.size)]
nonlocal_wh = np.hstack(answer_in)[np.argsort(np.hstack(sorting))]
coeffs, offsets = mpc.coefficients()
for i in range(mpc.num_local_slaves):
    islv = mpc.slaves[i]
    coeffs_of_slave_i = coeffs[offsets[islv]:offsets[islv+1]]
    masters_of_slave_i = mpc.masters.array[offsets[islv]:offsets[islv+1]]
    for j in range(masters_of_slave_i.size):
        if masters_of_slave_i[j] < n_local_dofs:
            wh.x.array[islv] += coeffs_of_slave_i[j] * wh.x.array[masters_of_slave_i[j]]
        else:
            wh.x.array[islv] += coeffs_of_slave_i[j] * nonlocal_wh[masters_of_slave_i[j] - n_local_dofs]
wh.x.scatter_forward()

1 Like

The input to the MPC uses global indices, as the dofs might be from a process that originally might not have a local index.

Internally, dolfinx_mpc creates new extended index maps, and adds the global dofs as a local dof on the relevant process.


Right, I got confused, it was the input…

Why are you creating a map with output equal to nan?

Do you have a code that illustrates your issue? It is quite hard to help debugging this, when I got plenty of passing backsubstitution tests/demos, and the example above is not complete, or uses np.nan in mappings, that to me doesn’t make sense

I do, but the mesh is not fully “built-in”. I will soon adapt it and share it so that the problem I experience can be easily reproduced. Thank you very much.

Just a note that in my case I do not use np.nan in mappings (with which I think tuderic seeks to avoid duplicate constraints, but I do it differently) and I still get that issue with parallel backsubstitution (but not if I use the aforementioned “workaround”).

Let me know when you have something I can look at.

The output should not be equal to nan. It should contain the DOFs at the corner that are not mapped before. Only the DOFs that are not at this corner are set to nan (out_x[0][~idx]). This is the same procedure as in the dolfinx_mpc demo.