Rotated free-slip: the custom-FMG velocity sub-KSP must converge, not just apply

[lmoresi]

Stacked on #458 (feature/rotated-datum-snes). Retarget to development when that merges — do not delete the base branch first, it would close this PR.

Why

Louis asked whether the rotated free-slip BC, which rotates the matrices after they are assembled, does that successfully for the preconditioner matrix as well — the concern being that an error there would be invisible under isotropy (a viscous operator is nearly rotation-invariant) but not under anisotropy, and the free-surface + transverse-isotropy fault model had flagged something.

The audit answer: the preconditioner is fully rotated

An algebraic identity probe reaches into the live PC after PCSetUp and compares every matrix it holds against two independently-formed references: d_rot (the rotated block, built by the opposite route to the solver — block the native Jacobian, then rotate) and d_unrot (the same block un-rotated, which serves as the negative control).

matrix	d_rot	d_unrot
velocity sub-KSP operator (Amat)	0.0	7.1e-2
velocity sub-PC operator (Pmat)	0.0	7.1e-2
MG level 1 (fine)	0.0	7.1e-2
MG level 0 (Galerkin coarse)	0.0	2.3e-1

Identical for the isotropic control. So:

• The rotated KSP has no Pmat distinct from Amat — ksp.setOperators(Ahat) is single-argument, so the velocity preconditioner is the rotated operator.
• The Galerkin coarse operators inherit the rotation correctly even though only the finest prolongation is rotated. This was the one place where "rotated operator, un-rotated preconditioner" could still have been true in effect.
• The Schur preconditioner is correctly un-rotated: Q is the identity on pressure, and the p-p block of the rotated operator is identically zero, so it must come from the native Pmat. UW3's Pmat differs from its Amat only in that p-p block (the velocity and coupling blocks are registered with identical pointwise functions, petsc_generic_snes_solvers.pyx:7803-7809), so discarding the rest costs nothing.

No un-rotated matrix reaches the preconditioner. The original hypothesis is not the defect.

The defect that is there

The custom-FMG velocity block was driven by preonly — "a full-MG cycle per Schur application, by design". But PCFieldSplit forms S = A11 - A10 A00^-1 A01 and applies A00^-1 through that same velocity KSP, so preonly replaces A00^-1 with a single multigrid cycle and hands the pressure Krylov a different system S~ != S, preconditioned by a 1/mu mass matrix built for S.

Pressure sub-KSP residual history (TI annulus, weak plane reaching the constrained boundary, eta_1/eta_0 = 1e-3):

preonly:  it 0 |r|=1.354e-02 -> it 16 |r|=3.112e-07 -> it 192 |r|=3.111e-07
          reduction 4.4e4, monotone=False -- 184 iterations for nothing
fgmres:   it 0 |r|=1.076e-01 -> it 17 |r|=9.757e-10
          reduction 1.1e8, monotone=True

It stagnates at a floor, on every outer iteration. Not slow convergence, and not a moving operator: applying the Schur operator twice to the same vector is bitwise identical under both settings. S~ is a fixed linear operator, just the wrong one. Why the floor sits at 3.1e-7 is not isolated here (a range/null-space inconsistency between S~ and the constant-pressure null space attached to S is the obvious suspect).

Cost, timed sequentially in-process, three repeats:

velocity sub-KSP	outer	pres	best (s)	median (s)
preonly (shipped)	9	200 (cap)	2.767	2.968
fgmres 0.1 x tol	1	17	0.672	0.705

4.1x. The isotropic control moves the same way (5 -> 1 outer), so this is the Schur application rather than the anisotropy — the anisotropy only makes it hurt more. On the test_1018 FMG configuration (SolCx box, 2-level barycentric hierarchy): 6 outer -> 1.

Is this a general problem with our FMG defaults?

No — the rotated custom-FMG branch was the sole exception. Verified at runtime on the native path with the identical hierarchy: velocity sub-KSP is fgmres with PC mg, outer 1 / vel 6 / pres 68. The native path sets that KSP type explicitly and custom_mg._configure_pcmg only replaces the PC; the rotated GAMG fallback was already FGMRES.

The transferable rule, now in docs/developer/subsystems/rotated-freeslip.md: a velocity sub-KSP that is preonly under pc_fieldsplit_schur_fact_type=full gives the pressure Krylov the wrong system to solve. Multigrid is fine as the preconditioner; it cannot be the whole solve when a Schur complement is applied through it.

Changes

• Wrap the custom-FMG velocity block in FGMRES. max_it 200 matches the GAMG fallback and the native path; rtol 0.1 x tol matches the GAMG fallback (the native path asks 0.033 x tol).
• Warn when the velocity sub-KSP exhausts its cap — it returns KSP_DIVERGED_ITS, which KSPCheckSolve deliberately does not escalate, so FGMRES can degrade silently where preonly could not fail.
• Report velocity_pc, schur_pre, velocity_pc_type (PETSc's own view of the sub-PC) and vel_its_last / pres_its_last in _rotated_freeslip_info. The outer count alone hides a degraded inner solve — a full Schur factorisation with a good pressure mass still reports ~1 outer iteration while the inner solve grinds underneath. That is how this went unnoticed. The _last names are deliberate: these are last-application samples, not the summed work axis solver_health uses. Wiring the rotated path into SolverInstrumentation.sub_reports() is the proper fix and is not done here.
• Guard the test_1018 FMG configuration against a return to preonly, and assert the velocity sub-PC really is multigrid rather than restating our own bookkeeping.
• Fix a pre-existing test_1064 failure: int(info["ksp_its"]) on a value the unified rotated loop of Rotated free-slip: prescribed-normal datum through the SNES path; one unified Newton loop (#438, #403) #458 made a list, so the 3D spherical partition-independence test errored at np>1 before reaching any of its assertions.
• Document the interaction in the rotated-freeslip subsystem page.

Also measured, no change made

The constraint-row diagonal (zeroRowsColumns(normal_rows, diag=s)) is a free conditioning parameter — any positive s is exact. Before the fix it was a real lever on the GAMG route (37 velocity iterations at s <= 1e-2 x mean, 44 at the shipped mean, 142 at 1e4 x mean), and the isotropic control had the same shape, so it was general rather than an anisotropy defect. After the fix, on the FMG route it is completely inert: 11 velocity / 17 pressure iterations at every factor from 1e-4 to 1e4. The sensitivity was a symptom of the degraded configuration, so nothing was changed here.

Tests

test_1018 18 passed; test_1064 + test_1066 10 passed at np2 and np4; the TI fault smoke passes (24 nonlinear iterations, datum relL2 4.4e-05).

Its 24 iterations are #457 (the TI consistent tangent is wrong when anisotropy is active), a native-path defect that reproduces with no rotation involved. It is not a preconditioner problem and should not be confounded with one.

Full study, scripts and numbers: ~/+Simulations/rotated_pmat_audit/README.md. Adversarial review in the comment below.

Filed separately while auditing: #463 (TI-VEP .viscosity discards the yield-limited eta_1), #464 (block-constrained Stokes + rotated free-slip is unguarded).

Underworld development team with AI support from Claude Code

[lmoresi]

Adversarial review

Findings against the first push (d9f76516), and what we did about each. Amended in 4a557b81.

1. The stated mechanism was inverted — accepted, corrected

We wrote that preonly makes the Schur operator "not fixed from one application to the next", so the pressure Krylov "cannot converge". That is wrong. The FMG cycle is richardson/sor smoothing over a redundant-LU coarse solve — a fixed linear operator — so S̃ = A₁₁ − A₁₀ M A₀₁ is itself a fixed linear operator.

Measured directly: applying the Schur operator twice to the same vector gives ||Sx − Sx||/||Sx|| = 0.000e+00 under both preonly and FGMRES.

What the residual history actually shows (pressure sub-KSP monitor, TI annulus):

preonly:  it 0 |r|=1.354e-02 -> it 16 |r|=3.112e-07 -> it 192 |r|=3.111e-07
          r0->rN reduction 4.4e4, monotone=False, 184 iterations for nothing
fgmres:   it 0 |r|=1.076e-01 -> it 17 |r|=9.757e-10
          reduction 1.1e8, monotone=True

So it is stagnation at a floor, not slow convergence and not non-convergence against a moving operator. preonly hands the pressure Krylov a different system S̃ ≠ S, preconditioned by a 1/μ mass built for S, and it stalls above its tolerance. We have not isolated why the floor sits at 3.1e-7 (a range/null-space inconsistency between S̃ and the constant-pressure null space attached to S is the obvious suspect). The commit message, code comment and docs now say exactly this and no more.

2. FGMRES can now fail silently where preonly could not — accepted, fixed

KSPCheckSolve deliberately exempts KSP_DIVERGED_ITS, so a velocity FGMRES that exhausts its cap does not escalate. The rotated path now calls _warn_if_ksp_diverged on the velocity sub-KSP.

3. The parity claim was numerically false — accepted, fixed

max_it was 100 against 200 in both the GAMG fallback and the native path; rtol is 0.1×tol against the native path's 0.033×tol (petsc_generic_snes_solvers.pyx:6311). max_it is now 200. rtol stays at 0.1×tol (matching the GAMG fallback) because the FMG cycle reaches it in ~11 iterations; the comment no longer claims sameness with both.

4. velocity_pc="direct-LU" mislabels a converged warm start — accepted, fixed

The defaults were keyed on ctx is None, but the iterative path also reaches the return with ctx unset when a warm start is already converged at increment 0. Now keyed on use_lu.

5. vel_its / pres_its are last-application samples, not work — accepted, renamed

KSPGetIterationNumber reports the most recent solve, and the velocity KSP is applied once per Schur MatMult. Reporting these as its contradicts the convention in solver_health.py where its is summed. Renamed vel_its_last / pres_its_last and documented as samples. Wiring the rotated path into SolverInstrumentation.sub_reports() is the proper fix and is not done here.

6. The new test assertions were tautological — accepted, fixed

velocity_pc == "custom-FMG" and schur_pre == "1/mu-mass" restate our own bookkeeping (custom_Pl is not None, Mp is not None) and cannot detect a failed MG install. We now also report and assert velocity_pc_type, which is vel_ksp.getPC().getType() — PETSc's own view.

7. Iteration counts do not establish cost — accepted, measured

9 outer x 200 pressure iterations at one MG cycle each, against 1 outer x 17 at ~11 cycles each, could in principle go either way. Timed sequentially in-process, three repeats, best/median:

velocity sub-KSP	outer	best (s)	median (s)
preonly	9	2.767	2.968
fgmres 0.1×tol	1	0.672	0.705

4.1x. The iteration-count story holds.

8. max(info["ksp_its"]) on an empty list — accepted as latent, not fixed

A zero-increment exit (warm start already converged) leaves ksp_its == [] and max raises. Not reachable from these cold-start tests, and the same pattern pre-dates this PR at test_1018:101. Left alone rather than widened here.

9. max(...) vs sum(...) in test_1064 — accepted, left as is

_capture_rotated_report sums; our one-line fix uses max, matching the pre-existing test_1018:101. It weakens its <= 25 from a total to a per-increment bound, with no practical difference for this ksponly single-increment case. Consistency across the two test files was the deciding factor; worth revisiting if ksp_its semantics are ever settled properly.

10. GOLDEN_ANNULUS_FMG sits on the changed path and was not re-recorded — checked

Correct that it is the only golden on this path. It passes unchanged at np2 and np4, which is the check; we did not re-record it.

11. Options-database growth — real, pre-existing, not ours to fix here

_configure_pcmg writes ~7 keys under the unique-per-solve rotfsN_fieldsplit_vel_ prefix and only mg_coarse_pc_type is removed, so the database grows under time-stepping. The "stays clean (and bounded under time-stepping)" comment in the block we touched is therefore false as written. Pre-existing and outside this change; flagging rather than fixing.

Not accepted, for the record: nothing. Every finding above was either a real defect or a real overclaim.

Underworld development team with AI support from Claude Code

[lmoresi]

Follow-up: sub-solve tolerances brought to native parity, and the Krylov axes written down

Louis's steer: defaults should err on the side of robust generality, with callers removing guardrails at their own risk. Applying that here found two more places where the rotated path was the looser of two paths solving the same problem, with no recorded reason:

sub-solve	rotated (was)	native
velocity rtol	0.1 × tol	0.033 × tol
pressure rtol	1.0 × tol	0.1 × tol

Measured cost of matching native, timed sequentially in-process, three repeats:

configuration	outer	best (s) TI/FMG	TI/GAMG	iso/FMG	iso/GAMG
rotated (was) vel .1 pres 1	1	0.577	0.258	0.999	0.732
native parity vel .033 pres .1	1	0.677	0.303	1.300	0.871

~17–30%, with identical outer iteration counts in every case. Adopted, because a too-loose inner solve does not degrade gracefully — it walks back toward the preonly failure this branch already fixes (velocity rtol 1e-1 → 6 outer, 1e-2 → 4, 0.033×tol → 1), and it fails by silent stagnation rather than loudly.

Unexpected bonus: the TI fault smoke dropped from 24 → 15 nonlinear iterations. Tighter inner solves feed a cleaner Newton step. Worth noting against #457 — some of that 24 was loose inner solves rather than the tangent alone, though the remaining Picard-pace behaviour is still consistent with the TI tangent defect.

Cheaper remains available to anyone who measures their own configuration; it is no longer the default.

Documentation

docs/developer/subsystems/solvers.md listed "preconditioner selection missing" as a known gap. Added a section there laying out the three axes that keep getting argued as one:

1. Flexible vs non-flexible (fgmres vs gmres/cg/fcg) — the question is whether the preconditioner is stationary, not whether the operator is symmetric. The velocity block is SPD so cg is admissible on symmetry grounds and still fails, because GAMG with mg_levels_ksp_converged_maxits is a non-linear application. Settled by Kramer free-slip Nitsche fieldsplit solve needs flexible subsolves on Gadi #147 at np=144 on Gadi.
2. preonly vs an iterative wrapper — positional, not a matter of taste. A preconditioner application may be preonly; an operator inverse underneath a Schur complement may not, because it defines the system the pressure Krylov solves rather than merely conditioning it.
3. Inner tolerance — the genuine judgment call, decided by the guardrail policy above.

Plus how to detect a degraded sub-solve, since all three failures are quiet: the outer count hides them, an exhausted cap returns KSP_DIVERGED_ITS which KSPCheckSolve deliberately does not escalate, and KSPGetIterationNumber on a sub-KSP is a last-application sample rather than work.

Tests: test_1018 18 passed; test_1064 + test_1066 10 passed at np2 and np4.

Underworld development team with AI support from Claude Code