From 6c4b0a4ae069bb79bdbe91b4630e7c99ddd3c879 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Sat, 22 Aug 2026 20:44:47 +1000 Subject: [PATCH] Count sub-solves that stop at their iteration cap, and record the FMG pairing PETSc's KSPCheckSolve deliberately does not treat DIVERGED_MAX_IT on a sub-KSP as a failure, so a block can run out of iterations on every application and the outer solve still reports CONVERGED. That silence is what #625 is made of. It is expensive because S = -B A^-1 B^T is applied THROUGH the velocity solve. Truncate those and the pressure Krylov is chasing an operator that moves between applications, so it cannot converge either and caps in turn. Measured on SolCx at h=1/30, changing only fieldsplit_velocity_ksp_max_it from 200 to 5000: 976 s to 25.6 s, outer 44 to 2, pressure 200/application to 30, and an identical velocity error. Thirty-eight times, from a number that is not in the pressure block. solve_report.sub[...].capped now counts the applications that ended that way, read from each application's converged reason via a post-solve hook rather than by comparing an iteration count against the cap -- a solve that converges on its last permitted iteration is converged, and counting would misreport it. The skill records the pairing this exposed: with FMG, penalty=10 improves every axis at once (21% faster, Schur count 3.3x smaller, total velocity work halved), because FMG absorbs grad-div augmentation where GAMG is driven into its cap by it. With GAMG the same penalty makes the solve slower, so the two must not be recommended independently. The accuracy cost is a consistent perturbation -- same convergence rate, gap shrinking 1.102 -> 1.087 -> 1.066 under refinement -- which is what makes the pairing safe to recommend at all. No numerical behaviour changes here: the counter is observation, and the defaults are untouched. Underworld development team with AI support from Claude Code --- .claude/skills/nonlinear-solver/SKILL.md | 78 ++++++++++++++++++++ src/underworld3/systems/solver_health.py | 33 ++++++++- tests/test_0205_subsolve_cap_detection.py | 90 +++++++++++++++++++++++ 3 files changed, 200 insertions(+), 1 deletion(-) create mode 100644 tests/test_0205_subsolve_cap_detection.py diff --git a/.claude/skills/nonlinear-solver/SKILL.md b/.claude/skills/nonlinear-solver/SKILL.md index 33b70a5a..04ea6a27 100644 --- a/.claude/skills/nonlinear-solver/SKILL.md +++ b/.claude/skills/nonlinear-solver/SKILL.md @@ -221,6 +221,84 @@ Companion skills: **`adapt-on-top-faults`** (building the child, engines, repair band sizing), **`adaptive-meshing`** (the mover, and `relax(pin_bands=...)` for relaxing a mesh that was refined onto an interface). +## The Schur complement: pair the penalty with FMG, never with GAMG + +**Symptom this is for**: the velocity block's iteration count is rock solid but +the pressure sub-solve wanders into the hundreds and eventually stops +converging. + +**First: it is probably not the pressure block.** `S = -B A^-1 B^T` is applied +*through* the velocity solve, so a velocity solve that exits at its iteration +cap makes the Schur operator inconsistent between applications — and no Krylov +method converges against an operator that moves under it. The pressure block +then caps too, and the outer flounders. Measured on SolCx (eta 1e6, P2-P0disc, +h=1/30), changing **only** `fieldsplit_velocity_ksp_max_it`: + +| velocity cap | sec | outer | pressure/app | velocity/app | +|---|---|---|---|---| +| 200 (default) | 976.0 | 44 | **200.0** | **200.0** | +| 5000 | **25.6** | **2** | **30.0** | 618.0 | + +**38x from a number that is not in the pressure block**, and the velocity error +is identical in both rows. Before tuning the Schur solve, check whether either +block sat at exactly its cap — `solve_report.sub` gives iterations and +applications per block, and a per-application count equal to the cap to the +digit is the tell. + +**Then: the penalty is the lever on the Schur count, and it needs FMG.** +`stokes.penalty = lambda` adds `lambda*mu*(div u)(div v)`, which makes the +eta-scaled mass matrix a better approximation to S. Matched on one mesh +(2592 cells), same discrete solve, only the velocity preconditioner differs: + +| lambda | velocity PC | sec | outer | Schur/app | velocity/app | velocity total | +|---|---|---|---|---|---|---| +| 0 | GAMG | 15.49 | 2 | 125.5 | 94.7 | 24802 | +| 0 | **FMG** | **3.88** | 1 | **59.0** | **8.8** | **546** | +| 10 | GAMG | 20.68 | 7 | 22.3 | **199.9 capped** | 33976 | +| 10 | **FMG** | **3.06** | 1 | **18.0** | **13.5** | **270** | + +- **With FMG, `penalty = 10` improves every axis at once**: 21% faster, Schur + count 3.3x smaller, total velocity work halved. FMG absorbs grad-div + augmentation (8.8 -> 13.5 iterations per application); GAMG does not + (94.7 -> capped). +- **With GAMG, do not use it at all.** The same `penalty = 10` makes the solve + *slower* (15.49 -> 20.68 s), because augmentation is exactly what drives GAMG + into its cap. Uncapping rescues it to 11.03 s but it still needs **833** + iterations per application, and FMG is 3.6x faster on the same mesh. + Feasible is not competitive. + +**The accuracy cost is consistent, so it is safe to pair by default.** The +penalty is grad-div, not a true augmented Lagrangian — `div(P2)` is not inside +`P0`, so the term does not vanish at the discrete solution and it does perturb +the answer. But the perturbation converges away: same rate, and the gap shrinks +under refinement. + +| cells | lambda=0 v err | rate | lambda=10 v err | rate | gap | +|---|---|---|---|---|---| +| 648 | 2.112e-1 | — | 2.327e-1 | — | 1.102 | +| 2592 | 1.266e-1 | 1.67 | 1.376e-1 | 1.69 | 1.087 | +| 10368 | 8.727e-2 | 1.45 | 9.305e-2 | 1.48 | **1.066** | + +For a pressure-dependent constitutive law use the mechanical pressure, +`p_mech = p - lambda*mu*(div u)`; the raw `p` is the multiplier. + +**Traps.** + +- **FMG needs a refined base or you silently get GAMG.** Measured: + `refinement=0` -> one hierarchy level -> default velocity PC is `gamg`; + `refinement=2` -> `mg`. So `penalty` set "with FMG" on an unrefined mesh is + actually the harmful GAMG pairing. Check + `snes.getKSP().getPC().getFieldSplitSubKSP()[0].getPC().getType()`, or read + `solver.pc_fallbacks`. +- **Scaling `saddle_preconditioner` by a constant does nothing** — it does not + change the Krylov subspace. `1/eta` and `101/eta` both give 28 iterations, + identical to every digit, so an "AL-matched" `1/(eta*(1+lambda))` cannot help. + The 1/eta *weighting* itself is worth 1.9x (28 vs 52 with a flat `1`). +- **Eisenstat-Walker is inert under `snes_type=ksponly`** — identical iterations + and error on or off. And `outer 1` is not an EW artefact: it is what a full + Schur factorisation gives when the Schur complement is solved well. +- Measurements: `~/+Simulations/pressure_schur_625/` (#625). + ## Gotchas - **`./uw build` → `amr-dev` env**; verify `uw.__file__` is the worktree site-packages. diff --git a/src/underworld3/systems/solver_health.py b/src/underworld3/systems/solver_health.py index 9c946615..577b3596 100644 --- a/src/underworld3/systems/solver_health.py +++ b/src/underworld3/systems/solver_health.py @@ -60,6 +60,10 @@ # Marking the sub-preconditioner failed is the whole point -- see the module docstring. _DEADLINE_REASON = PETSc.KSP.ConvergedReason.DIVERGED_BREAKDOWN +#: A solve that stopped because it ran out of iterations, not because it +#: converged. PETSc's KSPCheckSolve deliberately does NOT treat this as a +#: failure for a sub-KSP, which is exactly why it goes unnoticed. +_CAPPED_REASON = PETSc.KSP.ConvergedReason.DIVERGED_MAX_IT # KSP types whose behaviour changes when a monitor is attached (both are gated on # ksp->numbermonitors in PETSc): preonly.c computes norms it would otherwise skip, and @@ -85,6 +89,18 @@ class SubSolveReport: work axis. applications Number of times the sub-solve ran during the outer solve. + capped + How many of those applications ended at the block's iteration cap + (``KSP_DIVERGED_MAX_IT``) rather than at its tolerance. **Non-zero means + the block did not solve.** It matters far more than it looks: the Schur + operator ``S = -B A^-1 B^T`` is applied *through* the velocity solve, so + a capped velocity block hands the pressure Krylov an operator that moves + between applications, and no Krylov method converges against that. Every + expensive pathology measured in #625 was a cap — the pressure block, the + velocity block beneath it, and the velocity block again under + augmentation. Measured: raising only ``fieldsplit_velocity_ksp_max_it`` + from 200 to 5000 took one solve from 976 s to 25.6 s with an identical + answer. complete ``False`` when the block was instrumented *during* this solve rather than before it, which happens on a solver's very first solve: the fieldsplit blocks @@ -99,9 +115,11 @@ class SubSolveReport: its: int applications: int complete: bool = True + capped: int = 0 def __str__(self) -> str: return (f"{self.name}: {self.its} its / {self.applications} applications" + + (f" — {self.capped} AT THE ITERATION CAP" if self.capped else "") + ("" if self.complete else " (lower bound)")) @@ -129,20 +147,33 @@ def __init__(self, ksp, name, instrumentation, is_outer, mid_solve): self.name = name self.its = 0 self.applications = 0 + self.capped = 0 self.complete = not mid_solve self._instrumentation = instrumentation self.is_outer = is_outer self._test_installed = False ksp.setMonitor(self._monitor) + # The REASON, not an iteration count compared against the cap: a solve + # that happens to converge on its last permitted iteration is converged, + # and counting would call it capped. A post-solve hook is the only place + # the per-application reason can be read before the next solve resets it. + ksp.setPostSolve(self._after_application) def reset(self): self.its = 0 self.applications = 0 + self.capped = 0 self.complete = True # attached before this solve began def report(self): return SubSolveReport(name=self.name, its=self.its, - applications=self.applications, complete=self.complete) + applications=self.applications, + complete=self.complete, capped=self.capped) + + def _after_application(self, ksp, rhs, x): + """Record whether this application ran out of iterations.""" + if ksp.getConvergedReason() == _CAPPED_REASON: + self.capped += 1 def _monitor(self, ksp, its, rnorm): # PETSc calls the monitor once per iteration including iteration 0, which diff --git a/tests/test_0205_subsolve_cap_detection.py b/tests/test_0205_subsolve_cap_detection.py new file mode 100644 index 00000000..d33735c3 --- /dev/null +++ b/tests/test_0205_subsolve_cap_detection.py @@ -0,0 +1,90 @@ +"""A sub-solve that stops at its iteration cap must say so. + +PETSc's ``KSPCheckSolve`` deliberately does not treat ``DIVERGED_MAX_IT`` on a +sub-KSP as a failure — truncating an inner solve is normal — so a block can run +out of iterations on every application and the outer solve still reports +CONVERGED. That silence is expensive: ``S = -B A^-1 B^T`` is applied *through* +the velocity solve, so once those are truncated the pressure Krylov is chasing +an operator that moves between applications and cannot converge either. +Measured on SolCx at h=1/30 (#625): raising only +``fieldsplit_velocity_ksp_max_it`` from 200 to 5000 took the solve from 976 s to +25.6 s, with an identical answer. + +``solve_report.sub[...].capped`` counts the applications that ended that way. +""" + +import pytest + +import underworld3 as uw +from underworld3 import analytic as A + +pytestmark = [pytest.mark.level_1, pytest.mark.tier_a] + + +def _solcx_stokes(velocity_cap): + """SolCx with a viscosity contrast hard enough to need real iterations.""" + mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=1.0 / 10, qdegree=3 + ) + sol = A.SolCx(mesh, eta_A=1.0, eta_B=1.0e6, x_c=0.5, n=1) + v = uw.discretisation.MeshVariable("Ucap", mesh, mesh.dim, degree=2) + p = uw.discretisation.MeshVariable("Pcap", mesh, 1, degree=0, continuous=False) + + stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel + stokes.constitutive_model.Parameters.shear_viscosity_0 = sol.fn_viscosity + stokes.saddle_preconditioner = 1.0 / sol.fn_viscosity + stokes.bodyforce = sol.fn_bodyforce + for wall, condition in (("Left", (0.0, None)), ("Right", (0.0, None)), + ("Bottom", (None, 0.0)), ("Top", (None, 0.0))): + stokes.add_dirichlet_bc(condition, wall) + stokes.petsc_use_pressure_nullspace = True + stokes.petsc_options["snes_type"] = "ksponly" + stokes.tolerance = 1.0e-6 + stokes.petsc_options["fieldsplit_velocity_ksp_max_it"] = velocity_cap + stokes.solve() + return stokes + + +def _blocks(stokes): + return {entry.name: entry for entry in stokes.solve_report.sub.values()} + + +def test_a_capped_velocity_block_is_reported(): + """Cap the velocity block hard; every application must be counted.""" + blocks = _blocks(_solcx_stokes(velocity_cap=2)) + velocity = blocks["velocity"] + + assert velocity.applications > 0 + assert velocity.capped == velocity.applications, ( + f"{velocity.capped} of {velocity.applications} velocity applications " + "were reported as capped; a 2-iteration cap on a 1e6 viscosity " + "contrast cannot converge any of them" + ) + assert "AT THE ITERATION CAP" in str(velocity) + + +def test_an_untruncated_solve_reports_no_caps(): + """Negative control. Without this, a counter stuck at `applications` passes. + + The same problem with a cap it never reaches must come back clean, or the + check is reporting the cap rather than detecting it. + """ + for block in _blocks(_solcx_stokes(velocity_cap=5000)).values(): + assert block.capped == 0, ( + f"{block.name} reported {block.capped} capped application(s) with a " + f"5000-iteration cap it never approached: {block}" + ) + assert "AT THE ITERATION CAP" not in str(block) + + +def test_the_count_distinguishes_some_from_all(): + """`capped` counts applications, not a boolean. + + A block that truncates on some applications and converges on others is the + interesting middle case — reporting only "did it ever cap" would lose the + difference between an occasional truncation and a block that never solves. + """ + velocity = _blocks(_solcx_stokes(velocity_cap=2))["velocity"] + assert isinstance(velocity.capped, int) + assert 0 < velocity.capped <= velocity.applications