diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index 9aa2e32d..0fda2ebd 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/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. diff --git a/src/underworld3/systems/solve_report.py b/src/underworld3/systems/solve_report.py index 95b56efe..8cdd6358 100644 --- a/src/underworld3/systems/solve_report.py +++ b/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. diff --git a/src/underworld3/utilities/rotated_bc.py b/src/underworld3/utilities/rotated_bc.py index 1bb9c9c9..0bff8091 100644 --- a/src/underworld3/utilities/rotated_bc.py +++ b/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"} diff --git a/tests/test_1055_solve_report.py b/tests/test_1055_solve_report.py index 85c549ce..be3de21d 100644 --- a/tests/test_1055_solve_report.py +++ b/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")