diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index f0b4ef19d..048b52a65 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()) diff --git a/src/underworld3/utilities/rotated_bc.py b/src/underworld3/utilities/rotated_bc.py index e6afd60f1..8ccc1699a 100644 --- a/src/underworld3/utilities/rotated_bc.py +++ b/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 diff --git a/tests/parallel/test_1017_custom_mg_parallel_mpi.py b/tests/parallel/test_1017_custom_mg_parallel_mpi.py index 3a3e29a56..d3f925409 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 000000000..16f9ce7ee --- /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}")