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