underworldcode · lmoresi · Aug 23, 2026 · Aug 22, 2026
@@ -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.

@@ -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

@@ -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