Preconditioner defaults: PETSc's GAMG cycle, node aggregation, and a flexible outer Krylov

[lmoresi]

Two preconditioner defaults that deviate from PETSc's without a recorded reason, and one consequence of the first.

The GAMG cycle (#579)

pc_mg_type was additive. PETSc's default for GAMG is multiplicative, and additive is a strictly weaker cycle — the levels do not see each other's corrections.

Provenance, from git log -S: it arrived in December 2022 in example scripts, alongside pc_mg_levels = 2. With a single coarse level the two cycles very nearly coincide, so it cost nothing where it was written. It was promoted into the solver defaults, the hierarchies got deeper, and nobody revisited it. No comment ever gave a reason.

Measured on SolKz, Taylor–Hood P2–P1, 11 727 unknowns, cycles per velocity solve: 243–347 additive against 74–118 multiplicative.

The four open-coded copies

utilities/multigrid_options.py says of itself that "nobody writes a multigrid option value anywhere else", and that is the reason the module exists — hand-maintained copies are what let the smoother iteration count drift in #468. But four call sites in petsc_generic_snes_solvers.pyx open-coded the GAMG bundle key by key, which is how one value came to be wrong in four places at once.

Those four now read the bundle. The module's claim about itself is true again.

The outer Krylov (#576)

The Stokes outer ksp_type was never set, so it inherited PETSc's gmres, while both fieldsplit sub-blocks are fgmres run to a tolerance.

That is the wrong way round. Inside the Schur factorisation the velocity sub-solve computes the search direction, so the operator the outer method applies differs between outer iterations, and GMRES's residual recurrence assumes it does not. Measured on the Spiegelman notch at refinement 3: 983 velocity iterations per step and DIVERGED_LINEAR_SOLVE, against 58 for fgmres and 23 with a loosened inner tolerance. Raising the velocity iteration cap is inert — byte-identical residuals — so the failure is inconsistency, not a shortage of iterations.

Set explicitly in __init__ and in the strategy setter, which also resets the fieldsplit block.

Node aggregation

The velocity operator reached GAMG with block size 1 on a two-component field, so it aggregated scalars rather than nodes.

This is not an unconditional defect. PETSc infers the block size correctly wherever it can:

BCs on velocity	velocity IS	fieldsplit sub-matrix
none	2	2
full vector Dirichlet	2	2
component-wise (free slip)	1	1

It collapses to 1 under a component-wise Dirichlet condition because the field then genuinely carries one degree of freedom at a constrained wall node and two inside — block size 1 is the correct description. Component-wise Dirichlet is free slip, so the one case PETSc cannot infer is the standard geodynamics boundary condition.

Setting mat_block_size asks GAMG to aggregate on a pairing that does not align with nodes at the walls. That cannot change the answer — a preconditioner alters the route to the solution, not the solution — and measured on free-slip SolKz the two configurations agree to 2e-7 relative on a 1e-6 solve. What it changes is the cost:

free-slip SolKz, 11 727 unknowns	cycles/velocity solve	flops
block size 1	74–118	3.48e9
block size 2	34–47	1.60e9

With the same contrast concentrated in a band rather than spread smoothly, a plain stokes.solve() goes from 2.61e11 flops / 32.7 s to 4.44e9 / 1.8 s.

PETSc honours mat_block_size as an option, so this needs no hook and no post-hoc PC surgery — only that it is set before the matrix is built. It is recorded on the solver as _pc_block_size, not just at the __init__ call sites, because _apply_preconditioner_options re-applies the bundle on every build and a bundle built without it lists mat_block_size as a stale key and deletes what __init__ set.

SNES_MultiComponent keeps block size 1: its unknown is not mesh.dim components in general, and claiming the wrong node size would cost rather than save.

One test expectation changed

test_0203_solver_wallclock_guard::test_first_solve_admits_its_count_is_a_lower_bound asserted that a fresh solver's first velocity count strictly undercounts the second. The contract is that the first count is a lower bound and is flagged as one; how far short it falls depends on where PCSetUp lands relative to the first application of the block, which these defaults move. It is now tight on that problem (80 either way, still flagged incomplete), so the test asserts the bound rather than the size of the gap. Verified to pass on development before the change and to fail only because of it.

Still unaudited

pc_gamg_agg_nsmooths=2 (PETSc: 1), pc_gamg_repartition=True (PETSc: False), and the geometric bundle's pc_mg_type="full" (PETSc PCMG: multiplicative) remain deviations without recorded measurements.

Verification

./uw test: 1493 passed, 32 skipped, 2 xfailed. All three defaults confirmed live after ./uw build, including that strategy= no longer reverts the outer.

Underworld development team with AI support from Claude Code

[Copilot]

Copilot was unable to review this pull request because the user who requested the review has reached their quota limit.

[lmoresi]

Adversarial review — solver option-ownership cluster (#584, #597, #548)

Reviewed together: #584 and #597 both edit petsc_generic_snes_solvers.pyx in
the machinery that decides who owns a PETSc option, and #548 is the third open
change to the rotated solve. The regions are textually disjoint — #584 at or
below line 6027, #597 at or above 6783 — so they merge cleanly. The interaction
is semantic.

Cluster

C1. There are two ownership mechanisms, and these PRs touch one each.

	defined on	keyed by	question it answers
_managed_pc_options / _push_managed_option	SolverBaseClass (241, 363)	the global option name	"did we write this, or may the MG bundle back off?"
_owned_option_pushes / _owned_option_user / _resolve_owned_option	SNES_Stokes_SaddlePt (6005, 6774)	the prefixed key	"did we push this, or did the user set it?"

#584 extends the first (_pc_block_size flows into the GAMG bundle and out
through _push_managed_option). #597 fixes a latch defect in the second. Both
are answering the same question — did this value come from us or from the user —
with different state, different key spaces and different lifetimes.

The defect #597 fixes is that the latch outlived the option it was latched from.
The managed mechanism has the same shape: a dict recording what we wrote, read
later to decide whether to defer. It has no equivalent test, and nothing checks
what happens when a user deletes a key it has recorded. We are not asking either
PR to unify them, but whoever does should know they are two, and that only one
of them now has a regression test.

C2. #584 is red and the fix exists. Four failures, all
Local size N not compatible with block size d out of PetscLayoutSetBlockSize
via MatSetFromOptions on the velocity sub-matrix, plus
test_1021_mg_option_bundle asserting on a stale key set. Diagnosis and fix are
on bugfix/fix584-block-size: divisibility is a property of the particular
combination of boundary conditions, not of their kind, and it is mixed across
ranks ([127, 150] at np=2), so the gate has to be collective. That branch is
green at 1495 passed. #584 should not be read as "needs debugging" — it needs
that branch or an equivalent.

C3. #593 also edits rotated_bc.py. Docstring only, describing the new
rank-zero gather in boundary_flux. No functional overlap with #548; recorded so
the file collision is not mistaken for one.

#584

1. The measured win is narrower than the numbers imply. The 74–118 to 34–47
cycle improvement is from free-slip SolKz, whose velocity block divides. The
asymmetric boundary-condition mixes do not divide, and with the gate applied they
run without node aggregation — they still get the multiplicative cycle and the
flexible outer, but not the factor the section is named for. The PR text should
say which configurations get which.

2. SNES_Vector has the same exposure and no test. It sets
mat_block_size = mesh.dim on the main matrix, whose local size loses
constrained DOFs the same way. The gate on bugfix/fix584-block-size covers it;
nothing demonstrates it needed covering.

#597 — responses to our own earlier findings

Three findings were raised on this PR by us. Positions, so they are not left
open:

• Testing through the private resolve/push pair rather than solve() —
  keeping it. The pair is what solve() calls, the test runs in a second
  instead of three solves, and the alternative pins snes.getIterationNumber(),
  which is a weaker assertion about a stronger path. Recorded as a judgement, not
  an oversight.
• The sympy branch of the viscosity_min_rounding guard admits a negative
  constant — worth closing, and cheap: sympy.sympify(value).is_negative is
  True for Float(-1.0) and None for a free expression, so the check can cover
  the constant case without rejecting the δ atoms the property exists for. Not
  done in this PR; it is a one-line follow-up rather than a reason to hold it.
• Item 4 documented rather than rebound — standing. The orphaning is
  unreachable only because yield_continuation refuses anchor= alongside a
  ready-made control, which is a guard in another module. If that is relaxed the
  comment is all that remains. A setter that raises when a control is attached
  would close it by construction, and that is a change to yield_anchor: which side of exact Min the soft-min yield law sits on #475's design rather
  than to this follow-up batch.

#548

Its existing review stands; nothing in this cluster changes it. Note only that
it and #593 both touch rotated_bc.py, per C3.

Underworld development team with AI support from Claude Code

[lmoresi]

The block-size failures, fixed, and the branch brought current

Pushed 8ed2558. This branch was 57 commits behind development; it now merges
cleanly and carries the fix for both failure modes.

The four PETSc error 75 failures. mat_block_size is declared as a PETSc
option, so when PCFieldSplit extracts the velocity sub-matrix,
MatSetFromOptions hands it to PetscLayoutSetBlockSize — a hard error rather
than a hint it may decline:

Arguments are incompatible
Local size 67 not compatible with block size 2      (2-D)
Local size 1384 not compatible with block size 3    (3-D)

The commentary in _gamg_settings anticipated the misalignment under
component-wise Dirichlet but not the refusal. The reason the SolKz measurements
never met it is narrower than "component-wise Dirichlet": divisibility is a
property of the particular COMBINATION of boundary conditions, not of their
kind. On a 3x3 P2 velocity field —

velocity BCs	local size
none	98	even
full vector Dirichlet, every wall	50	even
component-wise free slip	70	even
(0, 0) Bottom, (0, None) elsewhere	67	odd, refused

Free slip is safe by an accident of counting, and it is the case the
measurements used. test_1010_stokesCart's asymmetric mix is not.

_withdraw_block_size_if_not_node_blocked removes the option when the block is
not node-blocked, immediately before setFromOptions — the first point the
field decomposition exists — and on every build, since
_apply_preconditioner_options re-pushes the bundle each time.

The decision is collective, because divisibility is mixed across ranks. On a
6x6 box the per-rank velocity sizes are [127, 150] at np=2 and
[72, 94, 111] at np=3, so a rank-local test would leave ranks disagreeing
about the contents of the options DB.

The fifth failure was a test defect.
test_1021_mg_option_bundle::test_bundles_clear_each_others_keys built owned
from a bundle list that never included gamg_bundle(block_size=2) — the one
bundle that sets mat_block_size — so every bundle's stale list carried a key
owned lacked. The sibling is now in the list, which also brings that bundle
under the test.

test_1022_velocity_block_size_gate.py pins both directions. The second
case is the one that matters: a gate that withdrew the option unconditionally
would fix the crash and silently undo everything #579 was for, and the crash
test alone cannot tell them apart.

Two things for the PR text. The measured 74-118 to 34-47 cycle win is from
free-slip SolKz, which divides; the asymmetric mixes now run without node
aggregation, so they get the multiplicative cycle and the flexible outer but not
that factor. And SNES_Vector has the same exposure — the gate covers it, but
nothing demonstrates it needed covering.

Verified: full ./uw test 1556 passed, 32 skipped, 2 xfailed, against
development at fe6b2e1.

Underworld development team with AI support from Claude Code