diff --git a/docs/advanced/figures/bench_fmg_vs_gamg_velocity_ksp.csv b/docs/advanced/figures/bench_fmg_vs_gamg_velocity_ksp.csv new file mode 100644 index 000000000..3b28a7d47 --- /dev/null +++ b/docs/advanced/figures/bench_fmg_vs_gamg_velocity_ksp.csv @@ -0,0 +1,51 @@ +adapt_step,fmg_vel_ksp_its,fmg_snes_reason,fmg_t_stokes_s,gamg_vel_ksp_its,gamg_snes_reason,gamg_t_stokes_s +1,5,3,2.97,113,3,4.02 +2,4,3,2.83,64,3,18.54 +3,3,3,5.85,111,3,6.18 +4,3,3,4.79,112,3,6.00 +5,4,3,2.93,71,3,15.83 +6,4,3,3.11,87,3,17.88 +7,4,3,2.98,90,3,17.81 +8,5,3,4.57,107,3,6.41 +9,4,3,10.12,89,3,15.76 +10,4,3,9.46,90,3,16.02 +11,4,3,9.42,112,3,12.17 +12,5,3,7.87,112,3,16.04 +13,5,3,5.01,79,3,17.48 +14,5,3,8.23,111,3,16.68 +15,4,3,8.09,113,3,18.80 +16,5,3,8.29,116,3,15.58 +17,5,3,8.23,116,3,18.02 +18,5,3,9.10,113,3,21.25 +19,5,3,10.11,115,3,21.99 +20,6,3,9.84,69,3,24.06 +21,6,3,10.44,123,3,22.41 +22,6,3,11.00,122,3,24.98 +23,5,3,12.64,117,3,25.87 +24,5,3,13.22,112,3,21.44 +25,5,3,16.53,87,3,37.25 +26,5,3,15.86,118,3,35.63 +27,5,3,20.67,122,3,29.64 +28,4,3,21.11,124,3,32.49 +29,5,3,32.59,131,3,31.71 +30,6,3,33.49,127,3,27.85 +31,5,3,21.81,127,3,21.79 +32,5,3,25.16,125,3,22.17 +33,5,3,20.27,122,3,22.54 +34,5,3,16.23,118,3,24.23 +35,6,3,15.01,117,3,20.54 +36,6,3,14.84,114,3,18.34 +37,6,3,12.57,114,3,18.48 +38,6,3,12.76,113,3,22.77 +39,5,3,13.56,108,3,21.46 +40,5,3,12.73,90,3,21.64 +41,6,3,12.76,115,3,18.44 +42,5,3,12.16,115,3,19.01 +43,5,3,11.32,115,3,17.96 +44,5,3,10.94,88,3,22.03 +45,5,3,10.85,118,3,22.92 +46,5,3,10.52,98,3,26.26 +47,5,3,9.12,108,3,16.03 +48,5,3,9.08,124,3,15.25 +49,5,3,8.91,124,3,15.04 +50,5,3,8.75,119,3,15.42 diff --git a/docs/advanced/figures/bench_fmg_vs_gamg_velocity_ksp.md b/docs/advanced/figures/bench_fmg_vs_gamg_velocity_ksp.md new file mode 100644 index 000000000..92157e4a9 --- /dev/null +++ b/docs/advanced/figures/bench_fmg_vs_gamg_velocity_ksp.md @@ -0,0 +1,121 @@ +--- +orphan: true +--- + +# Handover — FMG vs GAMG benchmark figure + adaptive-convection animation + +Produced 2026-06-11 from the boundary-slip / anisotropic-mover branch (PR #228). +Two assets: a solver-scaling **benchmark figure** (for `docs/advanced/`) and an +**animation** of the same run (for the PR description / a docs gallery). + +## Absolute locations (the worktree may differ — use these) + +- **Figure + CSV + this note** — currently *uncommitted* in the feature worktree: + `/Users/lmoresi/+Underworld/underworld3-pixi/.claude/worktrees/gmg-geometric-interp/docs/advanced/figures/` + Once PR #228 merges they live in the main checkout at + `/Users/lmoresi/+Underworld/underworld3-pixi/docs/advanced/figures/` + (repo-relative `docs/advanced/figures/` — stable in any checkout). +- **Animation, per-step PNGs, render scripts** — on disk, not in the repo: + `/Users/lmoresi/+Simulations/StagnantLid/anim_full_Ra1e7_dEta1e3_res32_R8_mode1/` + (frames + GIFs) and `/Users/lmoresi/+Simulations/StagnantLid/` (`_render_clean.py`, + `_render_watch.py`). +- **Harness** — `scripts/stagnant_lid_adapt_loop.py` in the repo (currently + `/Users/lmoresi/+Underworld/underworld3-pixi/.claude/worktrees/gmg-geometric-interp/scripts/stagnant_lid_adapt_loop.py`). + +The repo-relative paths used in the docs-reference snippets below are intentional +— image references resolve relative to the docs page, not the filesystem. + +--- + +## 1. Benchmark figure — `bench_fmg_vs_gamg_velocity_ksp.svg` + +**Files (here, `docs/advanced/figures/`):** +- `bench_fmg_vs_gamg_velocity_ksp.svg` — two-panel figure (text is `svg.fonttype=none`, so selectable/scalable). +- `bench_fmg_vs_gamg_velocity_ksp.csv` — the underlying per-step data: + `adapt_step, fmg_vel_ksp_its, fmg_snes_reason, fmg_t_stokes_s, gamg_vel_ksp_its, gamg_snes_reason, gamg_t_stokes_s`. + +**What it shows:** +- *Top* — inner velocity-block KSP iterations vs adaptation step. FMG (geometric full multigrid) holds a **mesh-independent ~5**; GAMG (algebraic) is **volatile ~64–131 (≈22×)** and does **not** cliff at R=8 over 50 steps. +- *Bottom* — Stokes-solve wall time. GAMG is only **~1.7× FMG**; FMG even spikes to match GAMG around the hardest adapts (steps 25–32) — that spike is the **cold-start Stokes solve** after an adapt (common to both engines), not the preconditioner. +- The **outer Schur-complement KSP converges in 1 iteration for both** — the entire difference lives in the inner velocity block. + +**Suggested caption:** +> **Velocity-block solver scaling under adaptive remeshing.** Inner velocity-block +> KSP iterations (top) and Stokes-solve wall time (bottom) versus adaptation step +> for a Ra = 10⁷, Δη = 10³ annulus convection model adapted every timestep with the +> MMPDE mover (res 32, resolution-ratio R = 8, np = 5). Geometric full multigrid +> (FMG) keeps a mesh-independent ≈ 5 inner iterations as the cells stretch, where +> algebraic multigrid (GAMG) runs a volatile ≈ 64–131 (≈ 22×) without cliffing at +> this anisotropy. The wall-clock gap is only ≈ 1.7×: each GAMG V-cycle is far +> cheaper than an FMG F-cycle, and the cold-start Stokes solve after each adapt +> (common to both) dominates the time. The outer Schur KSP converges in a single +> iteration for both — the difference is entirely in the inner velocity block. The +> value of geometric FMG here is *predictability and mesh-independence* (which widen +> with problem size and multigrid levels), not a large raw speed-up. + +**Docs wiring** (Markdown/MyST under the benchmark table): +```markdown +![Velocity-block solver scaling under adaptive remeshing](figures/bench_fmg_vs_gamg_velocity_ksp.svg) +``` +(Typst equivalent: `#image("figures/bench_fmg_vs_gamg_velocity_ksp.svg")`.) + +--- + +## 2. Animation companion — adaptive convection GIF + +**File (NOT in the repo — multi-MB — in the simulation directory):** +`/Users/lmoresi/+Simulations/StagnantLid/anim_full_Ra1e7_dEta1e3_res32_R8_mode1/` +- `anim_clean_Ra1e7_dEta1e3_res32_R8_mode1.gif` — **docs/PR size**: 8 MB, 63 frames (every 4th step), 480 px, clean style (T field + adaptive-mesh edges; no scalebar, streamlines, or text). +- `anim_clean_…_HQ.gif` — 24 MB, 125 frames, 600 px. +- `frame_clean_step0001…0250.png` — all 250 clean PNGs, kept for re-cutting cadence/size. + +**What it shows:** the 250-step run — a degree-1 perturbation breaking symmetry into +vigorous, irregular stagnant-lid convection (Nu 1 → 4.15), with the mesh +redistributing every step to track the inner thermal boundary layer and the plumes. +It is the *same run* whose Stokes solves stay flat under FMG in the figure above. + +**Suggested caption:** +> Adaptive-mesh convection in an annulus (Ra = 10⁷, Δη = 10³): temperature with the +> MMPDE-adapted mesh redistributing every timestep to follow the thermal boundary +> layer and plumes, mesh-owned tangent-slip keeping boundary nodes on the curved +> surface. 250 steps, np = 5 — the same run whose velocity solves stay +> mesh-independent under geometric FMG (see the benchmark figure). + +**Placement:** for the PR description, drag-drop the 8 MB GIF into the #228 comment +box via the GitHub web UI (under the 10 MB limit; `gh` can't upload binary +attachments). For docs, host/link it or commit a lighter cut — don't commit the +multi-MB GIF to the repo. + +--- + +## Reproduction + +Run (boundary-slip / mover branch; `scripts/stagnant_lid_adapt_loop.py`): +```bash +REFINE=2 PCVEL=gmg MG_TYPE=full MOVER=mmpde MMPDE_ACCEL=cg MOVER_SLIP=ring \ +mpirun -np 5 python scripts/stagnant_lid_adapt_loop.py \ + --from-perturbation --pert-mode 1 --Ra 1e7 --delta-eta 1e3 \ + --dt-cell-percentile 50 --skip-threshold 99 --adapt-every 1 \ + --resolution-ratio 8 --n-steps 250 --res 32 --snapshot-every 1 --out-tag +``` +- **FMG vs GAMG** (figure data): same config; `PCVEL=gmg MG_TYPE=full` (FMG) vs + `PCVEL=amg` (GAMG), 50 steps. Inner velocity iterations captured with a temporary + `FMG_DIAG` env-gated diagnostic in the harness + (`stokes.snes.getKSP().getPC().getFieldSplitSubKSP()[0].getIterationNumber()`, + plus `len(mesh.dm_hierarchy)` and the converged reason) — reverted after the runs. +- **Prerequisite** for FMG: the mesh must be built with `refinement` (so + `len(mesh.dm_hierarchy) > 1`) and adapted with a *coordinate-deforming* mover + (MMPDE preserves topology, so the hierarchy survives — confirmed `levels=3` on all + 50 steps). A checkpoint-resumed mesh has **no** hierarchy (it isn't persisted), so + `--resume` cannot be used for the FMG arm. +- **Render:** `_render_clean.py ` (clean — no scalebar/streamlines) or + `_render_watch.py ` (default — with scalebar + streamlines), both in + `/Users/lmoresi/+Simulations/StagnantLid/`. GIF assembled with Pillow (crop white margin → + resize → sample frames). + +## Status +- SVG + CSV + this note: **uncommitted** in `docs/advanced/figures/` — commit them + with the docs reference. +- The harness is at its committed state (the `FMG_DIAG` diagnostic was temporary). +- PR #228 (movers consume `mesh.boundary_slip` + the anisotropic-mover work) is + green and mergeable; this benchmark is its "why this works" companion. diff --git a/docs/advanced/figures/bench_fmg_vs_gamg_velocity_ksp.svg b/docs/advanced/figures/bench_fmg_vs_gamg_velocity_ksp.svg new file mode 100644 index 000000000..c670697c8 --- /dev/null +++ b/docs/advanced/figures/bench_fmg_vs_gamg_velocity_ksp.svg @@ -0,0 +1,946 @@ + + + + + + + + 2026-06-11T14:33:58.354179 + image/svg+xml + + + Matplotlib v3.10.8, https://matplotlib.org/ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 0 + + + + + + + + + + + + + 20 + + + + + + + + + + + + + 40 + + + + + + + + + + + + + 60 + + + + + + + + + + + + + 80 + + + + + + + + + + + + + 100 + + + + + + + + + + + + + 120 + + + + + + + + + + + + + 140 + + + + inner velocity-block + KSP iterations + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + GAMG: ~64–131 iters, volatile (no cliff at R=8) + + + FMG: flat ~5 (mesh-independent) + + + Velocity-block solver scaling under adaptive remeshing + + + + + + + + + + GAMG (algebraic) + + + + + + + + + FMG (geometric) + + + + + + + + + + + + + + + + + + + 10 + + + + + + + + + + + + + 20 + + + + + + + + + + + + + 30 + + + + + + + + + + + + + 40 + + + + + + + + + + + + + 50 + + + + adaptation step (anisotropy sharpens →) + + + + + + + + + + + + + + 0 + + + + + + + + + + + + + 10 + + + + + + + + + + + + + 20 + + + + + + + + + + + + + 30 + + + + + + + + + + + + + 40 + + + + Stokes solve + wall time (s) + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + wall: GAMG only ~1.7× FMG (cold-start-dominated; dilutes the 22× iteration gap) + + + + Annulus · Ra=1e7 · Δη=1e3 · res 32 · R=8 · mmpde + tangent-slip · np=5 | outer Schur KSP = 1 for both + + + + + + + + + + + diff --git a/docs/advanced/index.md b/docs/advanced/index.md index 48fab6fa6..252439272 100644 --- a/docs/advanced/index.md +++ b/docs/advanced/index.md @@ -18,6 +18,12 @@ Profile, optimize, and scale your simulations. **[→ Performance Guide](performance.md)** +### Multigrid Preconditioning (FMG vs GAMG) +Robust, anisotropy-tolerant solves on adapted meshes — build a mesh with +`refinement` and the solver uses geometric Full Multigrid automatically. + +**[→ Multigrid Preconditioning](multigrid-preconditioning.md)** + ### Complex Rheologies Implement advanced material models and constitutive laws. @@ -86,6 +92,7 @@ Ready to contribute to Underworld3? parallel-computing performance +multigrid-preconditioning complex-rheologies vep-transverse-isotropy-faults custom-meshes diff --git a/docs/advanced/multigrid-preconditioning.md b/docs/advanced/multigrid-preconditioning.md new file mode 100644 index 000000000..fd8369582 --- /dev/null +++ b/docs/advanced/multigrid-preconditioning.md @@ -0,0 +1,174 @@ +--- +title: "Multigrid Preconditioning (FMG vs GAMG)" +--- + +# Multigrid Preconditioning: FMG vs GAMG + +Underworld3 solvers are preconditioned with multigrid. There are two flavours, +and choosing the right one — or letting the solver choose for you — makes a large +difference to robustness and cost, especially on adapted or anisotropic meshes. + +- **GAMG** (*algebraic* multigrid) builds its coarse levels from the assembled + operator's connection graph. It is general and needs no mesh hierarchy, but it + is **sensitive to anisotropy**: stretched cells (exactly what mesh adaptation + produces) degrade the aggregates and the iteration count can cliff. +- **FMG** (geometric *Full Multigrid*) builds its coarse levels from a genuine + **mesh refinement hierarchy**. Because the hierarchy is geometric, it is + **inherently robust to anisotropy** — the coarse spaces are correct regardless + of how the operator is stretched. + +## The one-liner: build a mesh with `refinement` + +Geometric multigrid needs a refinement hierarchy. Build one by passing +`refinement=N` to any mesh constructor: + +```python +import underworld3 as uw + +# refinement=2 -> a 3-level geometric hierarchy (coarse -> medium -> fine) +mesh = uw.meshing.Annulus(radiusInner=0.5, radiusOuter=1.0, + cellSize=0.1, refinement=2) + +stokes = uw.systems.Stokes(mesh) +# ... constitutive model, body force, boundary conditions ... +stokes.solve() # velocity block is preconditioned with geometric FMG +``` + +That is the whole story for the common case. When the mesh carries a hierarchy, +the solver **automatically** uses geometric Full Multigrid on the velocity block +(for Stokes) or the top-level preconditioner (for scalar/vector solvers). When it +does not, the solver falls back to GAMG. You do not need to set any PETSc options. + +## The `preconditioner` knob + +Every Stokes / scalar / vector solver exposes a `preconditioner` property: + +```python +stokes.preconditioner = "auto" # default: FMG if the mesh has a hierarchy, else GAMG +stokes.preconditioner = "fmg" # force geometric multigrid (warns + falls back if no hierarchy) +stokes.preconditioner = "gamg" # force algebraic multigrid (the historical default) +``` + +`"mg"` is accepted as an alias for `"fmg"`. + +```{note} +`"auto"` is conservative. It only ever *adds* geometric multigrid on top of an +untouched default — it never rewrites a preconditioner you configured yourself. If +you set `pc_type` directly through `solver.petsc_options[...]`, or a solver applies +its own tuned options internally, `"auto"` leaves those settings alone. +``` + +## What the FMG bundle actually sets + +For reference (you should rarely need to set these by hand), selecting geometric +multigrid is equivalent to the following options on the relevant block +(`fieldsplit_velocity_` for Stokes, top-level for scalar/vector): + +```python +pc_type = "mg" +pc_mg_type = "full" # Full Multigrid (F-cycle) +pc_mg_galerkin = "both" # RAP coarse operators (see below) +mg_levels_ksp_type = "chebyshev" +mg_levels_pc_type = "sor" +mg_levels_ksp_max_it = 4 +mg_coarse_pc_type = "lu" # direct coarse solve +``` + +**Why Galerkin coarse operators?** Underworld3 does not install +residual/Jacobian callbacks on the coarse DMs, so the coarse operators must be +formed by Galerkin projection ($R A P$) from the fine operator rather than by +re-discretising on the coarse mesh. `pc_mg_galerkin = both` does this. + +## Hierarchy and mesh adaptation + +This is the key reason to prefer FMG when you adapt: + +- The coordinate-deforming adaptation movers (Winslow / anisotropic / + `OT_adapt` / `follow_metric`) **preserve mesh topology**. The refinement + hierarchy survives them, so geometric multigrid keeps working as the mesh + deforms — precisely where GAMG struggles with the resulting anisotropy. +- A **true remesh** (a topology change) collapses the hierarchy to a single + level. In `"auto"` mode the solver detects this and transparently reverts to + GAMG on the next solve, so nothing breaks — you simply lose the geometric path + until a hierarchy is available again. + +## Parallel coarse solve + +The default coarse solver is a serial direct solve (`mg_coarse_pc_type = "lu"`), +which is the fast, simple choice in serial and for modest core counts. For large +parallel partitions, replicate or gather the (small) coarse grid instead: + +```python +# redundant: copy the coarse grid to every rank and solve it there +stokes.petsc_options["fieldsplit_velocity_mg_coarse_pc_type"] = "redundant" +stokes.petsc_options["fieldsplit_velocity_mg_coarse_redundant_pc_type"] = "lu" + +# or telescope onto a sub-communicator for very large runs +# stokes.petsc_options["fieldsplit_velocity_mg_coarse_pc_type"] = "telescope" +``` + +Because `"auto"` does not overwrite options you set yourself, you can keep +`preconditioner = "auto"` and just override the coarse solver as above. + +## Benchmark: FMG vs GAMG on a deforming adaptive mesh + +The payoff is clearest on an aggressively adapted mesh. The figure below is a +Stokes convection run (annulus, Ra = 10⁷, Δη = 10³, res 32, resolution-ratio +R = 8, mode-1, `np = 5`) whose mesh is continuously deformed by the MMPDE +coordinate mover, adapted every timestep. The geometric hierarchy (3 levels) +survives every step — topology is preserved — and both engines converge cleanly +(`snes` reason 3) throughout. FMG (`PCVEL=gmg MG_TYPE=full`) and GAMG +(`PCVEL=amg`) ran the *identical* adapted-mesh sequence over 50 steps. + +```{figure} figures/bench_fmg_vs_gamg_velocity_ksp.svg +:alt: Inner velocity-block KSP iterations and Stokes-solve wall time versus adaptation step, for FMG and GAMG. + +**Velocity-block solver scaling under adaptive remeshing.** Inner velocity-block +KSP iterations (top) and Stokes-solve wall time (bottom) versus adaptation step. +Geometric full multigrid (FMG) keeps a mesh-independent ≈ 5 inner iterations as +the cells stretch, where algebraic multigrid (GAMG) runs a volatile ≈ 64–131 +(≈ 23×) without cliffing at this anisotropy. The wall-clock gap is only ≈ 1.8×: +each GAMG V-cycle is far cheaper than an FMG F-cycle, and the cold-start Stokes +solve after each adapt (common to both) dominates the time. The outer Schur KSP +converges in one iteration for both — the difference is entirely in the inner +velocity block. +``` + +The metric that matters is the **inner velocity-block KSP iteration count** — the +outer Schur KSP is one iteration for both engines, so the entire difference lives +in this inner block: + +- **FMG** holds a **mesh-independent ≈ 5** iterations (3–6) as the anisotropy + sharpens — the geometric-MG signature. +- **GAMG** runs a **volatile ≈ 64–131** (median ≈ 114, so ≈ 23×) as the algebraic + aggregates cope with the stretched cells. It does *not* cliff at R = 8 over + these 50 steps, but it is unpredictable from step to step. + +The wall-clock gap is only **≈ 1.8×**, not 23×: a single GAMG V-cycle is far +cheaper than an FMG F-cycle, and the cold-start Stokes solve after each adapt +(common to both engines) dominates the per-step time. So the value of geometric +FMG here is **predictability and mesh-independence** — properties that widen with +problem size and multigrid depth — rather than a large raw speed-up. The iteration +gap is also where GAMG eventually loses robustness on harder problems (higher +anisotropy, more levels). + +The per-step data is in +[`figures/bench_fmg_vs_gamg_velocity_ksp.csv`](figures/bench_fmg_vs_gamg_velocity_ksp.csv); +the reproduction command and full provenance are in the companion note +`figures/bench_fmg_vs_gamg_velocity_ksp.md`. + +## When to use which + +| Situation | Recommended | +|-----------|-------------| +| Adapted / deformed / anisotropic meshes | **FMG** (build with `refinement`) | +| Uniform mesh, want mesh-independent iteration counts | **FMG** | +| No refinement hierarchy available / quick prototype | GAMG (automatic fallback) | +| Reproducing historical results exactly | `preconditioner = "gamg"` | + +## See also + +- {doc}`mesh-adaptation` — the movers whose anisotropy FMG handles gracefully +- {doc}`performance` — profiling and scaling +- `docs/developer/design/solver-strategies-catalogue.md` — the solver-strategy + design notes (developer-facing) diff --git a/docs/developer/design/solver-strategies-catalogue.md b/docs/developer/design/solver-strategies-catalogue.md index 3558b36dc..f23a338f8 100644 --- a/docs/developer/design/solver-strategies-catalogue.md +++ b/docs/developer/design/solver-strategies-catalogue.md @@ -308,10 +308,34 @@ note (cf. `snes-atol-convergence-scale.md`) before implementation. anisotropy-robust (the hierarchy is built geometrically, not from the operator's connection graph). -**Status:** UW3 has `dm_hierarchy` / `refineHierarchy` -infrastructure. Pairs naturally with the error-estimator design -arc — the same multi-level structure yields both the -anisotropy-robust PC and the absolute error indicator. +**Status:** **Landed.** A mesh built with `refinement >= 1` carries a +`dm_hierarchy`, and the solvers now switch to geometric Full Multigrid +on it automatically. The user-facing control is the `preconditioner` +property (`"auto"` | `"fmg"` | `"gamg"`) on the Stokes / scalar / vector +solvers; `"auto"` (the default) selects FMG when a hierarchy is present +and falls back to GAMG otherwise. + +Implementation: `SolverBaseClass._apply_preconditioner_options()` in +`petsc_generic_snes_solvers.pyx`, invoked from `_build` so `"auto"` +re-resolves against the current mesh (a true remesh collapses the +hierarchy → automatic GAMG fallback). The auto path is deliberately +conservative — it only *adds* geometric MG on top of an untouched +default and never rewrites pc options a solver/user configured directly +(e.g. the tuned GAMG in the OT/φ-Poisson smoother). User guide: +`docs/advanced/multigrid-preconditioning.md`. Tests: +`tests/test_1014_stokes_multigrid.py`. + +Benchmark-validated on a deforming adaptive mesh (annulus res32, R=8, +mode-1, np=5, MMPDE mover, 50 steps): the 3-level hierarchy survives +every step and the inner velocity-block KSP stays flat at ~5 iters under +FMG where GAMG is a volatile ~64-131 (~23×) without cliffing at this +anisotropy; wall-clock gap only ~1.8× (the cold-start Stokes solve, common +to both, dominates the time). The value is predictability/mesh-independence, +not raw speed. Figure + data in `docs/advanced/multigrid-preconditioning.md`. + +Still pairs naturally with the error-estimator design arc — the same +multi-level structure yields both the anisotropy-robust PC and the +absolute error indicator. ## Mesh-quality / `mesh.quality()` API diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index 73257c4ec..4917ede7c 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -67,6 +67,203 @@ class SolverBaseClass(uw_object): self._pressure_is = None self._subdict = {} + # Preconditioner selection — see the `preconditioner` property. + # `_pc_option_prefix` is set by subclasses that participate in the easy + # FMG/GAMG switch ("" for scalar/vector, "fieldsplit_velocity_" for + # Stokes); None means "this solver manages its own PC options" and the + # helper below is a no-op. + self._preconditioner = "auto" + self._pc_option_prefix = None + # The pc_type value this helper last managed. Subclasses that opt in set + # their __init__ default ("gamg"); used in "auto" mode to tell an + # untouched framework default (eligible for FMG upgrade) apart from an + # explicit user override of pc_type, which must be respected. + self._pc_managed_value = "gamg" + # Latches once the user (or harness) is seen to have set the PC options + # themselves. _apply_preconditioner_options() runs on EVERY _build (so + # "auto" can re-resolve after a remesh), and without this latch the + # override detection only holds on the first build: once we adopted the + # user's pc_type as our managed value, a later rebuild could no longer + # tell "user set mg" from "we set mg" and would clobber their tuned + # smoother / coarse-solver options with the framework FMG bundle. + self._pc_user_override = False + + @property + def preconditioner(self): + """Preconditioner selection for the (velocity) block. + + One of: + + - ``"auto"`` (default) — use geometric Full Multigrid (FMG) when the + mesh carries a genuine refinement hierarchy + (``len(mesh.dm_hierarchy) > 1``, i.e. built with ``refinement >= 1``), + otherwise fall back to algebraic multigrid (GAMG). + - ``"fmg"`` (alias ``"mg"``) — force geometric multigrid. Requires a + refinement hierarchy; warns and falls back to GAMG if none exists. + - ``"gamg"`` — force algebraic multigrid (the historical default). + + Geometric multigrid is inherently robust to mesh anisotropy (it is built + from the geometric refinement hierarchy, not the operator connection + graph), which makes it the preferred choice on adapted/deformed meshes. + The hierarchy survives the coordinate-deforming adaptation movers and + only collapses under a true remesh — in which case ``"auto"`` + transparently reverts to GAMG. + + For Stokes this governs the velocity fieldsplit block; for scalar/vector + solvers it governs the top-level preconditioner. Other solvers ignore it. + """ + return self._preconditioner + + @preconditioner.setter + def preconditioner(self, value): + choice = str(value).lower() + if choice == "mg": + choice = "fmg" + if choice not in ("auto", "fmg", "gamg"): + raise ValueError( + f"preconditioner must be 'auto', 'fmg', or 'gamg' (got {value!r})" + ) + self._preconditioner = choice + # Force a full rebuild so the new option bundle is pushed to PETSc. + self.is_setup = False + + def _apply_preconditioner_options(self): + """Push the PETSc option bundle implied by ``self.preconditioner``. + + Called from ``_build`` so that ``"auto"`` re-resolves against the + *current* mesh — e.g. after a remesh that collapses the refinement + hierarchy. Solvers opt in by setting ``self._pc_option_prefix`` (``""`` + or ``"fieldsplit_velocity_"``); the default (None) makes this a no-op. + """ + prefix = self._pc_option_prefix + if prefix is None: + return + + opts = self.petsc_options + + # The mesh is the source of truth for hierarchy depth: the solver's own + # dm_hierarchy is a clone of it and may not be populated this early. + n_levels = len(getattr(self.mesh, "dm_hierarchy", []) or []) + + if self._preconditioner == "auto": + # Auto mode is deliberately conservative: it only ever *adds* + # geometric multigrid, never rewrites an existing GAMG/other + # configuration (internal utility solvers — mesh smoothing, OT, + # projections — set their own carefully tuned pc options that must + # be left intact). + if self._pc_user_override: + # User/harness owns the PC options — never re-apply our bundle, + # even across remesh rebuilds (see _pc_user_override init note). + return + current = opts.getString(f"{prefix}pc_type") + if current and current != self._pc_managed_value: + # Preconditioner was set explicitly elsewhere — respect it, and + # latch so later rebuilds keep respecting it. + self._pc_user_override = True + self._pc_managed_value = current + return + if n_levels > 1: + want_fmg = True + else: + # No hierarchy. Only act to revert a previous FMG upgrade of + # ours (e.g. after a remesh collapsed the hierarchy); otherwise + # leave the existing default/tuned options untouched. + if self._pc_managed_value != "mg": + return + want_fmg = False + elif self._preconditioner == "fmg": + want_fmg = n_levels > 1 + else: # "gamg" — explicit, always applied + want_fmg = False + + if want_fmg: + # Geometric Full Multigrid on the refinement hierarchy. Galerkin + # (RAP) coarse operators are required because UW3 does not install + # residual/Jacobian callbacks on the coarse DMs. + opts[f"{prefix}pc_type"] = "mg" + opts[f"{prefix}pc_mg_type"] = "full" # FMG (F-cycle) + opts[f"{prefix}pc_mg_galerkin"] = "both" # RAP coarse operators + # richardson+sor (not chebyshev): chebyshev needs eigenvalue + # estimates of the smoothed operator, which are fragile on the + # indefinite / variable-viscosity Stokes velocity block and diverge; + # richardson+sor is the benchmark-validated, mesh-independent choice. + opts[f"{prefix}mg_levels_ksp_type"] = "richardson" + opts[f"{prefix}mg_levels_pc_type"] = "sor" + opts[f"{prefix}mg_levels_ksp_max_it"] = 4 + opts[f"{prefix}mg_levels_ksp_converged_maxits"] = None + # redundant+lu, not bare lu: a bare serial LU cannot factor a + # distributed coarse matrix and fails at np>1 (DIVERGED_LINEAR_SOLVE + # after 0 iterations). redundant gathers the (small) coarse system to + # one rank and is identical to lu in serial — so it is np-safe by + # default without surprising small-np users. + opts[f"{prefix}mg_coarse_pc_type"] = "redundant" + opts[f"{prefix}mg_coarse_redundant_pc_type"] = "lu" + # Clear stale GAMG-only keys so toggling back and forth is clean. + for key in ("pc_gamg_type", "pc_gamg_repartition", "pc_gamg_agg_nsmooths"): + opts.delValue(f"{prefix}{key}") + self._pc_managed_value = "mg" + else: + if self._preconditioner == "fmg" and uw.mpi.rank == 0: + import warnings + warnings.warn( + f"[{self.name}] preconditioner='fmg' requested but the mesh " + f"has no refinement hierarchy; falling back to GAMG. Build the " + f"mesh with refinement >= 1 to enable geometric multigrid.", + stacklevel=2, + ) + opts[f"{prefix}pc_type"] = "gamg" + opts[f"{prefix}pc_gamg_type"] = "agg" + opts[f"{prefix}pc_gamg_repartition"] = True + opts[f"{prefix}pc_mg_type"] = "additive" + opts[f"{prefix}pc_gamg_agg_nsmooths"] = 2 + opts[f"{prefix}mg_levels_ksp_max_it"] = 3 + opts[f"{prefix}mg_levels_ksp_converged_maxits"] = None + # Clear stale geometric-MG-only keys. + for key in ("pc_mg_galerkin", "mg_levels_ksp_type", + "mg_levels_pc_type", "mg_coarse_pc_type", + "mg_coarse_redundant_pc_type"): + opts.delValue(f"{prefix}{key}") + self._pc_managed_value = "gamg" + + def _enforce_galerkin_for_geometric_mg(self): + """Geometric multigrid in UW3 REQUIRES Galerkin (RAP) coarse operators. + + UW3 installs no residual/Jacobian callbacks on the coarse DMs of the + refinement hierarchy, so PETSc cannot re-discretise the operator there. + With Galerkin RAP the coarse operators are assembled as R*A*P from the + fine operator and the coarse DMs are used only for interpolation — which + is the only mode that works. If ``pc_type=mg`` is selected on a block we + manage but Galerkin is unset (or explicitly ``none``), PETSc instead + tries to build coarse operators itself and fails cryptically: PETSc + error 73 ("KSPSetDM without ComputeOperators") in serial, or + ``DMCoarsen -> DMAdaptMetric -> ParMmg`` (which is 3D-only) in parallel. + + So rather than let users trip those, force ``pc_mg_galerkin=both`` + whenever a managed block uses geometric MG, regardless of who selected + it (user, harness, or our own auto/fmg bundle). A bare ``=None`` flag + already engages RAP (reads back as ``''`` but is *present*); only a + genuinely-unset or ``none`` value is overridden. + """ + prefix = self._pc_option_prefix + if prefix is None: + return + opts = self.petsc_options + if opts.getString(f"{prefix}pc_type") != "mg": + return + gkey = f"{prefix}pc_mg_galerkin" + if (not opts.hasName(gkey)) or opts.getString(gkey) == "none": + if uw.mpi.rank == 0: + import warnings + warnings.warn( + f"[{self.name}] geometric multigrid (pc_type=mg) requires " + f"Galerkin coarse operators in UW3 — forcing {gkey}=both. " + f"(UW3 installs no coarse-DM operator callbacks, so coarse " + f"re-discretisation is unsupported and fails as PETSc error " + f"73 / ParMmg.)", + stacklevel=2, + ) + opts[gkey] = "both" + def _check_expression_meshes(self): """Check that all MeshVariable symbols in solver expressions belong to this solver's mesh. @@ -597,6 +794,13 @@ class SolverBaseClass(uw_object): if self.is_setup: return + # Resolve the preconditioner choice (auto/fmg/gamg) against the current + # mesh hierarchy and push the option bundle before the per-class + # _setup_solver runs snes.setFromOptions(). No-op unless the solver + # opted in via _pc_option_prefix (Stokes / scalar / vector). + self._apply_preconditioner_options() + self._enforce_galerkin_for_geometric_mg() + # === Fast path 1: constants-only change === # If JIT cache key matches, the compiled code is identical — only # constant values (dt_elastic, scalar viscosity, BDF coeffs) differ. @@ -1616,6 +1820,13 @@ class SNES_Scalar(SolverBaseClass): # ROBUST and general GAMG, heavy-duty solvers in the suite + # Participate in the auto FMG/GAMG switch (see the `preconditioner` + # property). The pc/mg keys below are the GAMG default; + # _apply_preconditioner_options() upgrades this block to geometric FMG + # at build time when the mesh carries a refinement hierarchy, and + # otherwise leaves these defaults in place. + self._pc_option_prefix = "" + self.petsc_options["snes_type"] = "newtonls" self.petsc_options["ksp_type"] = "gmres" self.petsc_options["pc_type"] = "gamg" @@ -2469,6 +2680,13 @@ class SNES_Vector(SolverBaseClass): # options = PETSc.Options() # options["dm_adaptor"]= "pragmatic" + # Participate in the auto FMG/GAMG switch (see the `preconditioner` + # property). The pc/mg keys below are the GAMG default; + # _apply_preconditioner_options() upgrades this block to geometric FMG + # at build time when the mesh carries a refinement hierarchy, and + # otherwise leaves these defaults in place. + self._pc_option_prefix = "" + # Here we can set some defaults for this set of KSP / SNES solvers self.petsc_options["snes_type"] = "newtonls" self.petsc_options["ksp_rtol"] = 1.0e-3 @@ -4197,6 +4415,13 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): self._tolerance = 1.0e-4 self._strategy = "default" + # Participate in the auto FMG/GAMG switch on the velocity fieldsplit + # block (see the `preconditioner` property). The velocity pc/mg keys + # set below are the GAMG default; _apply_preconditioner_options() + # upgrades the velocity block to geometric FMG at build time when the + # mesh carries a refinement hierarchy, and otherwise leaves them. + self._pc_option_prefix = "fieldsplit_velocity_" + self.petsc_options["snes_rtol"] = self._tolerance self.petsc_options["snes_ksp_ew"] = None self.petsc_options["snes_ksp_ew_version"] = 3 diff --git a/tests/test_1014_stokes_multigrid.py b/tests/test_1014_stokes_multigrid.py new file mode 100644 index 000000000..e29c2be99 --- /dev/null +++ b/tests/test_1014_stokes_multigrid.py @@ -0,0 +1,156 @@ +import pytest +import sympy +import underworld3 as uw + +# Solves small Stokes/Poisson systems to exercise the auto FMG/GAMG switch. +pytestmark = [pytest.mark.level_2, pytest.mark.tier_a] + + +# A mesh built with refinement carries a real dm_hierarchy (FMG-capable); +# a plain mesh has a single level and must fall back to GAMG. +mesh_refined = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), + cellSize=0.25, refinement=2, qdegree=2, +) +mesh_plain = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), + cellSize=0.2, qdegree=2, +) + + +def _make_stokes(mesh): + x, y = mesh.X + stokes = uw.systems.Stokes(mesh) + stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel + stokes.constitutive_model.Parameters.shear_viscosity_0 = 1 + stokes.bodyforce = sympy.Matrix([0, 1.0e2 * x]) + stokes.add_dirichlet_bc((0.0, 0.0), "Bottom") + stokes.add_dirichlet_bc((0.0, 0.0), "Top") + stokes.add_dirichlet_bc((0.0, None), "Left") + stokes.add_dirichlet_bc((0.0, None), "Right") + return stokes + + +def _vel_pc(stokes): + return stokes.petsc_options.getString("fieldsplit_velocity_pc_type") + + +def test_refinement_mesh_has_hierarchy(): + # Sanity: refinement=2 builds a 3-level hierarchy; plain mesh has one level. + assert len(mesh_refined.dm_hierarchy) > 1 + assert len(mesh_plain.dm_hierarchy) == 1 + + +def test_auto_selects_geometric_fmg_on_refinement_mesh(): + stokes = _make_stokes(mesh_refined) + assert stokes.preconditioner == "auto" # default + stokes.solve() + assert _vel_pc(stokes) == "mg" + assert stokes.snes.getConvergedReason() > 0 + + +def test_auto_falls_back_to_gamg_without_hierarchy(): + stokes = _make_stokes(mesh_plain) + stokes.solve() + assert _vel_pc(stokes) == "gamg" + assert stokes.snes.getConvergedReason() > 0 + + +def test_explicit_gamg_override(): + stokes = _make_stokes(mesh_refined) + stokes.preconditioner = "gamg" + stokes.solve() + assert _vel_pc(stokes) == "gamg" + assert stokes.snes.getConvergedReason() > 0 + + +def test_toggle_back_to_fmg(): + # Switching gamg -> fmg must cleanly re-engage geometric multigrid. + stokes = _make_stokes(mesh_refined) + stokes.preconditioner = "gamg" + stokes.solve() + stokes.preconditioner = "fmg" + stokes.solve() + assert _vel_pc(stokes) == "mg" + assert stokes.snes.getConvergedReason() > 0 + + +def test_fmg_without_hierarchy_falls_back(): + # Asking for fmg on a mesh with no hierarchy warns and uses GAMG. + stokes = _make_stokes(mesh_plain) + stokes.preconditioner = "fmg" + with pytest.warns(UserWarning, match="falling back to GAMG"): + stokes.solve() + assert _vel_pc(stokes) == "gamg" + assert stokes.snes.getConvergedReason() > 0 + + +def test_invalid_preconditioner_raises(): + stokes = _make_stokes(mesh_plain) + with pytest.raises(ValueError): + stokes.preconditioner = "wibble" + + +def test_scalar_poisson_auto_geometric_mg(): + poisson = uw.systems.Poisson(mesh_refined) + poisson.constitutive_model = uw.constitutive_models.DiffusionModel + poisson.constitutive_model.Parameters.diffusivity = 1 + poisson.f = 0.0 + poisson.add_dirichlet_bc(0.0, "Bottom") + poisson.add_dirichlet_bc(1.0, "Top") + poisson.solve() + assert poisson.petsc_options.getString("pc_type") == "mg" + assert poisson.snes.getConvergedReason() > 0 + + +def test_explicit_velocity_options_survive_rebuild(): + # Regression: in "auto" mode the helper must NOT clobber user-set + # velocity-block PC options when it re-runs. _apply_preconditioner_options() + # is called on EVERY _build (so "auto" re-resolves after a remesh); a + # mesh-mover deform triggers exactly such a rebuild. Previously the override + # was respected only on the first call — the second overwrote the user's + # tuned coarse-solver / smoother with the framework FMG bundle. Drive the + # helper directly (twice) to test the option latch in isolation. + stokes = _make_stokes(mesh_refined) + vp = "fieldsplit_velocity_" + stokes.petsc_options[vp + "pc_type"] = "mg" + stokes.petsc_options[vp + "mg_coarse_pc_type"] = "redundant" + stokes.petsc_options[vp + "mg_coarse_redundant_pc_type"] = "lu" + stokes.petsc_options[vp + "mg_levels_ksp_type"] = "richardson" + + stokes._apply_preconditioner_options() # first build: adopt + respect + assert stokes.petsc_options.getString(vp + "mg_coarse_pc_type") == "redundant" + stokes._apply_preconditioner_options() # rebuild: must STILL respect them + assert stokes.petsc_options.getString(vp + "pc_type") == "mg" + assert stokes.petsc_options.getString(vp + "mg_coarse_pc_type") == "redundant" + assert stokes.petsc_options.getString(vp + "mg_levels_ksp_type") == "richardson" + + +def test_geometric_mg_without_galerkin_is_repaired(): + # UW3's geometric MG REQUIRES Galerkin RAP (no coarse-DM operator callbacks + # are installed, so PETSc cannot re-discretise on the coarse levels). A user + # who selects pc_type=mg but omits pc_mg_galerkin must be repaired (forced to + # "both") and warned — not left to fail as PETSc error 73 (serial) or + # DMCoarsen->ParMmg (parallel). + stokes = _make_stokes(mesh_refined) + vp = "fieldsplit_velocity_" + stokes.petsc_options[vp + "pc_type"] = "mg" + # deliberately DO NOT set pc_mg_galerkin + with pytest.warns(UserWarning, match="requires Galerkin"): + stokes.solve() + assert stokes.petsc_options.getString(vp + "pc_mg_galerkin") == "both" + assert stokes.snes.getConvergedReason() > 0 + + +def test_default_fmg_bundle_is_parallel_safe(): + # The property's OWN default FMG bundle must be usable at np>1 unaided: a + # parallel-safe coarse solver (redundant+lu, not bare serial lu) and a + # robust smoother (richardson+sor, not eigen-estimate-fragile chebyshev). + stokes = _make_stokes(mesh_refined) + stokes.preconditioner = "fmg" + stokes.solve() + vp = "fieldsplit_velocity_" + assert stokes.petsc_options.getString(vp + "mg_coarse_pc_type") == "redundant" + assert stokes.petsc_options.getString(vp + "mg_coarse_redundant_pc_type") == "lu" + assert stokes.petsc_options.getString(vp + "mg_levels_ksp_type") == "richardson" + assert stokes.snes.getConvergedReason() > 0