fix(constrained): guard empty stratum IS in interior-multiplier reduction (#291)

src/underworld3/cython/petsc_generic_snes_solvers.pyx

@@ -6997,8 +6997,13 @@ class SNES_Stokes_SaddlePt(SolverBaseClass):
             fid_h = cbc.lam._solver_field_id
             bvalue = self.mesh.boundaries[cbc.boundary].value
             bd_is = dm.getLabel("UW_Boundaries").getStratumIS(bvalue)
+            # Guard for parallel: a rank owning no points with this label value
+            # gets back a valid non-None PETSc.IS with size 0 (bool(bd_is) = True,
+            # bd_is.handle != 0). Calling .getIndices() on that IS segfaults.
+            # `if bd_is` catches the null-handle case; `.getSize() > 0` catches
+            # the valid-but-empty case that None / handle checks miss (issue #291).
             keep = set()
-            if bd_is is not None:
+            if bd_is and bd_is.getSize() > 0:
                 for bp in bd_is.getIndices().tolist():
                     keep.update(dm.getTransitiveClosure(bp)[0].tolist())
             iset = interior_pts.setdefault(fid_h, set())

src/underworld3/utilities/rotated_bc.py

@@ -31,6 +31,30 @@
 _ROT_SOLVE_COUNT = 0
 
 
+def _warn_if_ksp_diverged(ksp, kind):
+    """Emit a rank-0 warning if a KSP finished on a KSP_DIVERGED_* reason
+    (negative converged-reason). ``ksp.solve`` never raises on divergence, so
+    without this the linear rotated-freeslip solvers here would return a
+    partial answer to the caller silently — which is exactly what let the 3D
+    rotation-nullspace bug (#306) go unnoticed until it produced clearly-wrong
+    physics in a downstream test."""
+    reason = int(ksp.getConvergedReason())
+    # Convention: > 0 == converged, < 0 == diverged, 0 == KSP_CONVERGED_ITERATING
+    # (should not happen after ksp.solve() returns; warn if it does). Matches
+    # petsc_generic_snes_solvers.pyx:2979.
+    if reason > 0:
+        return
+    its = int(ksp.getIterationNumber())
+    try:
+        rnorm = float(ksp.getResidualNorm())
+    except Exception:
+        rnorm = float("nan")
+    from underworld3 import mpi
+    mpi.pprint(f"[rotated_bc] WARNING: {kind} KSP did NOT converge "
+               f"(reason={reason}, iterations={its}, |r|={rnorm:.3e}); "
+               f"proceeding with the last iterate.")
+
+
 # --------------------------------------------------------------------------- #
 #  Rotation construction
 # --------------------------------------------------------------------------- #
@@ -291,6 +315,7 @@ def solve_rotated_freeslip(solver, boundaries, remove_rotation_gauge=True, verbo
         pc = ksp.getPC(); pc.setType("lu"); pc.setFactorSolverType("mumps")
         Uhat = dm.createGlobalVec(); ksp.solve(bhat, Uhat)   # returned in info → own it
         ksp_reason = ksp.getConvergedReason(); ksp_its = ksp.getIterationNumber()
+        _warn_if_ksp_diverged(ksp, kind="rotated direct-LU")
     else:
         Mp = _pressure_mass_schur_pmat(solver)
         Uhat, ksp_reason, ctx = _solve_rotated_iterative(
@@ -797,6 +822,7 @@ def _solve_rotated_iterative(solver, Ahat, bhat, Q, Qt, normal_rows, verbose=Fal
 
     Uhat = Ahat.createVecRight(); Uhat.set(0.0)
     ksp.solve(bhat, Uhat)
+    _warn_if_ksp_diverged(ksp, kind="rotated fieldsplit-Schur")
     # An identity constraint row in an ITERATIVE solve only drives its residual
     # (= û_i) below tolerance, so û_i ~ tol, not exactly 0. Because zeroRowsColumns
     # made these DOFs fully decoupled (row AND column zeroed), û_i affects no other

tests/parallel/test_1017_custom_mg_parallel_mpi.py

@@ -126,14 +126,8 @@ def test_parallel_custom_fmg_vector():
     assert proj.snes.getConvergedReason() > 0
 
 
-@pytest.mark.skip(reason="Stokes_Constrained is not parallel-safe yet — it "
-                  "segfaults at np>1 independently of custom-P, in the "
-                  "interior-multiplier section reduction (issue #291; canonical "
-                  "test_1062_constrained_solcx also segfaults at np=2, plain GAMG). "
-                  "custom-P on the constrained velocity block works in SERIAL "
-                  "(see test_1017_custom_mg_stokes.test_custom_fmg_stokes_constrained); "
-                  "this auto-enables once #291 is fixed.")
 @pytest.mark.mpi(min_size=2)
+@pytest.mark.slow
 def test_parallel_custom_fmg_stokes_constrained():
     """Stokes_Constrained (free-slip multipliers, grouped [p,h] split) custom-P on
     the velocity block in parallel — velocity must match the analytic SolCx."""

tests/parallel/test_1062_constrained_stratum_guard_parallel.py

@@ -0,0 +1,75 @@
+"""Parallel regression for the empty-stratum guard in
+``SNES_Stokes_Constrained._constrain_interior_multipliers_in_section``
+(issue #291, PR #318).
+
+Before #318, ``Stokes_Constrained`` segfaulted at np>1 in the interior-
+multiplier section reduction: on a rank owning zero points with a given
+boundary label value, ``dm.getLabel("UW_Boundaries").getStratumIS(bvalue)``
+returned a valid non-None PETSc IS with ``getSize() == 0``, the existing
+``if bd_is is not None`` guard passed through, and the following
+``bd_is.getIndices()`` crashed with SIGSEGV.
+
+The canonical ``tests/test_1062_constrained_solcx.py`` uses the SolCx
+eta_B=1e6 benchmark which is far too slow to serve as a fast np≥2
+regression check (Constrained is ~34× slower than Nitsche per #244).
+This test uses the same solver family but at **unit viscosity** and a
+16×16 grid, so both the guard and a bit-identical parallel-vs-serial
+diagnostic land in ~10 seconds at np=2.
+
+Run:
+    mpirun -n 2 python -m pytest --with-mpi \\
+        tests/parallel/test_1062_constrained_stratum_guard_parallel.py
+    mpirun -n 4 python -m pytest --with-mpi \\
+        tests/parallel/test_1062_constrained_stratum_guard_parallel.py
+"""
+import numpy as np
+import sympy
+import pytest
+
+import underworld3 as uw
+
+pytestmark = [pytest.mark.mpi(min_size=2), pytest.mark.timeout(120)]
+
+# Serial reference (unit viscosity, sinusoidal body force, 16×16 quad box, free-
+# slip on all four walls, tol 1e-8). Re-derive with `python <thisfile>`.
+GOLDEN_U_INF = 2.533050e-2
+
+
+def _solve():
+    """Build + solve the constrained SolCx-lite problem; return |u|_inf."""
+    mesh = uw.meshing.StructuredQuadBox(
+        elementRes=(16, 16), minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), qdegree=3)
+    s = uw.systems.Stokes_Constrained(mesh)
+    s.constitutive_model = uw.constitutive_models.ViscousFlowModel
+    s.constitutive_model.Parameters.shear_viscosity_0 = 1.0
+    x, y = mesh.X
+    s.bodyforce = sympy.Matrix(
+        [[0.0, sympy.cos(sympy.pi * x) * sympy.sin(sympy.pi * y)]])
+    # Free-slip on all four walls via in-saddle multipliers. On a 2-way axis-
+    # aligned split at np=2 this puts one axis-aligned pair (Left/Right or
+    # Bottom/Top) with all points on one rank and zero on the other — the
+    # exact configuration that triggered #291's SIGSEGV.
+    s.add_constraint_bc("Left",   g=0.0, normal=sympy.Matrix([[-1.0, 0.0]]))
+    s.add_constraint_bc("Right",  g=0.0, normal=sympy.Matrix([[ 1.0, 0.0]]))
+    s.add_constraint_bc("Bottom", g=0.0, normal=sympy.Matrix([[ 0.0,-1.0]]))
+    s.add_constraint_bc("Top",    g=0.0, normal=sympy.Matrix([[ 0.0, 1.0]]))
+    s.petsc_use_pressure_nullspace = True
+    s.tolerance = 1.0e-8
+    s.solve()
+    return float(np.max(np.abs(s.Unknowns.u.data)))
+
+
+def test_constrained_parallel_no_segfault():
+    """Pre-#318 this test would SIGSEGV during solve() at np>1. Post-#318 it
+    solves and reproduces the serial |u|_inf to a tight tolerance."""
+    u_inf = _solve()
+    assert np.isclose(u_inf, GOLDEN_U_INF, rtol=1e-3, atol=0), (
+        f"|u|_inf differs serial vs np={uw.mpi.size}: "
+        f"golden={GOLDEN_U_INF:.6e} got={u_inf:.6e}")
+
+
+if __name__ == "__main__":
+    # Recompute the serial GOLDEN: `python <thisfile>`.
+    u_inf = _solve()
+    if uw.mpi.rank == 0:
+        print(f"DIAG |u|_inf = {u_inf:.6e}")