Skip to content

custom-P FMG smoother has drifted from the native FMG bundle (richardson vs gmres, and an accidental iteration count) #468

Description

@lmoresi

Two code paths write the geometric-multigrid option bundle for the Stokes velocity block, and they have drifted apart. The custom-P path — which is the mandatory route for rotated boundary conditions and for adapt() children, because native DMPlex MG cannot serve either — is running the configuration the native path measured and rejected.

Measured, read back off the live PC after setup

setting A native (DMPlex refinement) B custom-P (standard path) C custom-P (rotated path)
pc / levels / mg_type mg / 2 / FULL mg / 2 / FULL mg / 2 / FULL
smoother KSP gmres richardson richardson
smoother PC sor sor sor
smoother max_it 4 3 2
smoother norm type none none none
coarse KSP / PC preonly / redundant preonly / redundant preonly / svd

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

The coarse svd on the rotated path is deliberate and correct — the Galerkin-coarsened rotated velocity block inherits the rigid-rotation null space and redundant/LU hits a zero pivot (the #306 fix). Not drift. The other two rows are.

1. The smoother kind

petsc_generic_snes_solvers._configure_preconditioner sets gmres+sor with a recorded, measured justification:

Chebyshev needs eigenvalue estimates of the smoothed operator, which are fragile on the indefinite / variable-viscosity velocity block and diverge. Richardson is stationary and degrades on the NON-SYMMETRIC operator produced by the consistent-Newton tangent. Measured on the Spiegelman notch (Drucker-Prager, eta contrast 1e26) over a nested 4-level hierarchy: contraction per V-cycle rho = 0.75 (richardson) vs 0.56 (gmres) at the SAME four smoother iterations — and the gmres margin GROWS with depth (5% at 3 levels, 25% at 4).

custom_mg._configure_pcmg still writes richardson+sor. It never received that update.

2. The iteration count is accidental

_configure_pcmg never sets mg_levels_ksp_max_it at all, so the smoother inherits whatever is left in the options database under that prefix:

  • standard path → 3, left over from the GAMG bundle (opts[f"{prefix}mg_levels_ksp_max_it"] = 3), which _configure_pcmg does not clear;
  • rotated path → 2, PETSc's own PCMG default, because that prefix is fresh and never saw the GAMG bundle.

So the same function produces different smoothing depending on what ran before it. That is the drift mechanism, not just the symptom.

Cost

Same rotated operator and right-hand side, transversely isotropic annulus with a weak plane reaching the constrained boundary (eta_1/eta_0 = 1e-3), timed sequentially in-process, three repeats:

smoother velocity its median s
richardson/2 (rotated path today) 11 0.896
richardson/3 (standard custom-P path) 9 0.720
richardson/4 8 0.784
gmres/4 (native bundle) 5 0.525

1.7x, and the native path's Spiegelman-notch finding transfers directly to the anisotropic rotated case. Note this is a 2-level hierarchy — the native comment says the gmres margin grows with depth, so a deep adapt-child hierarchy should do better still.

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

Suggested fix

Give _configure_pcmg the same smoother bundle as the native path — mg_levels_ksp_type=gmres, mg_levels_pc_type=sor, mg_levels_ksp_max_it=4, mg_levels_ksp_norm_type=none, mg_levels_ksp_converged_maxits — and set the iteration count explicitly rather than inheriting it. Keep the rotated path's mg_coarse_pc_type=svd override.

Better still, and the reason to do it now rather than patch it: have one place own the bundle and both routes read it, so this cannot drift again. The two sites are _configure_preconditioner (petsc_generic_snes_solvers.pyx) and _configure_pcmg (utilities/custom_mg.py).

Validation this needs before merging (larger than the change): test_1017_custom_mg_stokes, the adaptivity suites that rely on mesh-owned custom-P hierarchies, and the rotated suites test_1018 / test_1064 / test_1066 — several carry iteration-count guards and goldens that this will move.

Found while auditing the rotated preconditioner (#465); related to #467 (mesh-owned custom-P hierarchy is not picked up on the rotated path at all).

Underworld development team with AI support from Claude Code

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions