fix(solve-report): honest rotated-path reports, KSP reason namespace, resume-anchor lifecycle

src/underworld3/cython/petsc_generic_snes_solvers.pyx

@@ -919,6 +919,9 @@ class SolverBaseClass(uw_object):
 
         self.natural_bcs = []
         self.essential_bcs = []
+        # A teardown means the next solve is a different discrete problem: a resume
+        # anchor from the old problem's ||F0|| must not survive into it.
+        self._resume_abs_target = None
 
         if self.snes is not None:
             self.snes.destroy()
@@ -985,6 +988,10 @@ class SolverBaseClass(uw_object):
         reason = int(self.snes.getConvergedReason())
         nl_its = int(self.snes.getIterationNumber())
         ksp_its = int(self.snes.getLinearSolveIterations())
+        # Sanctioned swallows (each optional-field read below): PETSc raises from these
+        # getters on SNES types/states that never computed the quantity (e.g. ksponly
+        # before a function evaluation, builds without the counter). The report field
+        # degrades to its honest "unavailable" value; nothing else is masked.
         try:
             fnorm = float(self.snes.getFunctionNorm())
         except Exception:
@@ -1022,27 +1029,47 @@ class SolverBaseClass(uw_object):
     def _capture_rotated_report(self, info):
         """Record a SolveReport for a rotated-free-slip solve, which runs its own
         ``ksp.solve()`` on the rotated operator (``utilities/rotated_bc.py``) and never
-        touches ``self.snes`` — so the generic reader would see stale SNES state. Build the
-        report from the rotated result dict (``ksp_reason`` / ``ksp_its``) so a rotated solve
-        also leaves a (read-only) report rather than a stale one. Best-effort, never raises."""
-        from underworld3.systems.solve_report import SolveReport, reason_string
-        try:
-            info = info or {}
-            reason = int(info.get("ksp_reason", 0))
-            ksp_its = int(info.get("ksp_its", 0))
-            nl_its = int(info.get("newton_its", 1))     # linear path = one outer solve
-            fnorm = float(info.get("rnorm", float("nan")))
-            report = SolveReport(
-                reason=reason, reason_str=reason_string(reason), converged=reason > 0,
-                nl_its=nl_its, ksp_its=ksp_its, fnorm=fnorm,
-                fnorm0=None, reduction=None, rho=None, fev=None, history=(), bounded=False,
-            )
-            self._solve_report = report
-            self._solve_history.append(report)
-            return report
-        except Exception:
-            self._solve_report = None
-            return None
+        touches ``self.snes`` — so the generic reader would see stale SNES state.
+
+        Handles BOTH rotated result dicts: the linear one-shot solve (``ksp_its`` an
+        int, one outer solve, no outer verdict) and the manual nonlinear Newton loop
+        (``nonlinear_iterations``, ``ksp_its`` a per-increment LIST, and an outer
+        ``converged`` flag from the residual/step-norm tests). ``ksp_reason`` is a KSP
+        code on both paths (the nonlinear dict carries the LAST increment's), so it is
+        rendered with the KSP reason table — the SNES table shares the integers but
+        names different outcomes. A malformed dict raises: this reader and the dicts in
+        ``rotated_bc.py`` are a contract, and a silent fallback here previously masked
+        a broken one."""
+        from underworld3.systems.solve_report import SolveReport, ksp_reason_string
+
+        info = info or {}
+        reason = int(info.get("ksp_reason", 0))
+        raw_its = info.get("ksp_its", 0)
+        if isinstance(raw_its, (list, tuple)):
+            ksp_its = int(sum(raw_its))                 # nonlinear: one count per Newton step
+        else:
+            ksp_its = int(raw_its)
+        nl_its = int(info.get("nonlinear_iterations", 1))   # linear path = one outer solve
+        if "converged" in info:
+            # The nonlinear loop's own verdict. The last KSP reason alone would mislabel
+            # a stalled Newton chain whose final linear solve happened to converge.
+            converged = bool(info["converged"])
+        else:
+            converged = reason > 0
+        rnorm = info.get("rnorm")
+        fnorm = float(rnorm) if rnorm is not None else float("nan")
+        rnorm0 = info.get("rnorm0")
+        fnorm0 = float(rnorm0) if rnorm0 else None
+        reduction = (fnorm / fnorm0) if (fnorm0 and fnorm == fnorm) else None
+        report = SolveReport(
+            reason=reason, reason_str=ksp_reason_string(reason), converged=converged,
+            nl_its=nl_its, ksp_its=ksp_its, fnorm=fnorm,
+            fnorm0=fnorm0, reduction=reduction, rho=None, fev=None, history=(),
+            bounded=False,
+        )
+        self._solve_report = report
+        self._solve_history.append(report)
+        return report
 
     def estimate_difficulty(self, max_nl_its, *, warm=True, **solve_kwargs):
         """Run a bounded, resumable solve to *estimate solver difficulty* and return its report.
@@ -1066,13 +1093,39 @@ class SolverBaseClass(uw_object):
             ``False`` starts a fresh chain from a zero initial guess.
         **solve_kwargs
             Forwarded to ``solve()`` (e.g. ``timestep``, ``picard`` for Stokes/VE).
+            ``zero_init_guess`` is owned by ``warm`` and ``divergence_retries`` by the
+            probe (always 0) — passing either raises.
 
         Returns
         -------
         SolveReport
             The (``bounded=True``) report for this capped solve; also available as
             ``self.solve_report``.
+
+        Notes
+        -----
+        Not available with rotated free-slip BCs: that path solves outside
+        ``self.snes``, so the cap and anchor cannot reach it. Under
+        ``consistent_jacobian="continuation"`` the cap applies to EACH stage (Picard
+        then Newton), so one probe may run up to twice ``max_nl_its``.
         """
+        if getattr(self, "_rotated_freeslip_bcs", None):
+            raise NotImplementedError(
+                "estimate_difficulty() is not available with rotated free-slip BCs: "
+                "that path solves outside self.snes (utilities/rotated_bc.py), so the "
+                "iteration cap and resume anchor cannot be applied to it."
+            )
+        if "zero_init_guess" in solve_kwargs:
+            raise TypeError(
+                "estimate_difficulty() sets the initial guess from warm= "
+                "(warm=False starts a fresh chain from zero); do not pass zero_init_guess."
+            )
+        if solve_kwargs.get("divergence_retries"):
+            raise TypeError(
+                "estimate_difficulty() does not accept divergence_retries: a capped probe "
+                "ends DIVERGED_MAX_IT by design, and retrying on it would run multiples of "
+                "the advertised cap."
+            )
         if not warm:
             self._resume_abs_target = None      # fresh chain: forget any prior anchor
         self._difficulty_probe = True
@@ -1083,17 +1136,25 @@ class SolverBaseClass(uw_object):
             self._difficulty_probe = False
             self._difficulty_max_it = None
 
-        # Anchor a NEW chain to the original ||F0||. Subsequent warm restarts via
-        # estimate_difficulty(warm=True) — including a large-cap call to run to
-        # completion — terminate at tolerance*||F0||, the same point an uninterrupted
-        # solve would, regardless of where the caps fell. The anchor is consulted ONLY
-        # under _difficulty_probe (see _snes_solve_with_retries), so plain solves are
-        # unaffected and cannot inherit a stale anchor.
-        if self._resume_abs_target is None and self._solve_report is not None:
-            f0 = self._solve_report.fnorm0
+        # Anchor lifecycle. A CONVERGED probe ends the chain: clear the anchor so it
+        # cannot leak into a later chain on different physics. This also covers the
+        # warm first call on an already-converged state — its ||F0|| is the *converged*
+        # residual, and anchoring to tolerance*||F_converged|| would set an unreachable
+        # target that silently degrades the chain to restart-relative semantics.
+        # An UNCONVERGED (capped) probe with no anchor yet establishes the chain,
+        # anchored to the original ||F0||: subsequent warm restarts — including a
+        # large-cap call to run to completion — terminate at tolerance*||F0||, the same
+        # point an uninterrupted solve would, regardless of where the caps fell. The
+        # anchor is consulted ONLY under _difficulty_probe (see
+        # _snes_solve_with_retries), so plain solves are unaffected.
+        report = self._solve_report
+        if report is not None and report.converged:
+            self._resume_abs_target = None
+        elif self._resume_abs_target is None and report is not None:
+            f0 = report.fnorm0
             if f0:
                 self._resume_abs_target = self.tolerance * float(f0)
-        return self._solve_report
+        return report
 
     def get_snes_diagnostics(self):
         """
@@ -1404,10 +1465,12 @@ class SolverBaseClass(uw_object):
             if _saved_tol is not None:
                 self.snes.setTolerances(rtol=_saved_tol[0], atol=_saved_tol[1],
                                         stol=_saved_tol[2], max_it=_saved_tol[3])
-
-        # Default-on difficulty report — single tail so every path is covered.
-        # bounded=True marks a report produced under an estimate_difficulty cap.
-        self._capture_solve_report(bounded=_probe)
+            # Default-on difficulty report — single tail so every path is covered,
+            # INCLUDING a solve that raises: the report then reflects the failed
+            # attempt's SNES state rather than leaving the previous solve's converged
+            # report behind for recovery logic to misread. bounded=True marks a report
+            # produced under an estimate_difficulty cap.
+            self._capture_solve_report(bounded=_probe)
 
     def _continuation_solve(self, gvec, verbose=False):
         """Picard -> Newton continuation via the constants[]-routed alpha.

src/underworld3/systems/solve_report.py

@@ -42,6 +42,40 @@ def reason_string(reason: int) -> str:
     return REASON_STRINGS.get(int(reason), f"UNKNOWN_{reason}")
 
 
+# PETSc KSPConvergedReason codes -> short names. This is a DIFFERENT namespace from the
+# SNES table above: the integers overlap but name different outcomes (e.g. -3 is
+# SNES DIVERGED_LINEAR_SOLVE but KSP DIVERGED_MAX_IT). Rotated free-slip solves report
+# KSP codes and must be rendered with this table. Names carry the KSP_ prefix so a
+# report string is unambiguous about which namespace it came from. Verified against
+# petsc4py's PETSc.KSP.ConvergedReason (see test_1055).
+KSP_REASON_STRINGS = {
+    1: "KSP_CONVERGED_RTOL_NORMAL_EQUATIONS",
+    2: "KSP_CONVERGED_RTOL",
+    3: "KSP_CONVERGED_ATOL",
+    4: "KSP_CONVERGED_ITS",
+    5: "KSP_CONVERGED_NEG_CURVE",
+    6: "KSP_CONVERGED_STEP_LENGTH",
+    7: "KSP_CONVERGED_HAPPY_BREAKDOWN",
+    9: "KSP_CONVERGED_ATOL_NORMAL_EQUATIONS",
+    0: "KSP_ITERATING",
+    -2: "KSP_DIVERGED_NULL",
+    -3: "KSP_DIVERGED_MAX_IT",
+    -4: "KSP_DIVERGED_DTOL",
+    -5: "KSP_DIVERGED_BREAKDOWN",
+    -6: "KSP_DIVERGED_BREAKDOWN_BICG",
+    -7: "KSP_DIVERGED_NONSYMMETRIC",
+    -8: "KSP_DIVERGED_INDEFINITE_PC",
+    -9: "KSP_DIVERGED_NANORINF",
+    -10: "KSP_DIVERGED_INDEFINITE_MAT",
+    -11: "KSP_DIVERGED_PCSETUP_FAILED",
+}
+
+
+def ksp_reason_string(reason: int) -> str:
+    """Human-readable name for a PETSc KSP converged-reason code."""
+    return KSP_REASON_STRINGS.get(int(reason), f"KSP_UNKNOWN_{reason}")
+
+
 def contraction(history) -> Optional[float]:
     """Geometric-mean per-iteration contraction factor ρ from a residual ladder.
 

src/underworld3/utilities/rotated_bc.py

@@ -277,7 +277,8 @@ def solve_rotated_freeslip(solver, boundaries, remove_rotation_gauge=True, verbo
         * ``"normal_rows"`` — global rows of the constrained normal components;
         * ``"boundaries"`` — the boundary specs the rotation was built from;
         * ``"rotation_gauge_removed"`` — whether a rigid-rotation gauge was projected out;
-        * ``"ksp_reason"``, ``"ksp_its"`` — outer KSP converged-reason and iteration count.
+        * ``"ksp_reason"``, ``"ksp_its"`` — outer KSP converged-reason and iteration count;
+        * ``"rnorm"`` — true rotated residual ‖·û − b̂‖ (feeds the solve report).
     """
     if getattr(solver, "snes", None) is None:
         solver._setup_pointwise_functions()
@@ -369,6 +370,15 @@ def solve_rotated_freeslip(solver, boundaries, remove_rotation_gauge=True, verbo
         ksp_its = ctx["ksp"].getIterationNumber()
         _destroy_rotated_ksp_ctx(ctx)
 
+    # True rotated residual ‖·û − b̂‖ for the solve report (one matvec). Computed
+    # explicitly rather than read off the KSP: preonly/LU never computes a norm, and
+    # the iterative norm can be the preconditioned one.
+    _res = bhat.duplicate()
+    Ahat.mult(Uhat, _res)
+    _res.axpy(-1.0, bhat)
+    rnorm = float(_res.norm())
+    _res.destroy()
+
     # rotate back u = Qᵀ û  (U is returned in the result dict → create, don't
     # borrow from the pool)
     U = dm.createGlobalVec()
@@ -379,7 +389,7 @@ def solve_rotated_freeslip(solver, boundaries, remove_rotation_gauge=True, verbo
     return {"Q": Q, "Qt": Qt, "A": Aorig, "b": b, "U": U, "Uhat": Uhat,
             "normal_rows": normal_rows, "boundaries": list(boundaries),
             "rotation_gauge_removed": removed, "ksp_reason": ksp_reason,
-            "ksp_its": ksp_its}
+            "ksp_its": ksp_its, "rnorm": rnorm}
 
 
 def _finalize_rotated_solution(solver, U, Q, normal_rows, remove_rotation_gauge):
@@ -558,6 +568,8 @@ def solve_rotated_freeslip_nonlinear(solver, boundaries, remove_rotation_gauge=T
         * ``"ksp_its"`` — list of linear iteration counts, one per Newton iteration;
         * ``"nonlinear_iterations"``, ``"converged"`` — outer-loop count (the number
           of Newton increments solved, ``== len(ksp_its)``) and status;
+        * ``"rnorm"``, ``"rnorm0"`` — final and initial rotated residual norms ‖F̂‖
+          (feed the solve report);
         * ``"continuation_switched"`` — whether the Picard→Newton tangent switch fired.
     """
     if getattr(solver, "snes", None) is None:
@@ -724,7 +736,7 @@ def rotated_residual(uvec, keep_cartesian=False):
             "normal_rows": normal_rows, "boundaries": list(boundaries),
             "rotation_gauge_removed": removed, "ksp_reason": last_reason,
             "nonlinear_iterations": newton_its, "converged": converged,
-            "ksp_its": lin_its,
+            "ksp_its": lin_its, "rnorm": rnorm, "rnorm0": r0,
             "continuation_switched": continuation and phase == "newton"}
 
 

tests/test_1055_solve_report.py

@@ -150,3 +150,114 @@ def test_solve_report_helpers():
     assert reason_string(999).startswith("UNKNOWN")
     assert contraction([50.0, 1e-3, 1e-6, 1e-9]) is not None
     assert contraction([50.0]) is None        # < 2 points -> undefined
+
+
+# --------------------------------------------------------------------------- #
+#  Rotated free-slip paths (solve outside self.snes -> report from the rotated
+#  result dicts; regression for the review findings on PR #377)
+# --------------------------------------------------------------------------- #
+
+def _rotated_stokes(nonlinear=False):
+    """Box Stokes with rotated strong free-slip on all four walls."""
+    mesh = uw.meshing.UnstructuredSimplexBox(
+        minCoords=(0, 0), maxCoords=(1, 1), cellSize=0.25
+    )
+    v = uw.discretisation.MeshVariable("Ur", mesh, 2, degree=2)
+    p = uw.discretisation.MeshVariable("Pr", mesh, 1, degree=1, continuous=True)
+    st = uw.systems.Stokes(mesh, velocityField=v, pressureField=p)
+    st.constitutive_model = uw.constitutive_models.ViscousFlowModel
+    if nonlinear:
+        edot = st.Unknowns.Einv2 + 1.0e-8
+        st.constitutive_model.Parameters.shear_viscosity_0 = (1.0e-2 + edot) ** (-0.4)
+    else:
+        st.constitutive_model.Parameters.shear_viscosity_0 = 1.0
+    st.tolerance = 1.0e-6
+    x, y = mesh.X
+    st.bodyforce = sympy.Matrix([0, sympy.sin(np.pi * x) * sympy.sin(np.pi * y)])
+    for w in ("Bottom", "Top", "Left", "Right"):
+        st.add_rotated_freeslip_bc(0, w)
+    st.petsc_use_pressure_nullspace = True
+    return st, v
+
+
+@pytest.mark.level_1
+@pytest.mark.tier_a
+def test_rotated_linear_solve_report():
+    """The linear rotated path leaves a populated report: KSP-namespace reason,
+    a real residual norm (not NaN), and one entry in the history trail."""
+    st, _ = _rotated_stokes(nonlinear=False)
+    st.solve()
+    r = st.solve_report
+    assert isinstance(r, SolveReport)
+    assert r.converged and r.ksp_its >= 1
+    assert r.nl_its == 1                       # one outer solve on the linear path
+    assert r.fnorm == r.fnorm                  # finite: rnorm is in the result dict
+    assert r.reason_str.startswith("KSP_")     # KSP namespace, not the SNES table
+    assert len(st.solve_history) == 1
+
+
+@pytest.mark.level_1
+@pytest.mark.tier_a
+def test_rotated_nonlinear_solve_report():
+    """The manual nonlinear rotated loop reports its OUTER verdict and effort:
+    nl_its from nonlinear_iterations, ksp_its summed over the per-Newton list."""
+    st, _ = _rotated_stokes(nonlinear=True)
+    st.solve()
+    info = st._rotated_freeslip_info
+    r = st.solve_report
+    assert isinstance(r, SolveReport), "nonlinear rotated solve left no report"
+    assert r.converged == info["converged"]
+    assert r.nl_its == info["nonlinear_iterations"] and r.nl_its >= 2
+    assert r.ksp_its == sum(info["ksp_its"])
+    assert r.fnorm == r.fnorm and r.fnorm0 is not None and r.fnorm < r.fnorm0
+    assert r.reason_str.startswith("KSP_")
+    assert len(st.solve_history) == 1
+
+
+@pytest.mark.level_1
+@pytest.mark.tier_a
+def test_estimate_difficulty_rejects_rotated_and_bad_kwargs():
+    st, _ = _rotated_stokes(nonlinear=False)
+    with pytest.raises(NotImplementedError):
+        st.estimate_difficulty(max_nl_its=2)
+    st2, _ = _nonlinear_stokes()
+    with pytest.raises(TypeError):
+        st2.estimate_difficulty(max_nl_its=2, zero_init_guess=True)
+    with pytest.raises(TypeError):
+        st2.estimate_difficulty(max_nl_its=2, divergence_retries=1)
+
+
+@pytest.mark.level_1
+@pytest.mark.tier_a
+def test_warm_probe_on_converged_state_does_not_poison_anchor():
+    """A warm estimate_difficulty on an already-converged state must not arm an
+    anchor at tolerance*||F_converged|| (an unreachable target that would make the
+    next probe grind and mislabel its exit)."""
+    st, _ = _nonlinear_stokes()
+    st.solve(zero_init_guess=True)
+    rep = st.estimate_difficulty(max_nl_its=10, warm=True)
+    assert rep.converged
+    assert st._resume_abs_target is None       # converged probe leaves no armed anchor
+    rep2 = st.estimate_difficulty(max_nl_its=10, warm=True)
+    assert rep2.converged and rep2.nl_its <= 1  # no grinding on a converged state
+    assert st._resume_abs_target is None
+
+
+@pytest.mark.level_1
+@pytest.mark.tier_a
+def test_ksp_reason_table_matches_petsc():
+    """The KSP reason table is hand-written (this module must import without petsc4py
+    at build-doc time) — pin it to the real enum so a PETSc renumbering is caught."""
+    from petsc4py import PETSc
+    from underworld3.systems.solve_report import KSP_REASON_STRINGS, ksp_reason_string
+
+    enum_names = {}
+    for name, value in vars(PETSc.KSP.ConvergedReason).items():
+        if isinstance(value, int) and not name.startswith("_"):
+            enum_names.setdefault(value, name)
+    for code, label in KSP_REASON_STRINGS.items():
+        if code == 0:
+            continue                           # petsc4py spells 0 CONVERGED_ITERATING
+        assert label == f"KSP_{enum_names[code]}", (code, label, enum_names[code])
+    assert ksp_reason_string(-3) == "KSP_DIVERGED_MAX_IT"
+    assert ksp_reason_string(999).startswith("KSP_UNKNOWN")