From 56d3d37fa9a7794c9a79d51eb3c309773e068e81 Mon Sep 17 00:00:00 2001 From: Louis Moresi Date: Sat, 4 Jul 2026 13:06:18 +1000 Subject: [PATCH 1/4] fix(rotated-bc): warn on non-converged linear KSP `ksp.solve` never raises on `KSP_DIVERGED_*`, so the two linear paths in `rotated_bc` (the default fieldsplit-Schur iterative solve at `_solve_rotated_iterative` and the opt-in direct MUMPS LU at `solve_rotated_freeslip`) both returned partial answers silently. That silence is what let the 3D rotation-nullspace bug fixed by #306 look like a working solve for as long as it did: the outer Krylov hit its `ksp_max_it` ceiling (300 in the pre-#306 code) with a residual well above tolerance, and neither the KSP nor the caller flagged it. The solution then propagated into downstream diagnostics that eventually exposed it as wrong physics rather than as a solver failure. Add a small `_warn_if_ksp_diverged` helper and call it immediately after each of the two `ksp.solve` sites. Rank-0 warning only (via `uw.mpi.pprint`), same style as the existing nonlinear-driver warning at the end of `solve_rotated_freeslip_nonlinear`. No behaviour change on the happy path; non-convergence is now loud instead of silent. test_1018_rotated_freeslip.py: 14/14 pass (2D behaviour unchanged). Underworld development team with AI support from Claude Code --- src/underworld3/utilities/rotated_bc.py | 23 +++++++++++++++++++++++ 1 file changed, 23 insertions(+) diff --git a/src/underworld3/utilities/rotated_bc.py b/src/underworld3/utilities/rotated_bc.py index e6afd60f..47278a9c 100644 --- a/src/underworld3/utilities/rotated_bc.py +++ b/src/underworld3/utilities/rotated_bc.py @@ -31,6 +31,27 @@ _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()) + 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 +312,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 +819,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 From 89cf405c4896ba915a8ab4e5869dc90ef64605c9 Mon Sep 17 00:00:00 2001 From: Louis Moresi Date: Sun, 5 Jul 2026 19:44:45 +1000 Subject: [PATCH 2/4] review(Copilot): treat KSP_CONVERGED_ITERATING (reason=0) as non-converged MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Copilot review on #315 noted the convention split: * KSP.getConvergedReason() > 0 → converged * KSP.getConvergedReason() == 0 → KSP_CONVERGED_ITERATING (initial state; unexpected after ksp.solve() returns) * KSP.getConvergedReason() < 0 → diverged The previous `reason >= 0` early-return silently accepted a hypothetical `reason == 0` outcome, undermining the guard's purpose. Change to `reason > 0` so a stray `KSP_CONVERGED_ITERATING` also fires the warning, matching the project convention documented in petsc_generic_snes_solvers.pyx:2979. test_1018_rotated_freeslip.py: 14/14 pass (unchanged). Underworld development team with AI support from Claude Code --- src/underworld3/utilities/rotated_bc.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/underworld3/utilities/rotated_bc.py b/src/underworld3/utilities/rotated_bc.py index 47278a9c..8ccc1699 100644 --- a/src/underworld3/utilities/rotated_bc.py +++ b/src/underworld3/utilities/rotated_bc.py @@ -39,7 +39,10 @@ def _warn_if_ksp_diverged(ksp, kind): rotation-nullspace bug (#306) go unnoticed until it produced clearly-wrong physics in a downstream test.""" reason = int(ksp.getConvergedReason()) - if reason >= 0: + # 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: From d50ab532503522960ab20049bc20ccf26643d2fb Mon Sep 17 00:00:00 2001 From: Louis Moresi Date: Sun, 5 Jul 2026 22:04:40 +1000 Subject: [PATCH 3/4] fix(constrained): guard empty stratum IS in interior-multiplier reduction (#291) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `Stokes_Constrained` segfaulted at np>1 in `_constrain_interior_multipliers_in_section`: on a rank owning zero points with a given boundary label value, `dm.getLabel("UW_Boundaries").getStratumIS(bvalue)` returns a VALID PETSc IS with `getSize() == 0` — not None, non-zero handle, `bool(IS) == True`. The existing `if bd_is is not None:` guard passed through, and the subsequent `bd_is.getIndices()` crashed with SIGSEGV in the PETSc IS binding. The failure was reliably reproduced by #291's repro (StructuredQuadBox + free-slip on all four walls, np=2) and by a lighter unit-viscosity variant which now solves bit-identically at np=1/2/4 with the fix. Fix: replace `if bd_is is not None:` with `if bd_is and bd_is.getSize() > 0:`. `bool(bd_is)` short-circuits the null-handle case (project convention, see `utilities/boundary_flux.py:70`); `.getSize() > 0` catches the valid-but-empty case that neither the None check nor the truthiness check protects against. ### Verification Lightweight unit-viscosity repro (16x16 quad, sinusoidal body force, free-slip on all four walls; `_reduce_interior_multiplier = True` default): np=1: |u|_inf = 2.533e-02 (6.6 s) np=2: |u|_inf = 2.533e-02 (8.5 s) [previously: SIGSEGV] np=4: |u|_inf = 2.533e-02 (5.0 s) [previously: SIGSEGV] Bit-identical parallel-vs-serial in this regime. `tests/test_1062_constrained_solcx.py` (canonical SolCx eta_B=1e6 test): serial passes (198 s); np=2 with-mpi rerun exceeds the default 300 s timeout under this session's compute — a longer-timeout parallel check is warranted as follow-up (probably belongs in `tests/parallel/`). Underworld development team with AI support from Claude Code --- src/underworld3/cython/petsc_generic_snes_solvers.pyx | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index f0b4ef19..048b52a6 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/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()) From 958f6f41efc14efc132b4a6ab3359dc6898d444e Mon Sep 17 00:00:00 2001 From: Louis Moresi Date: Sun, 5 Jul 2026 22:15:42 +1000 Subject: [PATCH 4/4] =?UTF-8?q?test(#291):=20add=20lightweight=20np?= =?UTF-8?q?=E2=89=A52=20regression=20+=20un-skip=20custom-FMG=20parallel?= =?UTF-8?q?=20test?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Copilot review on #318 rightly pointed out that a parallel-only SIGSEGV regression should not merge without a re-enabled parallel test. Two changes: ## 1. New lightweight parallel regression (tests/parallel/test_1062_constrained_stratum_guard_parallel.py) Unit-viscosity 16×16 quad box, sinusoidal body force, free-slip on all four walls at np≥2. On a 2-way axis-aligned partition (the default) this puts one full boundary on each rank and zero points on the other — the exact empty-stratum configuration that #291's guard now handles. Runs in ~10 s at np=2 (vs the canonical test_1062 SolCx eta_B=1e6 which timed out at 300 s in this session's parallel runs, per #244). Asserts `|u|_inf` matches the serial GOLDEN of 2.53305e-02 to rtol 1e-3 (actual serial-vs-parallel spread on this workload is < 1e-9). ## 2. Un-skip test_parallel_custom_fmg_stokes_constrained (tests/parallel/test_1017_custom_mg_parallel_mpi.py) Its skip reason explicitly cited "#291; auto-enables once #291 is fixed." The guard fix is now in place, so the test can run. Marked `slow` (SolCx eta_B=1e6 is not a fast benchmark) so CI selects it via `pytest -m slow` or `-m "not slow"` per environment. Underworld development team with AI support from Claude Code --- .../test_1017_custom_mg_parallel_mpi.py | 8 +- ...1062_constrained_stratum_guard_parallel.py | 75 +++++++++++++++++++ 2 files changed, 76 insertions(+), 7 deletions(-) create mode 100644 tests/parallel/test_1062_constrained_stratum_guard_parallel.py diff --git a/tests/parallel/test_1017_custom_mg_parallel_mpi.py b/tests/parallel/test_1017_custom_mg_parallel_mpi.py index 3a3e29a5..d3f92540 100644 --- a/tests/parallel/test_1017_custom_mg_parallel_mpi.py +++ b/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.""" diff --git a/tests/parallel/test_1062_constrained_stratum_guard_parallel.py b/tests/parallel/test_1062_constrained_stratum_guard_parallel.py new file mode 100644 index 00000000..16f9ce7e --- /dev/null +++ b/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 `. +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 `. + u_inf = _solve() + if uw.mpi.rank == 0: + print(f"DIAG |u|_inf = {u_inf:.6e}")