One owner for the geometric-MG option bundle; rotated free-slip picks up a mesh-owned hierarchy (#468, #467)

docs/developer/CHANGELOG.md

@@ -6,6 +6,49 @@ This log tracks significant development work at a conceptual level, suitable for
 
 ## 2026 Q3 (July – September)
 
+### One Owner for the Geometric-Multigrid Option Bundle (July 2026)
+
+**The PETSc option bundle that configures a Stokes velocity block's multigrid
+now lives in exactly one module, and all three routes that reach that block read
+it from there** (#468), **and rotated free-slip now picks up a mesh-owned
+multigrid hierarchy instead of silently discarding it** (#467).
+
+Three routes reach a multigrid velocity block: native (PETSc interpolation
+between refined DMPlex levels), custom-P on the standard solve path, and custom-P
+through the rotated free-slip path. They are the same preconditioner reached
+three ways, not alternatives — custom-P is *mandatory* wherever native cannot go,
+namely rotated boundary conditions and `adapt()` children. The bundle was
+written in two places and had drifted: the native path had been moved to a
+`gmres`+`sor` smoother on a recorded measurement, and the custom-P routes had
+not. Worse, the custom-P writer never *set* the smoother iteration count at all,
+so it inherited whatever had last written that options prefix — 3 left behind by
+the GAMG bundle on the standard path, PETSc's own default of 2 on the rotated
+path. The same function smoothed differently depending on what had run before
+it.
+
+Unifying the bundle recovers, on the same operator, right-hand side and coarse
+solve: rotated custom-P velocity-block iterations 11 → 5 (0.68 s → 0.39 s of
+linear solve, timed in isolation), standard custom-P 5 → 4, on a *two-level*
+hierarchy — the depth at which the native measurement says the gmres margin is
+smallest. The bundle also now derives which stale keys it must clear rather than
+carrying a hand-maintained list, which is what let the iteration count go unset
+in the first place.
+
+Separately, `mesh.adapt()` leaves a coarse tail on its refinement child so that
+every solver on an adapted mesh gets geometric multigrid with no per-solver call.
+The rotated path never consulted it — the standard path's injection hook runs
+after the rotated dispatch has already returned — so an adapt child under rotated
+free-slip fell back to algebraic multigrid, indistinguishable from having no
+hierarchy at all. That is the `adapt-on-top-faults` workflow's own configuration
+(a fault resolved by local refinement, with rotated free-slip chosen because it
+composes with transverse isotropy). Both paths now resolve the hierarchy through
+one shared rule, with the same opportunistic degrade-to-GAMG behaviour.
+
+The regression test reads the smoother configuration back off the **live PETSc
+objects** for all three routes and asserts they agree. An options-database
+assertion would not have caught the original drift, because the drift was
+precisely a key nobody wrote.
+
 ### One Rotated Free-Slip Path, Now With a Prescribed Wall-Normal Velocity (July 2026)
 
 **Rotated strong free-slip now takes a prescribed wall-normal velocity datum

docs/developer/subsystems/rotated-freeslip.md

@@ -41,11 +41,42 @@ There is a single driver, `solve_rotated_freeslip`: a manual outer
 Newton/Picard loop that rotates the residual and tangent every iteration
 (`F̂ = Q F`, `Ĵ = Q J Qᵀ`), imposes the constraint on the rotated normal rows,
 and solves each increment with a self-contained fieldsplit-Schur KSP (geometric
-FMG on the custom prolongation when a hierarchy is registered, else GAMG; the
+FMG on the custom prolongation when the solver has a hierarchy, else GAMG; the
 native 1/μ pressure mass as the Schur preconditioner). There is no separate
 linear path and no up-front nonlinearity probe: a linear model converges after
 its first increment and the loop self-terminates.
 
+### Where the multigrid hierarchy comes from
+
+Native geometric FMG cannot serve a rotated solve at all — the DM-coupled
+hierarchy has no way to express the per-node rotation — so custom-P multigrid is
+not an optimisation here, it is the only geometric option. The hierarchy is
+resolved by `custom_mg.build_transfers`, which is the **same rule the standard
+path uses**:
+
+1. an explicit `set_custom_fmg` registration on the solver wins; otherwise
+2. a **mesh-owned** coarse tail (`mesh._custom_mg_coarse_meshes`, what
+   `mesh.adapt()` leaves on a refinement child) is picked up opportunistically,
+   with the barycentric→RBF builder fallback, so a failed build degrades to GAMG
+   rather than crashing the solve.
+
+Step 2 used to be unreachable from here. The standard path's injection hook runs
+*after* the rotated dispatch has already returned, and the rotated builder
+consulted only `solver._custom_mg` — so an `adapt()` child under rotated
+free-slip reported `velocity_pc == "GAMG"`, indistinguishable from having no
+hierarchy at all, and lost its multigrid silently (#467). That is the
+`adapt-on-top-faults` workflow's own configuration: a fault resolved by local
+refinement on an adapt child, with rotated free-slip chosen over Nitsche because
+it composes with transverse isotropy.
+
+The option bundle for both the FMG and the GAMG velocity block comes from
+`utilities/multigrid_options.py`, shared with the native and standard custom-P
+routes, so the rotated block cannot be configured differently by accident. The
+single deliberate difference is the coarse solve: `coarse="svd"`, because the
+Galerkin-coarsened rotated block inherits the rigid-rotation null space where
+`redundant`/LU hits a zero pivot. See
+[the solvers subsystem doc](solvers.md) for the bundle itself.
+
 ### The velocity sub-KSP must converge, not just apply
 
 The velocity block is preconditioned by multigrid, but the KSP *around* that

docs/developer/subsystems/solvers.md

@@ -175,7 +175,7 @@ def validate_unknowns_sharing(multi_material_model):
             assert model.Unknowns.DFDt is reference_unknowns.DFDt, \
                 f"Model {i} $D\mathbf{{F}}/Dt$ not shared - stress history will be wrong"
 ```
-- ⚠️ Preconditioner *selection* is partly covered — see "Choosing the Krylov method for a fieldsplit sub-solve" below; multigrid and Schur preconditioner choice is still undocumented
+- ⚠️ Preconditioner *selection* is partly covered — see "Choosing the Krylov method for a fieldsplit sub-solve" and "The multigrid option bundle, and its one owner" below; Schur preconditioner choice is still undocumented
 - Could benefit from optimization examples
 
 ## Choosing the Krylov method for a fieldsplit sub-solve
@@ -327,6 +327,193 @@ application, and the velocity KSP is applied once per Schur `MatMult` — i.e. o
 pressure Krylov iteration. That number is a sample, not work. `SolverInstrumentation`
 (`systems/solver_health.py`) is the mechanism that sums properly.
 
+## The multigrid option bundle, and its one owner
+
+Three routes reach a multigrid velocity block, and they are **the same
+preconditioner reached three ways**, not three alternatives:
+
+| route | when | prolongation |
+|---|---|---|
+| native | mesh built with `refinement >= 1`, ordinary BCs | PETSc `DMCreateInterpolation` between refined DMPlex levels |
+| custom-P, standard path | `set_custom_fmg`, or an `adapt()` child's mesh-owned coarse tail | barycentric / RBF, Galerkin coarse operators |
+| custom-P, rotated path | rotated free-slip, via `rotated_bc` | as above, with the fine prolongation rotated, `P̂ = Q_v·P` |
+
+custom-P is **mandatory** wherever native cannot go: rotated boundary conditions
+(the DM-coupled hierarchy cannot express a per-node rotation) and non-nested
+grids (`adapt()` children have no DMPlex refinement relation). So the routes are
+not ranked — the one that matters most for adapted and curved-boundary work is
+custom-P.
+
+The option *values* live in one module, `utilities/multigrid_options.py`. Every
+writer reads a bundle from there and applies it to its own options object under
+its own prefix; nobody writes a multigrid option value anywhere else. That is
+structural rather than stylistic: when the bundle was written in two places, the
+native path was deliberately moved to a measured `gmres`+`sor` smoother and the
+custom-P path was not, and the custom-P routes ran `richardson` at an iteration
+count **nobody had set** — inherited from whatever last wrote that options
+prefix (3 left behind by the GAMG bundle on the standard path, PETSc's own PCMG
+default of 2 on the rotated path). The same function smoothed differently
+depending on what had run before it.
+
+### What a bundle carries
+
+A bundle is the settings it sets *and* the keys it must clear. The clear-list is
+derived, not hand-written: it is every key any sibling bundle sets that this one
+does not. These bundles share an options prefix, so switching a block from GAMG
+to geometric MG leaves the GAMG-only keys behind and `setFromOptions` will
+happily re-read them. Deriving the list means a key added to one bundle
+automatically becomes stale for the others.
+
+Two consequences worth stating outright:
+
+- **Set the smoother iteration count, never inherit it.** A bundle that omits
+  `mg_levels_ksp_max_it` is not "using the default" — it is using whatever the
+  last writer left.
+- **The measured smoother is `gmres`+`sor` at 4 iterations.** Chebyshev needs
+  eigenvalue estimates of the smoothed operator, which are fragile on the
+  indefinite / variable-viscosity velocity block. Against richardson, measured on
+  the Spiegelman notch (Drucker–Prager, η contrast 1e26, nested 4-level
+  hierarchy): per-cycle contraction ρ = 0.75 richardson against 0.56 gmres at the
+  *same* four iterations, the margin growing with depth (5% at 3 levels, 25% at 4).
+
+### Why richardson loses — the recorded reason is wrong
+
+Both this bundle's comment and the Layer 3 section of
+`design/nonlinear-solver-homotopy-warmstart.md` said richardson degrades because
+the consistent-Newton tangent makes the velocity block **non-symmetric**. That
+mechanism does not survive measurement.
+
+For an isotropic η(ε̇_II) the Newton term is
+
+$$2\,\frac{\eta'}{\dot\varepsilon_{II}}\big(\dot\varepsilon(u):\dot\varepsilon(\delta u)\big)\big(\dot\varepsilon(u):\dot\varepsilon(v)\big)$$
+
+— a rank-one outer product $a\otimes a$ with $a = \dot\varepsilon(u)$, which is
+symmetric. η depends on ∇v only *through* its symmetric part, so both factors
+project onto the same tensor and there is nothing left to break the symmetry.
+Measured as ‖A−Aᵀ‖_F/‖A‖_F on the assembled velocity block, against a linear
+calibration of 5e-17:
+
+| rheology | ‖A−Aᵀ‖/‖A‖ | Newton term ‖A_N−A_P‖/‖A‖ |
+|---|---|---|
+| linear (calibration) | 5.1e-17 | — |
+| shear-thinning, isotropic | 5.3e-17 | 2.1e-01 |
+| power-law, isotropic | 5.3e-17 | 4.0e-01 |
+| Drucker–Prager, pressure-dependent yield | 5.1e-17 | 5.1e-02 |
+| transverse isotropic, Picard tangent | **7.2e-02** | — |
+| transverse isotropic, Newton tangent | **8.9e-02** | 2.0e-01 |
+
+The third column is the control: if the consistent tangent had contributed
+nothing, the symmetry reading would be vacuous. It contributed 5–40% of ‖A‖ in
+every case, and the block stayed symmetric anyway.
+
+Non-symmetry *does* appear under transverse isotropy — but it appears under the
+**Picard** tangent too, where a frozen-coefficient form ∫ ε̇(δu):C:ε̇(v) is
+symmetric by construction. That is a defect signature rather than a property of
+the tangent (see #457, and the linear TI rows, which are symmetric at 6.3e-17
+even with anisotropy fully active).
+
+The mechanism that does fit everything measured is **operator conditioning**: at
+η contrast 1e26 the SOR-preconditioned spectrum is spread far enough that a
+stationary iteration stalls where a Krylov smoother adapts its polynomial. It
+predicts the notch result, and it predicts the cost below on well-conditioned
+problems.
+
+### Two regimes, and the named intent that selects them
+
+Neither smoother dominates, so there are two variants and `solver.strategy`
+chooses between them. This is the layering:
+
+| layer | what it is | where |
+|---|---|---|
+| `solver.strategy` | the **named intent** — `"default"` / `"robust"` / `"fast"` | public property on the solver |
+| `MGSettings` | the **option values** an intent resolves to | `utilities/multigrid_options.py` |
+| `solver.petsc_options` | the **escape hatch** — any key you set here outranks both | public, unchanged |
+
+`"robust"` (gmres/4) survives an operator a stationary smoother stalls on: the
+Spiegelman notch contraction above, and 11 → 5 velocity iterations on the
+transversely isotropic rotated annulus.
+
+`"fast"` (richardson/3) is cheaper per cycle and quicker where the operator is
+benign. Nested annulus, **linear (symmetric)** velocity block, η contrast 1e6,
+outer KSP timed in isolation — `"fast"` against `"robust"`:
+
+| depth | iterations | wall clock |
+|---|---|---|
+| 2 levels | ×0.69 | **×1.16** |
+| 3 levels | ×0.63 | **×1.30** |
+| 4 levels | ×0.83 | **×1.82** |
+
+The extra sweep costs ~3%; gmres itself costs ~40% per cycle. Custom-P builds its
+Galerkin coarse operators from barycentric transfers, denser than the native
+nested ones, so the smoother is a larger share of the cycle on that route.
+
+`"default"` is `"robust"`: the failure it avoids is worse than the cost it carries,
+and it carries that cost exactly where the problem is easy. `"fast"` is the
+documented opt-out.
+
+### Asking what you actually got
+
+`solver.strategy` reports it. The value is still the strategy *name* — it compares
+and formats as the plain string, so nothing that used it before changes — but
+displaying it shows the preconditioner that name resolved to:
+
+```python
+>>> stokes.strategy
+'default' — geometric multigrid (2 levels), full cycle, smoother gmresx4 + sor,
+            coarse redundant/lu
+>>> stokes.strategy == "default"
+True
+```
+
+Three properties of that report matter more than its formatting:
+
+- **It names a user override rather than absorbing it.** A key you set yourself is
+  listed as such (`; overridden by the user: mg_levels_ksp_max_it=6`), so the
+  summary cannot hide the difference between what the strategy asked for and what
+  is running.
+- **It refuses to report before it knows.** The preconditioner is resolved at the
+  first solve; until then the report says so instead of presenting the constructor
+  defaults as though they were the answer. A summary that looks authoritative and
+  is stale is worse than no summary.
+- **`solver.preconditioner_settings`** is the same information as a dict, so a test
+  can assert on it rather than parsing prose or inferring from timings.
+
+This is the general shape #484 asks for across all the managed fallbacks: report the
+resolved choice as readable state, and distinguish "could not" from "chose not to".
+
+```{warning}
+Do not fill a strategy by writing velocity-block options from the strategy setter.
+`_apply_preconditioner_options` runs later, at `_build`, and is the single writer
+of that block — anything written earlier is overwritten. That is how
+`pc_mg_type=kaskade` sat in the `strategy` setter for a long time carrying a
+comment warning against changing it, while never once taking effect (measured: the
+live PC was `mg`/FULL in every ordering). A strategy selects a bundle *variant*.
+```
+
+### The one legitimate per-route difference
+
+The coarse solve. The Galerkin-coarsened **rotated** velocity block inherits the
+rigid-rotation null space of the constrained problem (a closed circle: one mode;
+a spherical shell: three), and `redundant`/LU hits a zero pivot there —
+`SUBPC_ERROR`, outer reason −11. So the rotated route asks for
+`geometric_mg_bundle(coarse="svd")`, which is null-space robust and cheap on a
+small coarse level. This is a named variant of the shared bundle, not a
+call-site override, so it is visible in the same place as everything else.
+
+The other native/custom-P asymmetry — that native FMG is unusable for
+single-field solvers because `DMCreateInjection` cannot reliably be built on a
+refined DMPlex (#276) — is deliberately *not* in the bundle. It is a routing
+decision (which route a solver may take), not an option value, and lives with
+the route choice in `_apply_preconditioner_options`.
+
+### Testing it
+
+`tests/test_1021_mg_option_bundle.py` reads the smoother configuration back off
+the **live PETSc objects** after setup — not out of the options database — for
+all three routes and asserts they agree, with the coarse-solve difference
+asserted rather than tolerated. Options-database assertions would not have caught
+the original drift, because the drift was precisely a key that was never written.
+
 ## Critical Stability Note
 
 ```{warning} Solver Stability is Paramount
@@ -337,7 +524,7 @@ pressure Krylov iteration. That number is a sample, not work. `SolverInstrumenta
 
 ```{note} For Contributors
 This well-documented subsystem could benefit from:
-- Preconditioner selection guidance (Krylov sub-solve choice is now covered; multigrid and Schur preconditioner selection is not)
+- Preconditioner selection guidance (Krylov sub-solve choice and the multigrid option bundle are now covered; Schur preconditioner selection is not)
 - Performance tuning documentation  
 - Convergence analysis examples
 - Scaling studies and optimization