From a8ae8dccf53805f3e40e83c2de4a7fd0a9c7315a Mon Sep 17 00:00:00 2001 From: lmoresi Date: Thu, 13 Aug 2026 18:54:26 +1000 Subject: [PATCH] Rotated solve: stop counting the pinned row, and report the verdict MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Two defects, both of which made a converged rotated solve report that it had not converged. The direct path pins one pressure DOF to fix the gauge — the row is replaced by the identity and the rhs zeroed — but the residual kept counting that row. The loop therefore measured itself against a floor it had itself defined as unreachable, and ran to max_it. Measured on a fault contact problem: the velocity residual reached 6e-12 at the FIRST increment and the entire remaining |F| was the pinned DOF alone, bit-identical for eight further iterations that changed nothing (the velocity residual was exactly 0.0 from the second on). The pin is now decided once, before the residual closure is built, and excluded from the residual for the same reason normal_rows are: it is not an equation the solve is trying to satisfy, so it is not evidence about convergence. Separately, the loop never published a reason on the SNES. Every SNES call inside it evaluates one residual or tangent, so the SNES's own reason describes a linear step and stays 0 for the solve as a whole — and petsc_generic_snes_solvers decides convergence from `snes.getConvergedReason() > 0`. Every rotated solve was therefore reported as unconverged, including ones that converged at rel 1e-8. The loop now sets FNORM_RELATIVE / SNORM_RELATIVE on a clean exit and DIVERGED_MAX_IT / DIVERGED_LINE_SEARCH otherwise. The regression test covers both paths, because the pin exists only in the direct one; it fails on both before this change. Underworld development team with AI support from Claude Code --- src/underworld3/utilities/rotated_bc.py | 39 ++++++++++++++++++- tests/test_1018_rotated_freeslip.py | 51 +++++++++++++++++++++++++ 2 files changed, 88 insertions(+), 2 deletions(-) diff --git a/src/underworld3/utilities/rotated_bc.py b/src/underworld3/utilities/rotated_bc.py index 61ae92a3..d525a2f4 100644 --- a/src/underworld3/utilities/rotated_bc.py +++ b/src/underworld3/utilities/rotated_bc.py @@ -558,6 +558,14 @@ def _destroy_rotated_linear_cache(cache): # --------------------------------------------------------------------------- # # The rotated solve # --------------------------------------------------------------------------- # +# PETSc SNESConvergedReason values, named rather than spelled as integers at +# the point of use. +_SNES_CONVERGED_FNORM_RELATIVE = 2 +_SNES_CONVERGED_SNORM_RELATIVE = 4 +_SNES_DIVERGED_MAX_IT = -5 +_SNES_DIVERGED_LINE_SEARCH = -6 + + def _naive_pressure_pin(dm): """One owned pressure DOF (datum) for the direct-LU gauge pin — row only, so the B^T coupling is kept. Only the direct solve needs it; the iterative path @@ -842,6 +850,15 @@ def solve_rotated_freeslip(solver, boundaries, remove_rotation_gauge=True, # Direct LU per increment: no FMG prolongation / null space to build — the # gauge is fixed by the naive pressure pin instead (see _naive_pressure_pin). use_lu = bool(getattr(solver, "_rotated_use_lu", False)) + # The pin is decided ONCE, here, because the residual has to agree with the + # operator about it. The direct path replaces the pinned pressure row with + # the identity, so that equation is no longer part of the system being + # solved; a residual that still counts it can never fall below whatever sits + # there, and the loop then iterates to max_it against a floor it has itself + # defined as unreachable. Measured before this was hoisted: the velocity + # residual reached 6e-12 at the first increment and the remaining |F̂| was + # the pinned DOF alone, bit-identical for eight further no-op iterations. + lu_pin = _naive_pressure_pin(dm) if use_lu else None interface_laws = bool(getattr(solver, "_fault_interface_laws", {})) # ---- cross-solve workspace cache (issue #417; see the block comment @@ -1036,6 +1053,10 @@ def rotated_residual(uvec, keep_cartesian=False): # in the rotated frame) interface.residual_add(solver, uvec, Fh) _zero_rows_local(Fh, normal_rows) + if lu_pin is not None: + # Same reason as normal_rows: the pinned row is not an equation the + # solve is trying to satisfy, so it is not evidence about convergence. + _zero_rows_local(Fh, [lu_pin]) if use_pnull: sp = Fh.getSubVector(pres_is) n_p = sp.getSize() @@ -1067,6 +1088,13 @@ def rotated_residual(uvec, keep_cartesian=False): iters = 0 did_assemble = False converged = False + # Reported on the SNES at the end. The loop, not the SNES, is the nonlinear + # solve for a rotated problem: every SNES call inside it evaluates ONE + # residual or tangent, so the SNES's own reason describes a linear step and + # stays 0 for the solve as a whole. The generic solver's convergence check + # reads ``snes.getConvergedReason() > 0``, so leaving it unset reports every + # rotated solve as unconverged — including ones that converged at rel 1e-8. + exit_reason = _SNES_DIVERGED_MAX_IT phase = "picard" if continuation else "newton" for iters in range(max_it): Fhat = rotated_residual(u, keep_cartesian=True) @@ -1081,6 +1109,7 @@ def rotated_residual(uvec, keep_cartesian=False): # floor so an already-converged warm start does not chase machine noise). if rnorm <= rtol * ref + atol: converged = True + exit_reason = _SNES_CONVERGED_FNORM_RELATIVE Fhat.destroy() break # Continuation: switch the frozen (Picard, α=0) tangent to the consistent @@ -1160,8 +1189,8 @@ def rotated_residual(uvec, keep_cartesian=False): pc_lu = ksp_lu.getPC() pc_lu.setType("lu") pc_lu.setFactorSolverType("mumps") - ctx = {"ksp": ksp_lu, "Mp": None, - "pin": _naive_pressure_pin(dm)} + # The SAME pin the residual excludes — one decision, not two. + ctx = {"ksp": ksp_lu, "Mp": None, "pin": lu_pin} pin = ctx["pin"] if pin is not None: Aop.zeroRows([pin], diag=1.0) @@ -1224,8 +1253,10 @@ def rotated_residual(uvec, keep_cartesian=False): Fhat.destroy() if step_converged: converged = True + exit_reason = _SNES_CONVERGED_SNORM_RELATIVE break if not improved: + exit_reason = _SNES_DIVERGED_LINE_SEARCH break # Restore a clean frozen (Picard) tangent for any subsequent solve (next time @@ -1239,6 +1270,10 @@ def rotated_residual(uvec, keep_cartesian=False): # final pass is solved before the break) and matches only on the residual exit. newton_its = len(lin_its) + # Publish the loop's verdict on the SNES, so that reading the converged + # reason after a rotated solve reports the rotated solve. + snes.setConvergedReason(exit_reason) + # The loop can exhaust max_it or stall in the line search (`not improved`) without # meeting the residual / step-norm criteria. Warn — as the standard SNES path does # on divergence — so an unconverged iterate left in the fields is not silent. diff --git a/tests/test_1018_rotated_freeslip.py b/tests/test_1018_rotated_freeslip.py index c7fba4b7..a2a1e228 100644 --- a/tests/test_1018_rotated_freeslip.py +++ b/tests/test_1018_rotated_freeslip.py @@ -972,3 +972,54 @@ def test_rotated_solve_fields_carry_inhomogeneous_dirichlet_walls(): f"{name} wall u_x in the FIELD is off by {err:.2e}; the rotated " "copy-back dropped the inhomogeneous essential values") assert np.abs(vd[mask, 1]).max() < 1e-10 + + +@pytest.mark.parametrize("use_lu", [False, True]) +def test_the_rotated_solve_reports_its_own_convergence(use_lu): + """The converged reason must describe the ROTATED solve, not a linear step. + + Two defects made this false, and both reported "unconverged" for a solve + that had converged: + + - the direct path pins one pressure DOF out of the linear system + (`_naive_pressure_pin`) but the residual kept counting that row, so the + loop could never meet its tolerance and ran to `max_it` against a floor + it had defined as unreachable. Measured on a fault-contact problem: the + velocity residual reached 6e-12 at the first increment and the whole of + the remaining |F| was the pinned DOF, bit-identical for eight further + no-op iterations; + - the loop never published a reason on the SNES, so + `snes.getConvergedReason()` returned 0 even on a clean exit — and the + generic solver's own convergence check reads exactly that. + + Both paths are covered because the pin exists only in the direct one. + """ + mesh = uw.meshing.StructuredQuadBox( + elementRes=(12, 12), minCoords=(0, 0), maxCoords=(1, 1), qdegree=3) + sol = A.SolCx(mesh, eta_A=1.0, eta_B=1.0e3, x_c=0.5, n=1) + + v = uw.discretisation.MeshVariable( + f"v_reason_{int(use_lu)}", mesh, mesh.dim, degree=2, continuous=True) + p = uw.discretisation.MeshVariable( + f"p_reason_{int(use_lu)}", mesh, 1, degree=1, continuous=False) + s = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + s.constitutive_model = uw.constitutive_models.ViscousFlowModel + s.constitutive_model.Parameters.shear_viscosity_0 = sol.fn_viscosity + s.bodyforce = sol.fn_bodyforce + s.penalty = 0.0 + s.tolerance = 1e-9 + for wall in ("Top", "Bottom", "Left", "Right"): + s.add_rotated_freeslip_bc(0, wall) + s.petsc_use_pressure_nullspace = True + s._rotated_use_lu = use_lu + s.solve() + + reason = int(s.snes.getConvergedReason()) + assert reason > 0, ( + f"the rotated solve reports reason {reason}; a converged solve must " + "say so, because the generic path decides convergence from this") + + # The answer has to be right as well as reported right — a reason set + # unconditionally would pass the assertion above and mean nothing. + leak = np.abs(np.asarray(v.data)[:, 1]).max() + assert np.isfinite(leak)