Rotated solve: stop counting the pinned row, and report the verdict

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.

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)