Real-space Lagrange multiplier constraint silently violated (non-deterministically) when a periodic MPC slave/master pair spans different MPI ranks — reproducible on v0.10.5 and v0.11.0

dolfinx
,
1

I’m solving a periodic elasticity/Poisson-type problem where the periodic fluctuation field v (in a MultiPointConstraint-constrained space V) has its mean pinned to zero via a scalar Lagrange multiplier gamma living in a “Real” space R (scifem.create_real_functionspace). I expect that

sum(∫v dx) ≈ 1e-17 (machine precision) for any number of MPI ranks, since this is exactly what the Lagrange multiplier system is solving for. This is what happens:

  • -n 1: always correct (~1e-19).

  • -n 2 / -n 3 / -n 4: frequently wrong, but not consistently — the same script, same parameters, same rank count occasionally gives the correct machine-precision result and occasionally doesn’t, run to run.
    Here a minimal reproducer script:

    import numpy as np

import ufl

import dolfinx

import dolfinx.fem.petsc

from dolfinx.mesh import create_rectangle, CellType, GhostMode

from mpi4py import MPI

from petsc4py import PETSc

from scifem import create_real_functionspace

from dolfinx_mpc import (

    MultiPointConstraint,

    create_matrix_nest,

    assemble_matrix_nest,

    create_vector_nest,

    assemble_vector_nest,

)

comm = MPI.COMM_WORLD

Lx, Ly = 1.0, 0.3

nx, ny = 40, 12

mesh = create_rectangle(

    comm,

    [[0.0, 0.0], [Lx, Ly]],

    [nx, ny],

    cell_type=CellType.triangle,

    ghost_mode=GhostMode.shared_facet,

)

V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1))

R = create_real_functionspace(mesh)

vh = dolfinx.fem.Function(V, name="v")

bcs = []

def is_right(x):

    return np.isclose(x[0], Lx)

def is_top(x):

    return np.isclose(x[1], Ly)

mpc = MultiPointConstraint(V)

def slave_to_master_map_1(x):

    out_x = x.copy()

    out_x[1] = x[1] - Ly

    idx = is_right(x)

    out_x[1][idx] = np.nan

    return out_x

mpc.create_periodic_constraint_geometrical(V, is_top, slave_to_master_map_1, bcs)

def slave_to_master_map_0(x):

    out_x = x.copy()

    out_x[0] = x[0] - Lx

    idx = is_top(x)

    out_x[0][idx] = np.nan

    return out_x

mpc.create_periodic_constraint_geometrical(V, is_right, slave_to_master_map_0, bcs)

def slave_to_master_map_corner(x):

    out_x = x.copy()

    out_x[0] = x[0] - Lx

    out_x[1] = x[1] - Ly

    idx = np.logical_and(is_right(x), is_top(x))

    out_x[0][~idx] = np.nan

    out_x[1][~idx] = np.nan

    return out_x

mpc.create_periodic_constraint_geometrical(V, is_top, slave_to_master_map_corner, bcs)

mpc.finalize()

mpc_r = MultiPointConstraint(R)

mpc_r.finalize()

v = ufl.TrialFunction(V)

gamma = ufl.TrialFunction(R)

v_star = ufl.TestFunction(V)

dl = ufl.TestFunction(R)

Eps = dolfinx.fem.Constant(mesh, dolfinx.default_scalar_type(1.0))

zero = dolfinx.fem.Constant(mesh, dolfinx.default_scalar_type(0.0))

a00 = ufl.inner(ufl.grad(v), ufl.grad(v_star)) * ufl.dx

a01 = ufl.inner(gamma, v_star) * ufl.dx

a10 = ufl.inner(v, dl) * ufl.dx

a11 = None

a = [

    [dolfinx.fem.form(a00), dolfinx.fem.form(a01)],

    [dolfinx.fem.form(a10), dolfinx.fem.form(a11)],

]

L0 = -ufl.inner(Eps, v_star) * ufl.dx

L1 = ufl.inner(zero, dl) * ufl.dx

L = [dolfinx.fem.form(L0), dolfinx.fem.form(L1)]

A = create_matrix_nest(a, [mpc, mpc_r])

assemble_matrix_nest(A, a, [mpc, mpc_r], bcs)

A.assemble()

b = create_vector_nest(L, [mpc, mpc_r])

assemble_vector_nest(b, L, [mpc, mpc_r])

for bi in b.getNestSubVecs():

    bi.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)

bcs0 = dolfinx.fem.bcs_by_block(dolfinx.fem.extract_function_spaces(L), bcs)

dolfinx.fem.petsc.set_bc(b, bcs0)

ksp = PETSc.KSP().create(mesh.comm)

ksp.setOperators(A)

ksp.setType("preonly")

pc = ksp.getPC()

pc.setType("lu")

pc.setFactorSolverType("mumps")

ksp.setFromOptions()

Xi = b.copy()

ksp.solve(b, Xi)

Xi0 = Xi.getNestSubVecs()[0]

Xi0.copy(vh.x.petsc_vec)

vh.x.scatter_forward()

mpc.backsubstitution(vh)

vh.x.scatter_forward()

mean_local = dolfinx.fem.assemble_scalar(dolfinx.fem.form(vh * ufl.dx))

mean_global = comm.allreduce(mean_local, op=MPI.SUM)

if comm.rank == 0:

    print(f"MPI size={comm.size}  mean(v)_global={mean_global:.6e}")
2

Could you try the same thing with the basix.ufl.real_element which was introduced in 0.11?
I’ve refactored that a few times now to ensure that it is correct:

3

I modified the minimal example by simply replacing R with R = dolfinx.fem.functionspace(mesh, basix.ufl.real_element(mesh.basix_cell(), ())), but I am still experiencing the exact same issue.

4

Weird. I’ll have a look.

5

@czolesi1 there are a few things to improve in your code.
Instead of using a mapping for each side, you can map them all at once with:


atol = 1e-12


def is_right(x):

    return np.isclose(x[0], Lx, atol=atol)


def is_top(x):

    return np.isclose(x[1], Ly, atol=atol)


mpc = MultiPointConstraint(V)


def slave_boundary(x):
    return np.isclose(x[0], Lx, atol=atol) | np.isclose(x[1], Ly, atol=atol)


# mapping from slave -> master
def periodic_relation(x):
    out = np.zeros_like(x)
    out[0] = x[0] - np.isclose(x[0], Lx, atol=atol)
    out[1] = x[1] - np.isclose(x[1], Ly, atol=atol)
    return out


mpc.create_periodic_constraint_geometrical(V, slave_boundary, periodic_relation, [])

mpc.finalize()

furthermore.
The main error lies in the definition of:

which in turn is used in:

In parallel, you require additional information about ghosts that the MPC add in. Therefore, vh should be created as:

vh = dolfinx.fem.Function(mpc.function_space, name="vh")

after the finalization of mpc.
Full working script:

import numpy as np

import ufl

import dolfinx

import dolfinx.fem.petsc

from dolfinx.mesh import create_rectangle, CellType, GhostMode

from mpi4py import MPI

from petsc4py import PETSc


from dolfinx_mpc import (
    MultiPointConstraint,
    create_matrix_nest,
    assemble_matrix_nest,
    create_vector_nest,
    assemble_vector_nest,
)

comm = MPI.COMM_WORLD

Lx, Ly = 1.0, 0.3

nx, ny = 40, 12

mesh = create_rectangle(
    comm,
    [[0.0, 0.0], [Lx, Ly]],
    [nx, ny],
    cell_type=CellType.triangle,
    ghost_mode=GhostMode.shared_facet,
)

V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1))
import basix.ufl

r_el = basix.ufl.real_element(mesh.basix_cell(), value_shape=())
R = dolfinx.fem.functionspace(mesh, r_el)


bcs = []

atol = 1e-12


def is_right(x):

    return np.isclose(x[0], Lx, atol=atol)


def is_top(x):

    return np.isclose(x[1], Ly, atol=atol)


mpc = MultiPointConstraint(V)


def slave_boundary(x):
    return np.isclose(x[0], Lx, atol=atol) | np.isclose(x[1], Ly, atol=atol)


# mapping from slave -> master
def periodic_relation(x):
    out = np.zeros_like(x)
    out[0] = x[0] - np.isclose(x[0], Lx, atol=atol)
    out[1] = x[1] - np.isclose(x[1], Ly, atol=atol)
    return out


mpc.create_periodic_constraint_geometrical(V, slave_boundary, periodic_relation, [])

mpc.finalize()

mpc_r = MultiPointConstraint(R)

mpc_r.finalize()

vh = dolfinx.fem.Function(mpc.function_space, name="vh")
v = ufl.TrialFunction(V)

gamma = ufl.TrialFunction(R)

v_star = ufl.TestFunction(V)

dl = ufl.TestFunction(R)

Eps = dolfinx.fem.Constant(mesh, dolfinx.default_scalar_type(1.0))

zero = dolfinx.fem.Constant(mesh, dolfinx.default_scalar_type(0.0))

a00 = ufl.inner(ufl.grad(v), ufl.grad(v_star)) * ufl.dx

a01 = ufl.inner(gamma, v_star) * ufl.dx

a10 = ufl.inner(v, dl) * ufl.dx

a11 = None

a = [
    [dolfinx.fem.form(a00), dolfinx.fem.form(a01)],
    [dolfinx.fem.form(a10), dolfinx.fem.form(a11)],
]

L0 = -ufl.inner(Eps, v_star) * ufl.dx

L1 = ufl.inner(zero, dl) * ufl.dx

L = [dolfinx.fem.form(L0), dolfinx.fem.form(L1)]

A = create_matrix_nest(a, [mpc, mpc_r])

assemble_matrix_nest(A, a, [mpc, mpc_r], bcs)

A.assemble()

b = create_vector_nest(L, [mpc, mpc_r])

assemble_vector_nest(b, L, [mpc, mpc_r])

for bi in b.getNestSubVecs():
    bi.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)

bcs0 = dolfinx.fem.bcs_by_block(dolfinx.fem.extract_function_spaces(L), bcs)

dolfinx.fem.petsc.set_bc(b, bcs0)
for bi in b.getNestSubVecs():
    bi.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
ksp = PETSc.KSP().create(mesh.comm)

ksp.setOperators(A)

ksp.setType("preonly")

pc = ksp.getPC()

pc.setType("lu")

pc.setFactorSolverType("mumps")
ksp.setErrorIfNotConverged(True)
ksp.setFromOptions()

Xi = b.copy()

ksp.solve(b, Xi)
Xi0 = Xi.getNestSubVecs()[0]

Xi0.copy(vh.x.petsc_vec)

vh.x.scatter_forward()
mpc.backsubstitution(vh)

vh.x.scatter_forward()

mean_local = dolfinx.fem.assemble_scalar(dolfinx.fem.form(vh * ufl.dx))

mean_global = comm.allreduce(mean_local, op=MPI.SUM)

if comm.rank == 0:
    print(f"MPI size={comm.size}  mean(v)_global={mean_global:.6e}")
6

Thanks a lot! It works now