From 020352e0d9c28e00c04f351d1c2e10ad6e2f952d Mon Sep 17 00:00:00 2001 From: lmoresi Date: Thu, 11 Jun 2026 15:01:47 +1000 Subject: [PATCH 1/5] Add `preconditioner` property: automatic geometric FMG on refinement meshes MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit UW3 Stokes / scalar / vector solvers defaulted to GAMG (algebraic multigrid), which is sensitive to mesh anisotropy — exactly what mesh adaptation produces. Geometric Full Multigrid (FMG) is anisotropy-robust but required hand-setting ~8 cryptic `fieldsplit_velocity_*` PETSc options. Adds a `preconditioner` property to the SNES solver base class (inherited by Stokes, SNES_Scalar, SNES_Vector): - "auto" (default): geometric FMG when the mesh carries a refinement hierarchy (`len(mesh.dm_hierarchy) > 1`), otherwise GAMG. - "fmg" / "gamg": force either. So building a mesh with `refinement=N` and solving now uses FMG automatically — no PETSc options to set. The choice is resolved at build time (re-evaluated after a remesh), keyed by a per-class option prefix ("" for scalar/vector, "fieldsplit_velocity_" for Stokes). "auto" is deliberately conservative: it only adds geometric MG on top of an untouched framework default; it never rewrites pc options a solver or user set directly (e.g. the tuned GAMG in the OT/mesh-smoother, or an explicit `pc_type=lu`). Benchmark (annulus res32, R=8, mode-1, np=5, MMPDE mover, deforming mesh): the inner velocity-block KSP stays flat at ~4 iterations under FMG where GAMG climbs 75 -> ~147 as cells stretch (~25-35x the count, ~2x wall-clock). Tests: tests/test_1014_stokes_multigrid.py (8, tier_a). Docs: docs/advanced/multigrid-preconditioning.md. Underworld development team with AI support from Claude Code --- docs/advanced/index.md | 7 + docs/advanced/multigrid-preconditioning.md | 153 +++++++++++++++++ .../design/solver-strategies-catalogue.md | 30 +++- .../cython/petsc_generic_snes_solvers.pyx | 160 ++++++++++++++++++ tests/test_1014_stokes_multigrid.py | 102 +++++++++++ 5 files changed, 448 insertions(+), 4 deletions(-) create mode 100644 docs/advanced/multigrid-preconditioning.md create mode 100644 tests/test_1014_stokes_multigrid.py 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..98dafd457 --- /dev/null +++ b/docs/advanced/multigrid-preconditioning.md @@ -0,0 +1,153 @@ +--- +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 following is a +Stokes convection run (annulus, res ≈ 32, radius ratio 8, mode-1, `np=5`) whose +mesh is continuously deformed by the MMPDE coordinate mover. The geometric +hierarchy (3 levels) survives every step — topology is preserved — and the solve +converges cleanly (`snes` reason 3) throughout. Both configurations were run on +the *identical* adapted-mesh sequence. + +The metric that matters is the **inner velocity-block KSP iteration count** as the +cell anisotropy sharpens (the outer Schur KSP is 1 iteration in both cases, so the +entire cost lives in this inner block): + +| velocity-block KSP iters | step 1 (near-uniform) | steps 5–10 | steps 15–20 (sharp) | Stokes wall, step 20 | +|--------------------------|:---------------------:|:----------:|:-------------------:|:--------------------:| +| **FMG** (geometric) | 5 | ~4 | ~4 | 12.5 s | +| **GAMG** (algebraic) | 75 | ~85–115 | ~125–147 | 23.6 s | + +FMG stays **flat at ~4 iterations** — the textbook mesh-independent geometric-MG +signature — while GAMG climbs from 75 to ~147 (volatile and growing, ~25–35× FMG's +count and ~2× the wall-clock) as the algebraic aggregates degrade under the +stretched cells. A 250-step continuation confirms FMG holds the full horizon +without cliffing once the hierarchy is established. GAMG did not outright fail in +20 steps here, but the trend is unmistakable and on sharper problems it reaches +`DIVERGED_PC_FAILED`. + +## 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..33ee297c9 100644 --- a/docs/developer/design/solver-strategies-catalogue.md +++ b/docs/developer/design/solver-strategies-catalogue.md @@ -308,10 +308,32 @@ 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): the 3-level hierarchy survives every step +and the inner velocity-block KSP stays flat at ~4 iters under FMG where +GAMG climbs 75→~147 (~25-35×, ~2× wall) as cells stretch. Table 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..dbb347c2a 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -67,6 +67,139 @@ 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" + + @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). + current = opts.getString(f"{prefix}pc_type") + if current and current != self._pc_managed_value: + # Preconditioner was set explicitly elsewhere — respect it. + 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" + opts[f"{prefix}mg_levels_ksp_type"] = "chebyshev" + 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 + opts[f"{prefix}mg_coarse_pc_type"] = "lu" # direct coarse solve (serial) + # 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": + if uw.mpi.rank == 0: + print( + 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.", + flush=True, + ) + 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"): + opts.delValue(f"{prefix}{key}") + self._pc_managed_value = "gamg" + def _check_expression_meshes(self): """Check that all MeshVariable symbols in solver expressions belong to this solver's mesh. @@ -597,6 +730,12 @@ 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() + # === 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 +1755,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 +2615,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 +4350,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..907142f1a --- /dev/null +++ b/tests/test_1014_stokes_multigrid.py @@ -0,0 +1,102 @@ +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" + 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 From 108503098ced476bc5be288d0ce91d45434fa637 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Thu, 11 Jun 2026 16:53:52 +1000 Subject: [PATCH 2/5] docs(multigrid): add FMG vs GAMG benchmark figure, reconcile the numbers MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Replace the placeholder text table in the multigrid how-to with the two-panel benchmark figure (inner velocity-block KSP iterations + Stokes wall time vs adaptation step) and its per-step CSV + provenance note. Align the prose to the measured 50-step run: FMG holds a mesh-independent ~5 iterations; GAMG is a volatile ~64-131 (~23x) that does not cliff at R=8; the wall-clock gap is only ~1.8x because the cold-start Stokes solve after each adapt (common to both engines) dominates the time. The value of geometric FMG here is predictability and mesh-independence, not a large raw speed-up — this corrects the earlier overstated ~2x / climbing-to-failure framing. Figure produced on the anisotropic-mover branch (#228); it lives here with the how-to it illustrates. Underworld development team with AI support from Claude Code --- .../bench_fmg_vs_gamg_velocity_ksp.csv | 51 + .../figures/bench_fmg_vs_gamg_velocity_ksp.md | 117 +++ .../bench_fmg_vs_gamg_velocity_ksp.svg | 946 ++++++++++++++++++ docs/advanced/multigrid-preconditioning.md | 67 +- .../design/solver-strategies-catalogue.md | 10 +- 5 files changed, 1164 insertions(+), 27 deletions(-) create mode 100644 docs/advanced/figures/bench_fmg_vs_gamg_velocity_ksp.csv create mode 100644 docs/advanced/figures/bench_fmg_vs_gamg_velocity_ksp.md create mode 100644 docs/advanced/figures/bench_fmg_vs_gamg_velocity_ksp.svg 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..c68e54292 --- /dev/null +++ b/docs/advanced/figures/bench_fmg_vs_gamg_velocity_ksp.md @@ -0,0 +1,117 @@ +# 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/multigrid-preconditioning.md b/docs/advanced/multigrid-preconditioning.md index 98dafd457..fd8369582 100644 --- a/docs/advanced/multigrid-preconditioning.md +++ b/docs/advanced/multigrid-preconditioning.md @@ -112,29 +112,50 @@ Because `"auto"` does not overwrite options you set yourself, you can keep ## Benchmark: FMG vs GAMG on a deforming adaptive mesh -The payoff is clearest on an aggressively adapted mesh. The following is a -Stokes convection run (annulus, res ≈ 32, radius ratio 8, mode-1, `np=5`) whose -mesh is continuously deformed by the MMPDE coordinate mover. The geometric -hierarchy (3 levels) survives every step — topology is preserved — and the solve -converges cleanly (`snes` reason 3) throughout. Both configurations were run on -the *identical* adapted-mesh sequence. - -The metric that matters is the **inner velocity-block KSP iteration count** as the -cell anisotropy sharpens (the outer Schur KSP is 1 iteration in both cases, so the -entire cost lives in this inner block): - -| velocity-block KSP iters | step 1 (near-uniform) | steps 5–10 | steps 15–20 (sharp) | Stokes wall, step 20 | -|--------------------------|:---------------------:|:----------:|:-------------------:|:--------------------:| -| **FMG** (geometric) | 5 | ~4 | ~4 | 12.5 s | -| **GAMG** (algebraic) | 75 | ~85–115 | ~125–147 | 23.6 s | - -FMG stays **flat at ~4 iterations** — the textbook mesh-independent geometric-MG -signature — while GAMG climbs from 75 to ~147 (volatile and growing, ~25–35× FMG's -count and ~2× the wall-clock) as the algebraic aggregates degrade under the -stretched cells. A 250-step continuation confirms FMG holds the full horizon -without cliffing once the hierarchy is established. GAMG did not outright fail in -20 steps here, but the trend is unmistakable and on sharper problems it reaches -`DIVERGED_PC_FAILED`. +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 diff --git a/docs/developer/design/solver-strategies-catalogue.md b/docs/developer/design/solver-strategies-catalogue.md index 33ee297c9..f23a338f8 100644 --- a/docs/developer/design/solver-strategies-catalogue.md +++ b/docs/developer/design/solver-strategies-catalogue.md @@ -326,10 +326,12 @@ default and never rewrites pc options a solver/user configured directly `tests/test_1014_stokes_multigrid.py`. Benchmark-validated on a deforming adaptive mesh (annulus res32, R=8, -mode-1, np=5, MMPDE mover): the 3-level hierarchy survives every step -and the inner velocity-block KSP stays flat at ~4 iters under FMG where -GAMG climbs 75→~147 (~25-35×, ~2× wall) as cells stretch. Table in -`docs/advanced/multigrid-preconditioning.md`. +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 From 2e971fa9d008aab8e277c20120eef3e8798fa89c Mon Sep 17 00:00:00 2001 From: lmoresi Date: Thu, 11 Jun 2026 18:17:16 +1000 Subject: [PATCH 3/5] docs(multigrid): mark benchmark provenance note as an orphan page The figure handover note lives beside the SVG/CSV for provenance but is not a navigable doc page; `orphan: true` front-matter stops Sphinx warning that it is not in any toctree. `pixi run docs-build` is clean on the multigrid page (figure renders, CSV is a downloadable link). Underworld development team with AI support from Claude Code --- docs/advanced/figures/bench_fmg_vs_gamg_velocity_ksp.md | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/docs/advanced/figures/bench_fmg_vs_gamg_velocity_ksp.md b/docs/advanced/figures/bench_fmg_vs_gamg_velocity_ksp.md index c68e54292..92157e4a9 100644 --- a/docs/advanced/figures/bench_fmg_vs_gamg_velocity_ksp.md +++ b/docs/advanced/figures/bench_fmg_vs_gamg_velocity_ksp.md @@ -1,3 +1,7 @@ +--- +orphan: true +--- + # Handover — FMG vs GAMG benchmark figure + adaptive-convection animation Produced 2026-06-11 from the boundary-slip / anisotropic-mover branch (PR #228). From 77f2429405496eed3494334da3443be3a89af2fd Mon Sep 17 00:00:00 2001 From: lmoresi Date: Thu, 11 Jun 2026 18:48:16 +1000 Subject: [PATCH 4/5] Address review: warn (not print) on the FMG-without-hierarchy fallback - `preconditioner='fmg'` on a mesh with no refinement hierarchy now raises a `UserWarning` via `warnings.warn` instead of printing, so it is suppressible with the standard warning filters and capturable by pytest / user code. - `test_1014` asserts the warning with `pytest.warns`. Underworld development team with AI support from Claude Code --- .../cython/petsc_generic_snes_solvers.pyx | 16 ++++++++-------- tests/test_1014_stokes_multigrid.py | 3 ++- 2 files changed, 10 insertions(+), 9 deletions(-) diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index dbb347c2a..46b35f04b 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -179,14 +179,14 @@ class SolverBaseClass(uw_object): opts.delValue(f"{prefix}{key}") self._pc_managed_value = "mg" else: - if self._preconditioner == "fmg": - if uw.mpi.rank == 0: - print( - 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.", - flush=True, - ) + 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 diff --git a/tests/test_1014_stokes_multigrid.py b/tests/test_1014_stokes_multigrid.py index 907142f1a..1a78cc8e0 100644 --- a/tests/test_1014_stokes_multigrid.py +++ b/tests/test_1014_stokes_multigrid.py @@ -79,7 +79,8 @@ 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" - stokes.solve() + with pytest.warns(UserWarning, match="falling back to GAMG"): + stokes.solve() assert _vel_pc(stokes) == "gamg" assert stokes.snes.getConvergedReason() > 0 From 5aec2b2cf591c4c8d419f7defd1e62c2f7135477 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Fri, 12 Jun 2026 12:16:37 +1000 Subject: [PATCH 5/5] Harden geometric FMG preconditioner: durable override, parallel-safe defaults, enforce Galerkin The `preconditioner` property's auto/FMG path had three problems that broke geometric multigrid on adapted (mover-deformed) meshes in parallel: 1. Override clobber on rebuild. `_apply_preconditioner_options()` runs on every `_build` so "auto" can re-resolve after a remesh. The "respect an explicit user pc_type" back-off only held on the first build: once the user's `mg` was adopted as the managed value, a later rebuild could no longer tell "user set mg" from "we set mg" and overwrote the user's tuned smoother / coarse-solver with the framework bundle. Added a durable `_pc_user_override` latch. 2. Non-parallel-safe defaults. The bundle used a bare serial `mg_coarse_pc_type=lu` (fails at np>1 with DIVERGED_LINEAR_SOLVE after 0 iterations) and `chebyshev` level smoothing (eigen-estimate-fragile on the indefinite / variable-viscosity velocity block). Switched to `redundant`+`lu` (np-safe; identical to lu in serial) and `richardson`+`sor` (the benchmark-validated mesh-independent choice). 3. Missing-Galerkin footgun. UW3's geometric MG requires Galerkin RAP coarse operators (no coarse-DM operator callbacks are installed). Selecting pc_type=mg without pc_mg_galerkin failed cryptically as PETSc error 73 in serial or DMCoarsen->DMAdaptMetric->ParMmg (3D-only) in parallel. Added `_enforce_galerkin_for_geometric_mg()` to force pc_mg_galerkin=both (with a warning) on any managed block using geometric MG, regardless of who selected it. Regression tests added to test_1014_stokes_multigrid.py (now 11/11): override survives rebuild, default bundle is parallel-safe, missing-Galerkin is repaired. Underworld development team with AI support from Claude Code --- .../cython/petsc_generic_snes_solvers.pyx | 75 +++++++++++++++++-- tests/test_1014_stokes_multigrid.py | 53 +++++++++++++ 2 files changed, 123 insertions(+), 5 deletions(-) diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index 46b35f04b..4917ede7c 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -79,6 +79,14 @@ class SolverBaseClass(uw_object): # 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): @@ -143,9 +151,15 @@ class SolverBaseClass(uw_object): # 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. + # 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: @@ -168,12 +182,22 @@ class SolverBaseClass(uw_object): # 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" - opts[f"{prefix}mg_levels_ksp_type"] = "chebyshev" + 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 - opts[f"{prefix}mg_coarse_pc_type"] = "lu" # direct coarse solve (serial) + # 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}") @@ -196,10 +220,50 @@ class SolverBaseClass(uw_object): 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_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. @@ -735,6 +799,7 @@ class SolverBaseClass(uw_object): # _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 diff --git a/tests/test_1014_stokes_multigrid.py b/tests/test_1014_stokes_multigrid.py index 1a78cc8e0..e29c2be99 100644 --- a/tests/test_1014_stokes_multigrid.py +++ b/tests/test_1014_stokes_multigrid.py @@ -101,3 +101,56 @@ def test_scalar_poisson_auto_geometric_mg(): 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