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

[lmoresi]

Closes #468. Closes #467. Stacked on #465.

The problem

Three routes reach a multigrid Stokes velocity block, and they are the same preconditioner reached three ways, not alternatives:

route	when	prolongation
native	refinement >= 1, ordinary BCs	PETSc DMCreateInterpolation between refined DMPlex levels
custom-P, standard	set_custom_fmg, or an adapt() child's mesh-owned coarse tail	barycentric / RBF, Galerkin coarse operators
custom-P, rotated	rotated free-slip, via rotated_bc	as above, fine prolongation rotated P̂ = Q_v·P

custom-P is mandatory wherever native cannot go — rotated BCs (the DM-coupled hierarchy cannot express a per-node rotation) and adapt() children (no DMPlex refinement relation). So the drifted bundle was on the route that matters most.

#468 — the bundle was written in two places. The native path had been moved to a measured gmres+sor smoother; custom_mg._configure_pcmg never received that update, and never set mg_levels_ksp_max_it at all, so the smoother inherited whatever last wrote that prefix: 3 on the standard path (left over from the GAMG bundle) and 2 on the rotated path (PETSc's own PCMG default). The same function smoothed differently depending on what had run before it.

#467 — mesh.adapt() leaves a coarse tail on its refinement child so every solver on an adapted mesh gets geometric MG 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 reported velocity_pc == "GAMG" — indistinguishable from having no hierarchy at all. That is the adapt-on-top-faults workflow's own configuration.

The change

src/underworld3/utilities/multigrid_options.py owns both bundles (geometric MG and the GAMG fallback); all three writers read from it. Two design points worth review:

• The stale-key list is derived, not written. A bundle must clear every key a sibling bundle sets that it does not, or setFromOptions re-reads a leftover. Deriving it from the union reproduces all three previous hand-written lists exactly — which is the check that it is the right abstraction, and it is self-maintaining when a key is added.
• The coarse solve is a named variant (geometric_mg_bundle(coarse="svd")), not a call-site override, so the one legitimate per-route difference is visible in the same place as everything else. The other native/custom-P asymmetry — native FMG being unusable for single-field solvers because DMCreateInjection is unreliable on a refined DMPlex (FMG (preconditioner="fmg"/"auto") fails with PETSc err62 (injection) for scalar solvers on a refined Plex hierarchy #276) — deliberately stays out of the bundle: it is a routing decision, not an option value, and belongs with the route choice.

custom_mg.build_transfers is now the shared "which hierarchy does this solver get?" rule for both the standard injection hook and rotated_bc._build_rotated_custom_Pl, with the same opportunistic barycentric→RBF fallback and degrade-to-GAMG on failure.

Also: the rotated FMG branch now drops its option keys after setup, as the KSP's own keys already did. The per-solve prefix is unique, so those keys previously accumulated in the global options database once per timestep.

Measured

Same operator, RHS and coarse solve. Two levels — where the native measurement says the gmres margin is smallest.

route	pre-fix	unified	isolated linear solve
custom-P standard (isotropic annulus)	5 vel its (richardson/3)	4 vel its (gmres/4)	—
custom-P rotated (TI annulus, η₁/η₀ = 1e-3)	11 vel its (richardson/2)	5 vel its (gmres/4)	0.683 s → 0.392 s (1.74×)

Reproduce: ~/+Simulations/rotated_pmat_audit/bundle_unification_gain.py (iterations, cold solve each repeat) and smoother_drift_cost.py (the linear solve timed in isolation, setup warm). fmg_config_diff.py now shows the three routes agreeing on every row except the deliberate coarse_pc = svd.

Tests

tests/test_1021_mg_option_bundle.py reads the smoother configuration back off the live PETSc objects for all three routes and asserts they agree, with the coarse difference asserted rather than tolerated. An options-database assertion would not have caught the original drift — the drift was precisely a key nobody wrote. Verified to fail when the pre-fix bundle is reimposed.

Also covered: the mesh-owned pickup and its negative control; the rotated FMG PC surviving eight Newton increments (the failure mode the option cleanup introduces); and the same mesh-owned pickup at np>1 in tests/parallel/test_1064, where the transfers are built cross-partition.

Validation

Serial: test_1014, test_1015, test_1016, test_1017, test_1018, test_1020, test_0753, test_0835, test_0836, test_0840 all pass. Full level_1 and tier_a sweep: 526 passed, 2 failed — both pre-existing and unrelated (test_0641 still asserts a non-zero rotated datum raises NotImplementedError, which #458 implemented; filed as #470).

Parallel np2 and np4: test_1017_custom_mg_parallel_mpi 6 passed, test_1064_rotated_freeslip_parallel 9 passed, test_1066_rotated_datum_parallel 2 passed.

No goldens moved. The brief anticipated re-recording; none was needed. GOLDEN_ANNULUS_FMG is a converged-solution quantity at 1e-6 rtol, so a better preconditioner reaches the same answer, and the iteration guards (its <= 25, max(ksp_its) <= 3) are upper bounds a better smoother only helps.

Underworld development team with AI support from Claude Code

[lmoresi]

Adversarial review

We reviewed our own diff before marking this ready. Findings, most serious first.

1. The smoother we unified on costs wall clock on the standard route, at every depth we measured

This is the weakest point of the change and it deserves the top of the list rather than a footnote.

#468 asked for the native smoother bundle on the custom-P route. We measured what that costs. Outer Stokes KSP timed in isolation (setup warm), nested annulus hierarchy, linear solve, η contrast 1e6, varying only the velocity-block smoother on the live PC:

depth	richardson/3	gmres/4	its	time
2 levels	16 vel its, 0.752 s	11 vel its, 0.875 s	×1.45	×0.86
3 levels	8 vel its, 2.569 s	5 vel its, 3.346 s	×1.60	×0.77
4 levels	6 vel its, 10.870 s	5 vel its, 19.794 s	×1.20	×0.55

gmres/4 wins on iterations everywhere and loses on wall clock everywhere, by a margin that grows with depth — the opposite of the depth behaviour the native comment records. Decomposing at 2 levels: richardson/3 → richardson/4 is 2.382 → 2.456 s (the fourth sweep, ~3%); richardson/4 → gmres/4 is 2.456 → 3.437 s (gmres itself, ~40% per cycle). Custom-P uses Galerkin RAP coarse operators built from barycentric transfers, denser than the native nested ones, so the smoother is a larger share of the cycle here than on the route the original measurement came from.

We do not think this refutes that measurement. It was taken on the Spiegelman notch with the consistent-Newton tangent — a non-symmetric operator, where richardson is stationary and degrades. The sweep above is linear and symmetric: richardson's best case. The complementary evidence is the rotated TI case in the PR body, which is anisotropic and where gmres/4 wins 11 → 5 iterations and 0.683 → 0.392 s.

The gap we could not close: an attempt to measure the non-symmetric case on the standard route (smoother_depth_scaling.py -uw_nonlinear 1) did not converge (DIVERGED_LINE_SEARCH), so there are no numbers for it. The script now refuses to time a failed solve rather than report whatever operator the line search abandoned. Someone should close this with a standard-path non-symmetric case that converges.

The call we took: unify on gmres/4 anyway. Richardson's failure on a non-symmetric operator is a robustness failure; gmres's cost on a symmetric one is a speed cost, and defaults err toward robust generality. The consistent-Newton tangent is the recommended setting for nonlinear rheologies, so the symmetric regime is the linear case only. But this is a defaults decision with a measured price on the linear path, and it is @lmoresi's to confirm, not ours.

2. We widened the scope past what #468 asked for

#468 is about the geometric-MG bundle. We also brought the GAMG fallback into the owner, because rotated_bc carried a hand-copied duplicate of the native GAMG bundle (commented "tuned to native parity") — the identical two-writer structure, one drift away from the same failure. Ten lines, and it removes a third copy.

We did not touch the pc_mg_type=additive / mg_levels_ksp_max_it=3 blocks in the per-solver __init__ defaults (petsc_generic_snes_solvers.pyx, five sites). Those are per-solver defaults rather than the managed bundle, and folding them in is a behavioural change well beyond this PR. They remain duplicated.

3. The rotated path now deletes its whole option bundle after setup — a new silent failure mode

Previously the FMG branch deleted one key (mg_coarse_pc_type) and leaked the rest. The prefix is unique per rotated solve, so those keys accumulated in the global options database once per timestep — the module docstring already claims the cleanup keeps it bounded, and the FMG branch did not honour that.

Deleting the rest is only safe if nothing re-reads them: a later setFromOptions on the velocity sub-PC would find pc_type gone and silently abandon the multigrid. We added test_rotated_fmg_survives_repeated_newton_increments in the same commit — eight Newton increments on a power-law rotated solve, all keeping velocity_pc_type == "mg" and vel_its_last flat at 5. It is a check, not a proof: a future PETSc that re-reads options inside PCSetUp would break it, and the test is the thing that would say so.

4. _build_rotated_custom_Pl mutates solver._custom_mg after an opportunistic pickup

Registering a mesh-owned hierarchy on the solver promotes it from opportunistic to solver-set, so a later build failure would raise instead of degrading to GAMG. auto_inject_custom_mg has done exactly this since #424, so this is parity rather than a new semantic — but it is a semantic worth being deliberate about, and we were.

5. The parity test's coverage is narrower than the claim

test_1021 compares the three routes on a 2-level hierarchy in serial. It would catch the drift this PR fixes (verified: it fails when the pre-fix bundle is reimposed) but it does not exercise depth, and the native route it compares against is the only one we can build natively on this geometry. The mesh-owned pickup is covered at np>1 in test_1064; the bundle parity is not.

6. What we checked and found clean

• The derived stale-key lists reproduce all three previous hand-written lists exactly — that equality is the evidence the derivation is the right abstraction, not a re-implementation with a different footprint.
• Value types crossing into the shared bundle: pc_gamg_repartition moved from the string "true" to True, and mg_levels_ksp_converged_maxits from "true" to a bare flag. Both reach PETSc identically (setValue maps True → "true", None → a set flag, and a bare bool flag is TRUE). The GAMG-route tests in test_1064 pass at np2 and np4 against unchanged goldens, which is the behavioural check.
• _configure_pcmg's new coarse= keyword defaults to the previous behaviour; the two _install_transfers call sites are unaffected.
• No goldens moved. We expected to re-record and did not have to: GOLDEN_ANNULUS_FMG is a converged-solution quantity at 1e-6 rtol, so a better preconditioner reaches the same answer, and the iteration guards are upper bounds.
• Two failures in the level_1 and tier_a sweep are pre-existing and unrelated (Stale rotated free-slip datum guard tests in test_0641 (nonzero datum is now implemented, #458) #470).

Underworld development team with AI support from Claude Code

[lmoresi]

Follow-up: the smoother's recorded mechanism is wrong (2550d9a)

@lmoresi queried the "non-symmetric" framing in finding 1 above. He was right, and it does not survive measurement. Correcting it is now part of this PR — comments and docs only, no behavioural change.

The claim. Both the bundle comment we moved into multigrid_options.py and the Layer 3 section of design/nonlinear-solver-homotopy-warmstart.md state that richardson degrades because the consistent-Newton tangent makes the velocity block non-symmetric. We inherited it without checking.

On paper it is wrong. For an isotropic η(ε̇_II) the Newton term is 2 (η'/ε̇_II) (ε̇(u):ε̇(δu)) (ε̇(u):ε̇(v)) — a rank-one outer product a⊗a with a = ε̇(u), symmetric under swapping δu and v. η depends on ∇v only through its symmetric part, so both factors project onto the same tensor.

Measured, ||A - A^T||_F / ||A||_F on the assembled velocity block:

| rheology | asymmetry | 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 — a null symmetry reading is vacuous if the consistent tangent contributed nothing. It contributed 5-40% of ||A|| every time, and the block stayed symmetric anyway.

Non-symmetry appears only under transverse isotropy, and there it appears under the Picard tangent too, where int ε̇(δu):C:ε̇(v) is symmetric by construction. Linear TI with anisotropy fully active is symmetric at 6.3e-17. That is a defect signature, not a property of the tangent — reported on #457, whose scope as written ("the anisotropy part of the consistent tangent") now looks too narrow.

A trap worth recording: petsc4py's Mat.transpose() transposes IN PLACE and returns self. Called on the matrix rather than a copy it measures A - A and reports an exact zero for every operator. Our first run "confirmed" everything was symmetric that way, and the tell was that it was exactly zero in cases that had no business agreeing.

What this does and does not change.

The Spiegelman-notch contraction measurement (ρ 0.75 richardson vs 0.56 gmres, margin growing with depth) is untouched — it is a measurement, not a derivation. Only the explanation attached to it was wrong. The mechanism that fits both it and the cost in finding 1 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. That predicts the notch result and predicts that gmres loses on the well-conditioned problems we measured.

It also dissolves the gap we flagged in finding 1. We were looking for a "standard path × non-symmetric operator" case to justify gmres there. That cell does not exist for isotropic rheologies — the standard-path block is symmetric whether the tangent is Picard or Newton. So the discriminator is not symmetry but conditioning, and the honest position is:

• gmres/4 is the right default because it is the one that does not stall on high-contrast operators — the regime this code exists for;
• it costs up to ×1.8 in the linear solve on well-conditioned problems, at every depth we measured;
• both facts are now recorded in multigrid_options.geometric_mg_bundle and subsystems/solvers.md so removing the guardrail is a deliberate act.

The design doc's Layer 3 claim is marked with a TODO(DESIGN) rather than rewritten — it is #441's document, and the correction is bigger than a drive-by edit.

Underworld development team with AI support from Claude Code

[lmoresi]

Scope note: this does not un-gate geometric MG for single-field solvers (#478)

A parallel session flagged that geometric FMG is currently gated off. Checking what this PR does and does not close: the gate is real, it is the #276 single-field lockout, and #471 does not remove it — deliberately, as the one design question the brief asked us to settle.

Measured off the live PC after setup:

solver	mesh	preconditioner=	live PC
Stokes	refinement=1	auto	mg(2)
Stokes	refinement=2	auto	mg(3)
Poisson	refinement=1	auto	gamg
Poisson	refinement=2	fmg (explicit)	gamg + warning

Stokes is fine. Every single-field solver gets algebraic multigrid at any depth, even when geometric is requested explicitly, via if want_fmg and prefix == "": want_fmg = False.

We kept that out of the shared bundle owner on the grounds that it is a routing decision — which route a solver may take — not an option value, and so belongs with the route choice rather than with the settings. That reasoning still holds. But it does mean the headline "all three routes now agree" is scoped to the routes a solver is allowed to reach, and for scalar and vector solvers there is only one.

What this PR does change for them: the sanctioned workaround (set_custom_fmg, which the lockout warning points at) previously installed the drifted bundle — richardson at an iteration count nobody set. It now installs the same bundle as everything else. Same for the mesh-owned auto-pickup on adapt children, which already worked for single-field solvers. So single-field geometric MG got better where it was reachable; it did not become reachable.

Filed as #478 with a proposed fix, because #471 makes it small: mesh.dm_hierarchy is already a coarse tail, barycentric transfers with Galerkin RAP need no injection (which is the whole of err62), and custom_mg.build_transfers is now the single hierarchy-resolution rule with a safe degrade. It wants its own PR — it is a default change across every scalar and vector solver.

Repro: ~/+Simulations/rotated_pmat_audit/fmg_gate_audit.py.

Underworld development team with AI support from Claude Code