`dolfinx_mpc.NonlinearProblem` with periodic geometrical constraints diverges for mixed hyperelasticity

[bhaveshshrimali]

Hi Jorgen, many thanks for the fix in #237. I was trying to implement a homogenization example which fails on the nonlinear solve and wanted to check what could be going wrong in the problem setup (since the same material model with normal dirichlet BC's converges just fine).

Summary

dolfinx_mpc.NonlinearProblem with create_periodic_constraint_geometrical on a
mixed displacement–pressure space diverges with SNES reason = -6
(DIVERGED_LINE_SEARCH) at every load step. The identical problem with standard
Dirichlet BCs converges in 2–3 Newton iterations per step.

---

Problem description

Kinematics

Unit cube $\Omega = [0,1]^3$, periodic homogenization. Total deformation gradient:

$$\mathbf{F} = \mathbf{I} + \overline{\mathbf{H}} + \nabla\tilde{\mathbf{u}},$$

where $\overline{\mathbf{H}}$ is a prescribed macro gradient (a fem.Constant) and
$\tilde{\mathbf{u}}$ is the periodic fluctuation.

Constitutive model

Near-incompressible Neo-Hookean,
mixed space $[\mathbb{P}_2]^3 \times \mathbb{P}_1$:

$$\psi = p,(J - J^\star) + \tfrac{\mu}{2}(I_1 - 3) - \mu\ln J^\star + \tfrac{\kappa}{2}(J^\star-1)^2,$$

$$I_1 = \text{tr}(\mathbf{F}^T\mathbf{F}), \quad J = \det\mathbf{F}, \quad J^\star = \frac{(\kappa+p) + \sqrt{(\kappa+p)^2 + 4\kappa\mu}}{2\kappa}.$$

Parameters: $\mu = 1$, $\kappa = 200$.

Periodic constraint Abaqus *EQUATION style

For every opposite-face node pair $(\mathbf{x}^+, \mathbf{x}^-)$ with
$\mathbf{x}^+ = \mathbf{x}^- + \mathbf{e}_k$:

$$\tilde{u}_i(\mathbf{x}^+) - \tilde{u}_i(\mathbf{x}^-) = 0.$$

In Abaqus CTENODS notation (with corner reference nodes $C_0, C_k$):

$$u_i^{\text{slave}} - u_i^{\text{master}} - u_i^{C_k} + u_i^{C_0} = 0,$$

*EQUATION
4
slave_node,  i,  1.0
master_node, i, -1.0
Ck_node,     i, -1.0
C0_node,     i,  1.0

Here $\overline{\mathbf{H}}$ enters via fem.Constant directly in UFL, so no
macroscopic strain is carried by the corner nodes; the constraint reduces to
slave = master (zero RHS). The origin is pinned for rigid-body modes; pressure
is pinned at the origin.

Periodic face-pairing (axis 0):

  y                              y
  ^  slave (x=1)                 ^  master (x=0)
  |  ●──●──●                     |  ●──●──●
  |  │  │  │  ← ũᵢˢ = ũᵢᵐ →    |  │  │  │
  |  ●──●──●                     |  ●──●──●
  └─────────→ z                  └─────────→ z

dolfinx_mpc implementation

for comp in range(3):
    for axis in range(3):
        mpc.create_periodic_constraint_geometrical(
            W.sub(0).sub(comp),        # displacement component
            _periodic_indicator(axis), # slave: high face of axis
            _periodic_relation(axis),  # mapping: x[axis] -= 1
            [bc_origin],
            1.0,
        )

---

Reproducer

issue_mpc_periodic_nonlinear.py

from __future__ import annotations
import importlib.metadata
import basix.ufl
import dolfinx_mpc
import numpy as np
import ufl
from dolfinx import default_scalar_type, fem, mesh
from dolfinx.fem import petsc as fem_petsc
from mpi4py import MPI
from petsc4py import PETSc

EXX    = 0.03   # applied macro strain
NCELLS = 3      # uniform NCELLS³ mesh
NSTEPS = 6      # load steps
MU     = 1.0    # shear modulus
BULK   = 200.0  # bulk modulus  (bulk/mu = 200 → near-incompressible)
TOL    = 1e-10

_SNES_OPTS = {
    "snes_type": "newtonls",
    "snes_linesearch_type": "bt",
    "snes_max_it": 30,
    "snes_rtol": 1e-8,
    "snes_atol": 1e-8,
    "ksp_type": "preonly",
    "pc_type": "lu",
}

def _face(axis, value):
    return lambda x: np.isclose(x[axis], value, atol=TOL)

def _origin(x):
    return np.isclose(x[0], 0.) & np.isclose(x[1], 0.) & np.isclose(x[2], 0.)

def _periodic_indicator(axis):
    def f(x):
        mask = np.isclose(x[axis], 1., atol=TOL)
        for later in range(axis + 1, 3):
            mask &= x[later] < 1. - TOL
        return mask
    return f

def _periodic_relation(axis):
    def f(x):
        y = x.copy(); y[axis] -= 1.; return y
    return f

def _zero_fn(W_sub):
    V, _ = W_sub.collapse()
    fn = fem.Function(V); fn.x.array[:] = 0.; fn.x.scatter_forward()
    return V, fn

def _build_forms(msh, W, up, macro_H):
    mu    = fem.Constant(msh, default_scalar_type(MU))
    kappa = fem.Constant(msh, default_scalar_type(BULK))
    du, v = ufl.TrialFunction(W), ufl.TestFunction(W)
    u, p  = ufl.split(up)
    I      = ufl.Identity(3)
    F      = ufl.variable(I + macro_H + ufl.grad(u))
    I1     = ufl.tr(F.T * F)
    J      = ufl.det(F)
    J_star = ((kappa + p) + ufl.sqrt((kappa + p)**2 + 4.*kappa*mu)) / (2.*kappa)
    psi    = (p*(J - J_star) + .5*mu*(I1 - 3.)
              - mu*ufl.ln(J_star) + .5*kappa*(J_star - 1.)**2)
    R = ufl.derivative(psi * ufl.dx, up, v)
    K = ufl.derivative(R, up, du)
    return R, K

# ── Case 1: Dirichlet BCs (converges) ────────────────────────────────────────
def run_nonperiodic():
    comm = MPI.COMM_WORLD
    msh  = mesh.create_unit_cube(comm, NCELLS, NCELLS, NCELLS)
    u_el = basix.ufl.element(basix.ElementFamily.P, msh.basix_cell(), 2,
                              shape=(3,), dtype=default_scalar_type)
    p_el = basix.ufl.element(basix.ElementFamily.P, msh.basix_cell(), 1,
                              dtype=default_scalar_type)
    W = fem.functionspace(msh, basix.ufl.mixed_element([u_el, p_el]))
    V_ux, ux0 = _zero_fn(W.sub(0).sub(0))
    V_uy, uy0 = _zero_fn(W.sub(0).sub(1))
    V_uz, uz0 = _zero_fn(W.sub(0).sub(2))
    Q_p,  p0  = _zero_fn(W.sub(1))
    ux_load   = fem.Function(V_ux)
    bcs = [
        fem.dirichletbc(ux0,     fem.locate_dofs_geometrical((W.sub(0).sub(0), V_ux), _face(0, 0.)), W.sub(0).sub(0)),
        fem.dirichletbc(uy0,     fem.locate_dofs_geometrical((W.sub(0).sub(1), V_uy), _face(1, 0.)), W.sub(0).sub(1)),
        fem.dirichletbc(uz0,     fem.locate_dofs_geometrical((W.sub(0).sub(2), V_uz), _face(2, 0.)), W.sub(0).sub(2)),
        fem.dirichletbc(ux_load, fem.locate_dofs_geometrical((W.sub(0).sub(0), V_ux), _face(0, 1.)), W.sub(0).sub(0)),
        fem.dirichletbc(p0,      fem.locate_dofs_geometrical((W.sub(1), Q_p),          _origin),     W.sub(1)),
    ]
    up      = fem.Function(W)
    macro_H = fem.Constant(msh, default_scalar_type(np.zeros((3, 3))))
    R, K    = _build_forms(msh, W, up, macro_H)
    problem = fem_petsc.NonlinearProblem(R, up, bcs=bcs, J=K,
                                         petsc_options_prefix="nonperiodic_",
                                         petsc_options=_SNES_OPTS)
    ok = True
    for step, alpha in enumerate(np.linspace(0., 1., NSTEPS + 1)):
        ux_load.x.array[:] = alpha * EXX; ux_load.x.scatter_forward()
        if step == 0:
            continue
        try:
            problem.solve()
            reason = problem.solver.getConvergedReason()
            iters  = problem.solver.getIterationNumber()
            print(f"[nonperiodic] step={step} alpha={alpha:.4f} reason={reason} iters={iters}")
        except PETSc.Error as exc:
            print(f"[nonperiodic] step={step} alpha={alpha:.4f} FAILED: {exc}")
            ok = False; break
    return ok

# ── Case 2: dolfinx_mpc periodic BCs (diverges — this is the bug) ────────────
def run_periodic():
    comm = MPI.COMM_WORLD
    msh  = mesh.create_unit_cube(comm, NCELLS, NCELLS, NCELLS)
    u_el = basix.ufl.element(basix.ElementFamily.P, msh.basix_cell(), 2,
                              shape=(3,), dtype=default_scalar_type)
    p_el = basix.ufl.element(basix.ElementFamily.P, msh.basix_cell(), 1,
                              dtype=default_scalar_type)
    W = fem.functionspace(msh, basix.ufl.mixed_element([u_el, p_el]))
    mpc  = dolfinx_mpc.MultiPointConstraint(W)
    bcs  = []
    for comp in range(3):
        V_ui, ui0 = _zero_fn(W.sub(0).sub(comp))
        dofs = fem.locate_dofs_geometrical((W.sub(0).sub(comp), V_ui), _origin)
        bc   = fem.dirichletbc(ui0, dofs, W.sub(0).sub(comp))
        bcs.append(bc)
        for axis in range(3):
            mpc.create_periodic_constraint_geometrical(
                W.sub(0).sub(comp),
                _periodic_indicator(axis),
                _periodic_relation(axis),
                [bc], default_scalar_type(1.), default_scalar_type(TOL),
            )
    Q_p, p0 = _zero_fn(W.sub(1))
    bcs.append(fem.dirichletbc(p0,
                               fem.locate_dofs_geometrical((W.sub(1), Q_p), _origin),
                               W.sub(1)))
    mpc.finalize()
    up      = fem.Function(mpc.function_space)
    macro_H = fem.Constant(msh, default_scalar_type(np.zeros((3, 3))))
    R, K    = _build_forms(msh, W, up, macro_H)
    problem = dolfinx_mpc.NonlinearProblem(
        R, up, mpc=mpc, bcs=bcs, J=K,
        petsc_options={**_SNES_OPTS, "snes_error_if_not_converged": False},
    )
    ok = True
    for step, alpha in enumerate(np.linspace(0., 1., NSTEPS + 1)):
        macro_H.value[...] = default_scalar_type(np.diag([alpha * EXX, 0., 0.]))
        if step == 0:
            continue
        _, reason, iters = problem.solve()
        converged = reason > 0
        print(f"[periodic]    step={step} alpha={alpha:.4f} reason={reason} iters={iters} converged={converged}")
        if not converged:
            ok = False
    return ok

def _print_versions():
    import dolfinx, basix, ufl, petsc4py, mpi4py, numpy as np
    libs = [
        ("dolfinx",  getattr(dolfinx,  "__version__", "n/a")),
        ("basix",    getattr(basix,    "__version__", "n/a")),
        ("ufl",      getattr(ufl,      "__version__", "n/a")),
        ("petsc4py", getattr(petsc4py, "__version__", "n/a")),
        ("mpi4py",   getattr(mpi4py,   "__version__", "n/a")),
        ("numpy",    np.__version__),
    ]
    for name in ("dolfinx_mpc",):
        try:
            libs.append((name, importlib.metadata.version(name)))
        except Exception:
            import dolfinx_mpc as _m
            libs.append((name, getattr(_m, "__version__", "n/a")))
    print("Library versions:")
    for name, ver in libs:
        print(f"  {name:<14} {ver}")
    print()

if __name__ == "__main__":
    _print_versions()
    print("=" * 60)
    print("Case 1: non-periodic (Dirichlet BCs)")
    print("=" * 60)
    ok_np = run_nonperiodic()
    print("=" * 60)
    print("Case 2: periodic (dolfinx_mpc geometrical constraints)")
    print("=" * 60)
    ok_p = run_periodic()
    print("=" * 60)
    print(f"nonperiodic converged : {ok_np}")
    print(f"periodic    converged : {ok_p}")

---

Output

Library versions:
  dolfinx        0.11.0.dev0
  basix          0.11.0.dev0
  ufl            2025.3.0.dev0
  petsc4py       3.24.5
  mpi4py         4.1.1
  numpy          2.4.3
  dolfinx_mpc    0.10.0.dev0

============================================================
Case 1: non-periodic (Dirichlet BCs)
============================================================
[nonperiodic] step=1 alpha=0.1667 reason=2 iters=3
[nonperiodic] step=2 alpha=0.3333 reason=2 iters=2
[nonperiodic] step=3 alpha=0.5000 reason=2 iters=2
[nonperiodic] step=4 alpha=0.6667 reason=2 iters=2
[nonperiodic] step=5 alpha=0.8333 reason=2 iters=2
[nonperiodic] step=6 alpha=1.0000 reason=2 iters=2
============================================================
Case 2: periodic (dolfinx_mpc geometrical constraints)
============================================================
[periodic]    step=1 alpha=0.1667 reason=-6 iters=1 converged=False
[periodic]    step=2 alpha=0.3333 reason=-6 iters=0 converged=False
[periodic]    step=3 alpha=0.5000 reason=-6 iters=0 converged=False
[periodic]    step=4 alpha=0.6667 reason=-6 iters=0 converged=False
[periodic]    step=5 alpha=0.8333 reason=-6 iters=0 converged=False
[periodic]    step=6 alpha=1.0000 reason=-6 iters=0 converged=False
============================================================
nonperiodic converged : True
periodic    converged : False

reason = -6 → SNES_DIVERGED_LINE_SEARCH at every load step.

---

Expected behaviour

The periodic and non-periodic formulations are mathematically equivalent for a
homogeneous single-material cube (no inclusions). Both should converge in a small
number of Newton iterations.

Actual behaviour

dolfinx_mpc.NonlinearProblem diverges immediately with
SNES_DIVERGED_LINE_SEARCH (reason = -6) at every load step,
even at very small applied strains ($\overline{H}_{11} \approx 0.005$).

Any clues what is wrong with the setup ? Maybe the way the mpc is setup?