From 99761e022ab24b9e3ad43b8b0adae24302a480c5 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Wed, 29 Jul 2026 17:03:44 +1000 Subject: [PATCH 1/5] Give the geometric-MG option bundle one owner; pick up mesh-owned hierarchies on the rotated path (#468, #467) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Three routes reach a multigrid Stokes velocity block — native (PETSc interpolation between refined DMPlex levels), custom-P on the standard solve path, and custom-P through the rotated free-slip path. They are the same preconditioner reached three ways, not alternatives: custom-P is mandatory wherever native cannot go, namely rotated BCs (the DM-coupled hierarchy cannot express a per-node rotation) and adapt() children (no DMPlex refinement relation). The option bundle was written in two places and had drifted (#468). The native path had been moved to a measured gmres+sor smoother; custom_mg._configure_pcmg never received that update, and never SET mg_levels_ksp_max_it at all, so the smoother inherited whatever last wrote that options prefix — 3 left behind by the GAMG bundle on the standard path, PETSc's own PCMG default of 2 on the rotated path. The same function smoothed differently depending on what had run before it. utilities/multigrid_options.py now owns both bundles (geometric MG and the GAMG fallback) and all three writers read from it. Each bundle DERIVES the stale keys it must clear — every key a sibling bundle sets that it does not — instead of carrying a hand-maintained delete list; the derivation reproduces all three previous hand-written lists exactly. The rotated route's svd coarse solve is a named variant of the shared bundle (coarse="svd") rather than a call-site override, because the Galerkin-coarsened rotated block inherits the rigid-rotation null space where redundant/LU hits a zero pivot (#306). The #276 single-field injection fragility deliberately stays out of the bundle: it is a routing decision, not an option value. Separately (#467), mesh.adapt() leaves a coarse tail on its refinement child so every solver on an adapted mesh gets geometric MG with no per-solver call. The rotated path never consulted it — the standard path's injection hook runs after the rotated dispatch has already returned — so an adapt child under rotated free-slip fell back to GAMG, indistinguishable from having no hierarchy at all. That is the adapt-on-top-faults workflow's own configuration. custom_mg.build_transfers is now the shared "which hierarchy does this solver get?" rule for both paths, with the same opportunistic barycentric->RBF fallback and degrade-to-GAMG on failure. Measured (same operator, RHS and coarse solve; two-level hierarchy, i.e. where the native measurement says the gmres margin is SMALLEST): route pre-fix unified custom-P standard (isotropic annulus) 5 vel its 4 vel its custom-P rotated (TI, eta1/eta0 = 1e-3) 11 vel its 5 vel its 0.683 s 0.392 s (1.74x, linear solve timed in isolation) Also here: the rotated FMG branch now drops its option keys after setup, as the KSP's own keys already were. The per-solve prefix is unique, so those keys previously accumulated in the global database once per timestep. A new test covers the failure mode that introduces — a later setFromOptions finding pc_type gone and abandoning the multigrid — across eight Newton increments. Tests: tests/test_1021_mg_option_bundle.py reads the smoother configuration back off the LIVE PCs for all three routes and asserts they agree, with the coarse difference asserted rather than tolerated (an options-database assertion would not have caught the drift, which was a key nobody wrote). Verified to fail when the pre-fix bundle is reimposed. The mesh-owned pickup is covered serially there and at np>1 in tests/parallel/test_1064, where the transfers are built cross-partition. Validation: test_1014/1015/1016/1017/1018/1020, test_0753/0835/0836/0840 pass; full "level_1 and tier_a" sweep 526 passed with 2 pre-existing failures (#470). Parallel np2 and np4: test_1017_custom_mg_parallel 6 passed, test_1064_rotated_freeslip_parallel 9 passed, test_1066_rotated_datum_parallel 2 passed. No goldens moved — GOLDEN_ANNULUS_FMG is a converged-solution quantity, so a better preconditioner reaches the same answer, and the iteration guards are upper bounds a better smoother only helps. Underworld development team with AI support from Claude Code --- docs/developer/CHANGELOG.md | 43 ++++ docs/developer/subsystems/rotated-freeslip.md | 33 ++- docs/developer/subsystems/solvers.md | 79 +++++- .../cython/petsc_generic_snes_solvers.pyx | 61 +---- src/underworld3/utilities/custom_mg.py | 117 ++++++--- .../utilities/multigrid_options.py | 198 +++++++++++++++ src/underworld3/utilities/rotated_bc.py | 85 ++++--- .../test_1064_rotated_freeslip_parallel.py | 38 ++- tests/test_1021_mg_option_bundle.py | 240 ++++++++++++++++++ 9 files changed, 770 insertions(+), 124 deletions(-) create mode 100644 src/underworld3/utilities/multigrid_options.py create mode 100644 tests/test_1021_mg_option_bundle.py diff --git a/docs/developer/CHANGELOG.md b/docs/developer/CHANGELOG.md index 2a8f28f2..036d63ce 100644 --- a/docs/developer/CHANGELOG.md +++ b/docs/developer/CHANGELOG.md @@ -6,6 +6,49 @@ This log tracks significant development work at a conceptual level, suitable for ## 2026 Q3 (July – September) +### One Owner for the Geometric-Multigrid Option Bundle (July 2026) + +**The PETSc option bundle that configures a Stokes velocity block's multigrid +now lives in exactly one module, and all three routes that reach that block read +it from there** (#468), **and rotated free-slip now picks up a mesh-owned +multigrid hierarchy instead of silently discarding it** (#467). + +Three routes reach a multigrid velocity block: native (PETSc interpolation +between refined DMPlex levels), custom-P on the standard solve path, and custom-P +through the rotated free-slip path. They are the same preconditioner reached +three ways, not alternatives — custom-P is *mandatory* wherever native cannot go, +namely rotated boundary conditions and `adapt()` children. The bundle was +written in two places and had drifted: the native path had been moved to a +`gmres`+`sor` smoother on a recorded measurement, and the custom-P routes had +not. Worse, the custom-P writer never *set* the smoother iteration count at all, +so it inherited whatever had last written that options prefix — 3 left behind by +the GAMG bundle on the standard path, PETSc's own default of 2 on the rotated +path. The same function smoothed differently depending on what had run before +it. + +Unifying the bundle recovers, on the same operator, right-hand side and coarse +solve: rotated custom-P velocity-block iterations 11 → 5 (0.68 s → 0.39 s of +linear solve, timed in isolation), standard custom-P 5 → 4, on a *two-level* +hierarchy — the depth at which the native measurement says the gmres margin is +smallest. The bundle also now derives which stale keys it must clear rather than +carrying a hand-maintained list, which is what let the iteration count go unset +in the first place. + +Separately, `mesh.adapt()` leaves a coarse tail on its refinement child so that +every solver on an adapted mesh gets geometric multigrid with no per-solver call. +The rotated path never consulted it — the standard path's injection hook runs +after the rotated dispatch has already returned — so an adapt child under rotated +free-slip fell back to algebraic multigrid, indistinguishable from having no +hierarchy at all. That is the `adapt-on-top-faults` workflow's own configuration +(a fault resolved by local refinement, with rotated free-slip chosen because it +composes with transverse isotropy). Both paths now resolve the hierarchy through +one shared rule, with the same opportunistic degrade-to-GAMG behaviour. + +The regression test reads the smoother configuration back off the **live PETSc +objects** for all three routes and asserts they agree. An options-database +assertion would not have caught the original drift, because the drift was +precisely a key nobody wrote. + ### One Rotated Free-Slip Path, Now With a Prescribed Wall-Normal Velocity (July 2026) **Rotated strong free-slip now takes a prescribed wall-normal velocity datum diff --git a/docs/developer/subsystems/rotated-freeslip.md b/docs/developer/subsystems/rotated-freeslip.md index 14737fe1..7c47271e 100644 --- a/docs/developer/subsystems/rotated-freeslip.md +++ b/docs/developer/subsystems/rotated-freeslip.md @@ -41,11 +41,42 @@ There is a single driver, `solve_rotated_freeslip`: a manual outer Newton/Picard loop that rotates the residual and tangent every iteration (`F̂ = Q F`, `Ĵ = Q J Qᵀ`), imposes the constraint on the rotated normal rows, and solves each increment with a self-contained fieldsplit-Schur KSP (geometric -FMG on the custom prolongation when a hierarchy is registered, else GAMG; the +FMG on the custom prolongation when the solver has a hierarchy, else GAMG; the native 1/μ pressure mass as the Schur preconditioner). There is no separate linear path and no up-front nonlinearity probe: a linear model converges after its first increment and the loop self-terminates. +### Where the multigrid hierarchy comes from + +Native geometric FMG cannot serve a rotated solve at all — the DM-coupled +hierarchy has no way to express the per-node rotation — so custom-P multigrid is +not an optimisation here, it is the only geometric option. The hierarchy is +resolved by `custom_mg.build_transfers`, which is the **same rule the standard +path uses**: + +1. an explicit `set_custom_fmg` registration on the solver wins; otherwise +2. a **mesh-owned** coarse tail (`mesh._custom_mg_coarse_meshes`, what + `mesh.adapt()` leaves on a refinement child) is picked up opportunistically, + with the barycentric→RBF builder fallback, so a failed build degrades to GAMG + rather than crashing the solve. + +Step 2 used to be unreachable from here. The standard path's injection hook runs +*after* the rotated dispatch has already returned, and the rotated builder +consulted only `solver._custom_mg` — so an `adapt()` child under rotated +free-slip reported `velocity_pc == "GAMG"`, indistinguishable from having no +hierarchy at all, and lost its multigrid silently (#467). That is the +`adapt-on-top-faults` workflow's own configuration: a fault resolved by local +refinement on an adapt child, with rotated free-slip chosen over Nitsche because +it composes with transverse isotropy. + +The option bundle for both the FMG and the GAMG velocity block comes from +`utilities/multigrid_options.py`, shared with the native and standard custom-P +routes, so the rotated block cannot be configured differently by accident. The +single deliberate difference is the coarse solve: `coarse="svd"`, because the +Galerkin-coarsened rotated block inherits the rigid-rotation null space where +`redundant`/LU hits a zero pivot. See +[the solvers subsystem doc](solvers.md) for the bundle itself. + ### The velocity sub-KSP must converge, not just apply The velocity block is preconditioned by multigrid, but the KSP *around* that diff --git a/docs/developer/subsystems/solvers.md b/docs/developer/subsystems/solvers.md index 32eaff42..3f492df1 100644 --- a/docs/developer/subsystems/solvers.md +++ b/docs/developer/subsystems/solvers.md @@ -175,7 +175,7 @@ def validate_unknowns_sharing(multi_material_model): assert model.Unknowns.DFDt is reference_unknowns.DFDt, \ f"Model {i} $D\mathbf{{F}}/Dt$ not shared - stress history will be wrong" ``` -- ⚠️ Preconditioner *selection* is partly covered — see "Choosing the Krylov method for a fieldsplit sub-solve" below; multigrid and Schur preconditioner choice is still undocumented +- ⚠️ Preconditioner *selection* is partly covered — see "Choosing the Krylov method for a fieldsplit sub-solve" and "The multigrid option bundle, and its one owner" below; Schur preconditioner choice is still undocumented - Could benefit from optimization examples ## Choosing the Krylov method for a fieldsplit sub-solve @@ -327,6 +327,81 @@ application, and the velocity KSP is applied once per Schur `MatMult` — i.e. o pressure Krylov iteration. That number is a sample, not work. `SolverInstrumentation` (`systems/solver_health.py`) is the mechanism that sums properly. +## The multigrid option bundle, and its one owner + +Three routes reach a multigrid velocity block, and they are **the same +preconditioner reached three ways**, not three alternatives: + +| route | when | prolongation | +|---|---|---| +| native | mesh built with `refinement >= 1`, ordinary BCs | PETSc `DMCreateInterpolation` between refined DMPlex levels | +| custom-P, standard path | `set_custom_fmg`, or an `adapt()` child's mesh-owned coarse tail | barycentric / RBF, Galerkin coarse operators | +| custom-P, rotated path | rotated free-slip, via `rotated_bc` | as above, with the fine prolongation rotated, `P̂ = Q_v·P` | + +custom-P is **mandatory** wherever native cannot go: rotated boundary conditions +(the DM-coupled hierarchy cannot express a per-node rotation) and non-nested +grids (`adapt()` children have no DMPlex refinement relation). So the routes are +not ranked — the one that matters most for adapted and curved-boundary work is +custom-P. + +The option *values* live in one module, `utilities/multigrid_options.py`. Every +writer reads a bundle from there and applies it to its own options object under +its own prefix; nobody writes a multigrid option value anywhere else. That is +structural rather than stylistic: when the bundle was written in two places, the +native path was deliberately moved to a measured `gmres`+`sor` smoother and the +custom-P path was not, and the custom-P routes ran `richardson` at an iteration +count **nobody had set** — inherited from whatever last wrote that options +prefix (3 left behind by the GAMG bundle on the standard path, PETSc's own PCMG +default of 2 on the rotated path). The same function smoothed differently +depending on what had run before it. + +### What a bundle carries + +A bundle is the settings it sets *and* the keys it must clear. The clear-list is +derived, not hand-written: it is every key any sibling bundle sets that this one +does not. These bundles share an options prefix, so switching a block from GAMG +to geometric MG leaves the GAMG-only keys behind and `setFromOptions` will +happily re-read them. Deriving the list means a key added to one bundle +automatically becomes stale for the others. + +Two consequences worth stating outright: + +- **Set the smoother iteration count, never inherit it.** A bundle that omits + `mg_levels_ksp_max_it` is not "using the default" — it is using whatever the + last writer left. +- **The measured smoother is `gmres`+`sor` at 4 iterations.** Chebyshev needs + eigenvalue estimates of the smoothed operator, which are fragile on the + indefinite / variable-viscosity velocity block. Richardson is stationary and + degrades on the non-symmetric operator the consistent-Newton tangent produces: + measured on the Spiegelman notch (Drucker–Prager, η contrast 1e26, nested + 4-level hierarchy) the per-cycle contraction is ρ = 0.75 richardson against + 0.56 gmres at the *same* four iterations, and the gmres margin **grows with + depth** (5% at 3 levels, 25% at 4). + +### The one legitimate per-route difference + +The coarse solve. The Galerkin-coarsened **rotated** velocity block inherits the +rigid-rotation null space of the constrained problem (a closed circle: one mode; +a spherical shell: three), and `redundant`/LU hits a zero pivot there — +`SUBPC_ERROR`, outer reason −11. So the rotated route asks for +`geometric_mg_bundle(coarse="svd")`, which is null-space robust and cheap on a +small coarse level. This is a named variant of the shared bundle, not a +call-site override, so it is visible in the same place as everything else. + +The other native/custom-P asymmetry — that native FMG is unusable for +single-field solvers because `DMCreateInjection` cannot reliably be built on a +refined DMPlex (#276) — is deliberately *not* in the bundle. It is a routing +decision (which route a solver may take), not an option value, and lives with +the route choice in `_apply_preconditioner_options`. + +### Testing it + +`tests/test_1021_mg_option_bundle.py` reads the smoother configuration back off +the **live PETSc objects** after setup — not out of the options database — for +all three routes and asserts they agree, with the coarse-solve difference +asserted rather than tolerated. Options-database assertions would not have caught +the original drift, because the drift was precisely a key that was never written. + ## Critical Stability Note ```{warning} Solver Stability is Paramount @@ -337,7 +412,7 @@ pressure Krylov iteration. That number is a sample, not work. `SolverInstrumenta ```{note} For Contributors This well-documented subsystem could benefit from: -- Preconditioner selection guidance (Krylov sub-solve choice is now covered; multigrid and Schur preconditioner selection is not) +- Preconditioner selection guidance (Krylov sub-solve choice and the multigrid option bundle are now covered; Schur preconditioner selection is not) - Performance tuning documentation - Convergence analysis examples - Scaling studies and optimization diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index 5ff3ce6b..21c3520b 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -12,6 +12,7 @@ from underworld3.utilities._jitextension import getext, JITCallbackSet import underworld3.timing as timing from underworld3.utilities._api_tools import uw_object +from underworld3.utilities import multigrid_options from underworld3.function import expression as public_expression expression = lambda *x, **X: public_expression(*x, _unique_name_generation=True, **X) @@ -618,46 +619,14 @@ class SolverBaseClass(uw_object): ) want_fmg = False + # The option VALUES live in utilities.multigrid_options, which is the + # single owner shared with the custom-P routes (custom_mg, rotated_bc) — + # the routes are the same preconditioner reached three ways and must not + # be configured from three places (#468). Coarse solve: the native + # hierarchy is not rotated, so its coarse operator carries no inherited + # null space and redundant+LU is right here. 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 - # gmres+sor, sized for a DEEP hierarchy (the only kind worth having: - # a two-level cycle is a coarse-grid correction, not a V-cycle, and is - # not worth special-casing). Chebyshev needs eigenvalue estimates of the - # smoothed operator, which are fragile on the indefinite / - # variable-viscosity velocity block and diverge. Richardson is - # stationary and degrades on the NON-SYMMETRIC operator produced by the - # consistent-Newton tangent. Measured on the Spiegelman notch (Drucker- - # Prager, eta contrast 1e26) over a nested 4-level hierarchy: contraction - # per V-cycle rho = 0.75 (richardson) vs 0.56 (gmres) at the SAME four - # smoother iterations -- and the gmres margin GROWS with depth (5% at 3 - # levels, 25% at 4), because deeper cycles apply the smoother on more - # coarse operators. Four iterations, not more: per unit work gmres/4 - # (rho^(1/4) = 0.87) beats gmres/8 (0.91). - opts[f"{prefix}mg_levels_ksp_type"] = "gmres" - opts[f"{prefix}mg_levels_pc_type"] = "sor" - opts[f"{prefix}mg_levels_ksp_max_it"] = 4 - # Run EXACTLY max_it smoother iterations: no residual-norm computation - # and no convergence test, so every V-cycle costs the same. A Krylov - # smoother makes the cycle non-stationary, which is why the velocity - # block is fgmres (flexible) rather than gmres -- see the fieldsplit - # defaults in the Stokes __init__. - opts[f"{prefix}mg_levels_ksp_norm_type"] = "none" - 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}") + multigrid_options.geometric_mg_bundle().apply(opts, prefix) self._pc_managed_value = "mg" else: if self._preconditioner == "fmg" and n_levels <= 1 and uw.mpi.rank == 0: @@ -668,19 +637,7 @@ class SolverBaseClass(uw_object): 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_levels_ksp_norm_type", - "mg_coarse_pc_type", - "mg_coarse_redundant_pc_type"): - opts.delValue(f"{prefix}{key}") + multigrid_options.gamg_bundle().apply(opts, prefix) self._pc_managed_value = "gamg" def _enforce_galerkin_for_geometric_mg(self): diff --git a/src/underworld3/utilities/custom_mg.py b/src/underworld3/utilities/custom_mg.py index d4de0371..5d96d272 100644 --- a/src/underworld3/utilities/custom_mg.py +++ b/src/underworld3/utilities/custom_mg.py @@ -43,6 +43,8 @@ import numpy as np from petsc4py import PETSc +from underworld3.utilities import multigrid_options + __all__ = ["barycentric_prolongation", "rbf_prolongation", "inject_custom_mg", "CustomMGHierarchy", "set_custom_fmg", "sbr_refine", "sbr_refine_where", "nvb_refine"] @@ -593,30 +595,29 @@ def _assert_no_zero_columns_serial(P_csr, level): f"operator would be singular.") -def _configure_pcmg(pc, Ps): +def _configure_pcmg(pc, Ps, coarse="redundant"): """Reconfigure ``pc`` as a fresh PCMG (FMG F-cycle) driven by the supplied reduced->reduced prolongations ``Ps``, Galerkin RAP for coarse operators. - Writes the MG bundle into the options DB under the PC's OWN options prefix + The option VALUES come from :func:`multigrid_options.geometric_mg_bundle` — + the same bundle the native (DMPlex-refinement) route applies — so custom-P + and native multigrid cannot be configured differently (#468). ``coarse`` + selects the coarse-solve variant: ``"svd"`` on the rotated path, whose + Galerkin-coarsened operator inherits the rigid-rotation null space. + + The bundle is written into the options DB under the PC's OWN options prefix (``pc.getOptionsPrefix()`` — e.g. ``Solver_N_`` for the scalar top-level PC, - ``Solver_N_fieldsplit_velocity_`` for the Stokes velocity sub-PC) and removes - any gamg keys BEFORE ``setFromOptions`` — otherwise ``setFromOptions`` re-reads - a lingering ``pc_type=gamg`` and reverts ``setType("mg")``. ``setMGInterpolation`` - persists through ``setFromOptions``; the first ``PCSetUp`` builds the coarse - operators from our P (no ``MatProductReplaceMats`` shape bug, since the PCMG is - fresh and P's size is fixed).""" + ``Solver_N_fieldsplit_velocity_`` for the Stokes velocity sub-PC) BEFORE + ``setFromOptions``, and it clears the gamg keys — otherwise ``setFromOptions`` + re-reads a lingering ``pc_type=gamg`` and reverts ``setType("mg")``. + ``setMGInterpolation`` persists through ``setFromOptions``; the first + ``PCSetUp`` builds the coarse operators from our P (no + ``MatProductReplaceMats`` shape bug, since the PCMG is fresh and P's size is + fixed).""" nlev = len(Ps) + 1 prefix = pc.getOptionsPrefix() or "" - opts = PETSc.Options() - opts.setValue(prefix + "pc_type", "mg") - opts.setValue(prefix + "pc_mg_type", "full") - opts.setValue(prefix + "pc_mg_galerkin", "both") - opts.setValue(prefix + "mg_levels_ksp_type", "richardson") - opts.setValue(prefix + "mg_levels_pc_type", "sor") - opts.setValue(prefix + "mg_coarse_pc_type", "redundant") - opts.setValue(prefix + "mg_coarse_redundant_pc_type", "lu") - for key in ("pc_gamg_type", "pc_gamg_repartition", "pc_gamg_agg_nsmooths"): - opts.delValue(prefix + key) + multigrid_options.geometric_mg_bundle(coarse=coarse).apply( + PETSc.Options(), prefix) pc.setType("mg") pc.setMGLevels(nlev) pc.setMGType(PETSc.PC.MGType.FULL) @@ -963,25 +964,51 @@ def set_custom_fmg(solver, coarse_meshes, *, builder="barycentric", solver.is_setup = False -def auto_inject_custom_mg(solver, field_id=None): - """Solve-hook entry: inject custom-P FMG from either a solver-set hierarchy - (``set_custom_fmg``) or a **mesh-owned** one (``mesh.adapt`` refinement child). +def build_transfers(solver, field_id=None): + """The custom-P prolongations this solver should drive, built and ready to + install — from either a solver-set hierarchy (``set_custom_fmg``) or a + **mesh-owned** one (a ``mesh.adapt`` refinement child). + + This is the shared resolution rule for every route that can drive custom-P + multigrid: the standard solve path via :func:`auto_inject_custom_mg`, and the + rotated free-slip path, which builds its own IS-based fieldsplit inside + ``utilities.rotated_bc`` and so never reaches the standard injection hook + (#467). Both must answer "which hierarchy does this solver get?" the same way. A refinement child carries ``mesh._custom_mg_coarse_meshes`` (the static - coarse tail). The first time a solver on such a mesh solves, we lazily build a - :class:`CustomMGHierarchy` ``[*coarse, solver.mesh]`` targeting ``field_id`` - (0 for the Stokes velocity block, None for scalar/vector) and register it on - the solver — so every solver on an adapted mesh drives geometric MG with no - per-solver call. A solver-set hierarchy (if present) always wins. + coarse tail), so a :class:`CustomMGHierarchy` ``[*coarse, solver.mesh]`` + targeting ``field_id`` (0 for the Stokes velocity block, None for + scalar/vector) is built lazily on first solve — every solver on an adapted + mesh drives geometric MG with no per-solver call. A solver-set hierarchy (if + present) always wins. + + Parameters + ---------- + solver : SolverBaseClass + Built solver (``_build`` has run) whose hierarchy is wanted. + field_id : int or None + Field index for multi-field solvers (0 = Stokes velocity), None = single + field. + + Returns + ------- + (CustomMGHierarchy, list of Mat) or (None, None) + ``(None, None)`` means "no hierarchy here" — either none is registered or + an OPPORTUNISTIC mesh-owned build failed and the caller should keep its + default preconditioner. A build failure on a solver-set hierarchy raises + instead: the user asked for it explicitly. """ - # Solver-set hierarchy (set_custom_fmg): the user asked for it explicitly — - # build + install directly and let any error surface. - if solver._custom_mg is not None: - inject_custom_mg(solver) - return + cfg = getattr(solver, "_custom_mg", None) + if cfg is not None: + if not (isinstance(cfg, dict) and cfg.get("mode") == "hierarchy"): + raise NotImplementedError( + "the legacy set_custom_mg registration has no hierarchy to " + "resolve; use set_custom_fmg().") + h = cfg["hierarchy"] + return h, h.build(solver) # explicit request: errors surface # Mesh-owned hierarchy (adapt() child): OPPORTUNISTIC auto-pickup. It must never - # crash a solve, so build the transfers (which now validate the finest reduced + # crash a solve, so build the transfers (which validate the finest reduced # map against the assembled operator — see CustomMGHierarchy.build) inside a # try/except and fall back to the solver's default preconditioner on any failure. # The finest map is derived from the DM section AFTER snes.setUp() finalizes it, @@ -989,7 +1016,7 @@ def auto_inject_custom_mg(solver, field_id=None): # semi-Lagrangian advection-diffusion (which earlier had to be skipped). coarse = getattr(solver.mesh, "_custom_mg_coarse_meshes", None) if coarse is None: - return # nothing to inject + return None, None # nothing to inject builder = getattr(solver.mesh, "_custom_mg_builder", "barycentric") # Retry with the RBF builder before abandoning geometric MG. The @@ -1029,7 +1056,29 @@ def auto_inject_custom_mg(solver, field_id=None): warnings.warn( f"custom_mg: mesh-owned FMG build failed ({exc}); using the " "solver's default preconditioner.") - return + return None, None + + return h, Ps + + +def auto_inject_custom_mg(solver, field_id=None): + """Solve-hook entry on the STANDARD path: resolve this solver's custom-P + hierarchy (:func:`build_transfers`) and install it on the managed PC. + + The rotated free-slip path does not come through here — it builds its own + IS-based fieldsplit and calls :func:`build_transfers` directly. + """ + # Solver-set hierarchy (set_custom_fmg, or the deprecated set_custom_mg): + # the user asked for it explicitly — build + install directly and let any + # error surface. inject_custom_mg is the one place that still understands the + # legacy registration. + if solver._custom_mg is not None: + inject_custom_mg(solver) + return + + h, Ps = build_transfers(solver, field_id=field_id) + if h is None: + return # Dimensional guard (checkable for the monolithic operator, field_id is None): # the finest transfer must chain to the operator PCMG will Galerkin against. diff --git a/src/underworld3/utilities/multigrid_options.py b/src/underworld3/utilities/multigrid_options.py new file mode 100644 index 00000000..2245a371 --- /dev/null +++ b/src/underworld3/utilities/multigrid_options.py @@ -0,0 +1,198 @@ +r""" +The multigrid option bundles for a solver's managed preconditioner block. + +Three routes reach a multigrid velocity block in UW3, and all three must be +configured the same way: + +============================ ============================================== +route prolongation +============================ ============================================== +native PETSc ``DMCreateInterpolation`` between refined + DMPlex levels (mesh built with ``refinement>=1``) +custom-P, standard path barycentric / RBF, Galerkin coarse operators + (:mod:`underworld3.utilities.custom_mg`) +custom-P, rotated path as above, with the FINE prolongation rotated + (:mod:`underworld3.utilities.rotated_bc`) +============================ ============================================== + +custom-P is mandatory wherever native cannot go — rotated boundary conditions +(the DM-coupled hierarchy cannot express the per-node rotation) and non-nested +grids (``adapt()`` children have no DMPlex refinement relation) — so the routes +are not alternatives, they are the same preconditioner reached three ways. + +This module is where the bundles live so that they cannot drift apart. Each +writer reads a bundle from here and applies it to its own options object under +its own prefix; nobody writes a multigrid option value anywhere else. + +A bundle carries two things: the settings it *sets*, and the **stale** keys it +must *remove*. The stale list is derived, not hand-written: it is every key any +other bundle sets that this one does not. That matters because these bundles are +written into a shared options database under a shared prefix — toggling a block +from GAMG to geometric MG leaves the GAMG-only keys behind, and ``setFromOptions`` +will happily re-read them. Hand-maintained delete lists are exactly what let the +smoother iteration count go unset and be inherited from whatever ran before +(issue #468). + +Examples +-------- +Apply the geometric bundle to a solver's managed block:: + + from underworld3.utilities import multigrid_options + multigrid_options.geometric_mg_bundle().apply(self.petsc_options, prefix) + +Read the settings without writing them (the rotated path stages its options in a +dict so it can drop them again after ``setUp``):: + + cfg.update({f"fieldsplit_vel_{k}": v + for k, v in multigrid_options.gamg_bundle().settings.items()}) +""" + +from typing import NamedTuple + +__all__ = ["MGBundle", "geometric_mg_bundle", "gamg_bundle", + "GEOMETRIC_MG_COARSE_SOLVERS"] + + +#: Coarse-solve variants of the geometric bundle. ``"redundant"`` (redundant+LU) +#: is the default and is np-safe; ``"svd"`` is required whenever the coarse +#: operator inherits a null space — see :func:`geometric_mg_bundle`. +GEOMETRIC_MG_COARSE_SOLVERS = ("redundant", "svd") + + +class MGBundle(NamedTuple): + """A preconditioner option bundle: the keys to set, and the keys to clear. + + ``settings`` maps an option suffix (no prefix, no leading ``-``) to its + value; ``None`` means a bare PETSc flag. ``stale`` lists suffixes this + bundle does not own but a sibling bundle does, which must be removed before + ``setFromOptions`` so a previously-applied bundle cannot leak through. + """ + + settings: dict + stale: tuple + + def apply(self, opts, prefix=""): + """Write this bundle into ``opts`` (a ``PETSc.Options``) under ``prefix``. + + ``prefix`` is the managed block's own prefix — ``""`` for a top-level PC, + ``"fieldsplit_velocity_"`` for the Stokes velocity sub-block — and is + applied on top of whatever prefix ``opts`` itself carries. + """ + for key, value in self.settings.items(): + opts.setValue(prefix + key, value) + for key in self.stale: + opts.delValue(prefix + key) + + +def _geometric_mg_settings(coarse): + """The geometric-MG settings for one coarse-solve variant.""" + settings = { + "pc_type": "mg", + "pc_mg_type": "full", # FMG (F-cycle) + # Galerkin (RAP) coarse operators are REQUIRED: UW3 installs no + # residual/Jacobian callbacks on the coarse DMs, so PETSc cannot + # re-discretise the operator there. + "pc_mg_galerkin": "both", + # gmres+sor, sized for a DEEP hierarchy (the only kind worth having: a + # two-level cycle is a coarse-grid correction, not a V-cycle, and is not + # worth special-casing). Chebyshev needs eigenvalue estimates of the + # smoothed operator, which are fragile on the indefinite / + # variable-viscosity velocity block and diverge. Richardson is stationary + # and degrades on the NON-SYMMETRIC operator produced by the + # consistent-Newton tangent. Measured on the Spiegelman notch + # (Drucker-Prager, eta contrast 1e26) over a nested 4-level hierarchy: + # contraction per V-cycle rho = 0.75 (richardson) vs 0.56 (gmres) at the + # SAME four smoother iterations -- and the gmres margin GROWS with depth + # (5% at 3 levels, 25% at 4), because deeper cycles apply the smoother on + # more coarse operators. Four iterations, not more: per unit work gmres/4 + # (rho^(1/4) = 0.87) beats gmres/8 (0.91). + "mg_levels_ksp_type": "gmres", + "mg_levels_pc_type": "sor", + # SET the count, never inherit it. PCMG's own default is 2 and the GAMG + # bundle below leaves 3 under the same prefix, so a bundle that omits + # this key smooths differently depending on what ran before it (#468). + "mg_levels_ksp_max_it": 4, + # Run EXACTLY max_it smoother iterations: no residual-norm computation and + # no convergence test, so every V-cycle costs the same. A Krylov smoother + # makes the cycle non-stationary, which is why the velocity block is + # fgmres (flexible) rather than gmres. + "mg_levels_ksp_norm_type": "none", + "mg_levels_ksp_converged_maxits": None, + } + if coarse == "redundant": + # 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 — np-safe by default without surprising + # small-np users. + settings["mg_coarse_pc_type"] = "redundant" + settings["mg_coarse_redundant_pc_type"] = "lu" + elif coarse == "svd": + # The Galerkin-coarsened ROTATED velocity block inherits every + # rigid-rotation null-space mode of the constrained problem (a closed + # circle: one; a spherical shell: three), and redundant/LU hits a zero + # pivot there (SUBPC_ERROR, outer reason -11 — the #306 fix). SVD is + # null-space robust and the coarse level is small. Same choice as the + # native spherical FMG setups. + settings["mg_coarse_pc_type"] = "svd" + else: + raise ValueError( + f"coarse must be one of {GEOMETRIC_MG_COARSE_SOLVERS} (got {coarse!r})") + return settings + + +def _gamg_settings(): + """The algebraic-multigrid (GAMG) settings — the fallback whenever no + geometric hierarchy is available.""" + return { + "pc_type": "gamg", + "pc_gamg_type": "agg", + "pc_gamg_repartition": True, + "pc_mg_type": "additive", + "pc_gamg_agg_nsmooths": 2, + "mg_levels_ksp_max_it": 3, + "mg_levels_ksp_converged_maxits": None, + } + + +def _all_keys(): + """Every option key any bundle here owns — the basis for the stale lists.""" + keys = set(_gamg_settings()) + for coarse in GEOMETRIC_MG_COARSE_SOLVERS: + keys |= set(_geometric_mg_settings(coarse)) + return keys + + +def _bundle(settings): + return MGBundle(settings=settings, + stale=tuple(sorted(_all_keys() - set(settings)))) + + +def geometric_mg_bundle(coarse="redundant"): + """Geometric multigrid (FMG F-cycle) on an explicit hierarchy. + + Parameters + ---------- + coarse : {"redundant", "svd"} + The coarse-level solve. ``"redundant"`` (redundant+LU) is the default. + Use ``"svd"`` when the coarse operator inherits a null space — which the + Galerkin-coarsened *rotated* velocity block always does, because the + rigid rotations survive the constraint. + + Returns + ------- + MGBundle + Settings and stale keys; see :meth:`MGBundle.apply`. + """ + return _bundle(_geometric_mg_settings(coarse)) + + +def gamg_bundle(): + """Algebraic multigrid — the fallback when no geometric hierarchy exists. + + Returns + ------- + MGBundle + Settings and stale keys; see :meth:`MGBundle.apply`. + """ + return _bundle(_gamg_settings()) diff --git a/src/underworld3/utilities/rotated_bc.py b/src/underworld3/utilities/rotated_bc.py index f755374d..b83bad4d 100644 --- a/src/underworld3/utilities/rotated_bc.py +++ b/src/underworld3/utilities/rotated_bc.py @@ -9,9 +9,12 @@ no separate linear path and no up-front nonlinearity probe. Each increment is solved by a self-contained fieldsplit-Schur KSP by default: the -velocity block is geometric FMG on the custom prolongation (``set_custom_fmg``) when a -hierarchy is registered (the PREFERRED route), else GAMG tuned to the native path's -settings; the Schur complement is preconditioned by the native 1/mu pressure mass +velocity block is geometric FMG on the custom prolongation whenever this solver has a +hierarchy — registered by ``set_custom_fmg`` or owned by an ``adapt()`` child mesh — +which is the PREFERRED route, else GAMG; the option bundle for both is the shared one +in ``underworld3.utilities.multigrid_options``, so the rotated velocity block is +configured identically to the native and standard custom-P routes. +The Schur complement is preconditioned by the native 1/mu pressure mass (the Pmat p-p block, exactly the standard path's ``schur_precondition=a11``), with a constant-pressure nullspace on the inner Schur solve for enclosed domains; direct MUMPS LU is opt-in via ``solver._rotated_use_lu``. @@ -32,6 +35,11 @@ from underworld3.utilities.boundary_flux import ( _boundary_stratum_is, _desmear, write_boundary_scalar_field) +# The multigrid option bundles are owned in ONE place and shared with the native +# and standard custom-P routes, so the rotated velocity block cannot be +# configured differently from the others (#468). +from underworld3.utilities import multigrid_options + # Monotonic counter so each rotated solve gets a UNIQUE PETSc options prefix. With a # fixed prefix, sequential rotated solves (e.g. two solvers in one script, or one # solver in a time-stepping loop) share and re-set the same global-options keys; the @@ -847,19 +855,33 @@ def rotated_residual(uvec, keep_cartesian=False): def _build_rotated_custom_Pl(solver, Q, normal_rows): """The rotated custom-FMG prolongation list [*coarse, Q_v·P_fine] for the - velocity block, or None if no hierarchy is registered. Depends only on Q and + velocity block, or None if this solver has no hierarchy. Depends only on Q and the mesh (NOT the solution), so the rotated Newton loop builds it ONCE and reuses - it across Newton iterations (the prolongation build is the expensive part).""" - if getattr(solver, "_custom_mg", None) is None: + it across Newton iterations (the prolongation build is the expensive part). + + The hierarchy is resolved by ``custom_mg.build_transfers``, which is the same + rule the standard path uses: an explicit ``set_custom_fmg`` registration wins, + otherwise a MESH-OWNED coarse tail (what ``mesh.adapt()`` leaves on a + refinement child) is picked up opportunistically. That mesh-owned pickup used + to be unreachable from here — the standard path's injection hook runs after + the rotated dispatch has already returned — so an adapt child under rotated + free-slip silently lost its multigrid and solved on GAMG (#467). This is the + ``adapt-on-top-faults`` workflow's own configuration.""" + from underworld3.utilities import custom_mg + h, Ps = custom_mg.build_transfers(solver, field_id=0) + if h is None: return None vel_is = solver._subdict["velocity"][0] vis = np.asarray(vel_is.getIndices()) g2blk = {int(g): k for k, g in enumerate(vis)} Qv = Q.createSubMatrix(vel_is, vel_is) nrows_blk = sorted({g2blk[g] for g in normal_rows if g in g2blk}) - Ps = solver._custom_mg["hierarchy"].build(solver) Pfine = Qv.matMult(Ps[-1]) Pfine.zeroRows(nrows_blk, diag=0.0) + # Remember an opportunistic mesh-owned pickup, so a later solve on this solver + # (rotated or not) reuses the hierarchy instead of re-resolving it. + if getattr(solver, "_custom_mg", None) is None: + solver._custom_mg = {"mode": "hierarchy", "hierarchy": h, "verbose": False} return list(Ps[:-1]) + [Pfine] @@ -909,8 +931,9 @@ def _solve_rotated_iterative(solver, Ahat, bhat, Q, Qt, normal_rows, verbose=Fal custom_Pl=None, nsp=None, Mp=None, ctx=None): """Solve the rotated saddle with a SELF-CONTAINED fieldsplit-Schur KSP on the rotated operator. The velocity block is geometric FMG on the CUSTOM prolongation - (PR#290, rotated) when a hierarchy is registered (``set_custom_fmg``), else GAMG - (tuned to the native path's settings). + (PR#290, rotated) whenever the solver has a hierarchy — ``set_custom_fmg`` or a + mesh-owned ``adapt()`` tail — else GAMG. Both bundles come from + ``utilities.multigrid_options``, shared with the native and standard routes. A plain rotated Mat has no DM field info, so UW3's DM-coupled fieldsplit cannot split it — we build the split from EXPLICIT velocity/pressure index sets. For the @@ -988,21 +1011,18 @@ def _solve_rotated_iterative(solver, Ahat, bhat, Q, Qt, normal_rows, verbose=Fal cfg["pc_fieldsplit_schur_precondition"] = "selfp" cfg["fieldsplit_pres_pc_type"] = "jacobi" if custom_Pl is None: - # GAMG fallback velocity block, tuned to native parity (pyx Stokes - # defaults). NOTE: the custom-FMG route is the preferred velocity - # block — this applies only when no hierarchy is registered. + # GAMG fallback velocity block. The bundle is the SHARED one (see + # utilities.multigrid_options), so it cannot drift from what the + # standard path applies. NOTE: the custom-FMG route is the preferred + # velocity block — this applies only when no hierarchy is available. + # No stale keys to clear: `pfx` is unique per rotated solve. cfg.update({ "fieldsplit_vel_ksp_type": "fgmres", "fieldsplit_vel_ksp_rtol": str(tol * 0.033), "fieldsplit_vel_ksp_max_it": "200", - "fieldsplit_vel_pc_type": "gamg", - "fieldsplit_vel_pc_gamg_type": "agg", - "fieldsplit_vel_pc_gamg_repartition": "true", - "fieldsplit_vel_pc_mg_type": "additive", - "fieldsplit_vel_pc_gamg_agg_nsmooths": "2", - "fieldsplit_vel_mg_levels_ksp_max_it": "3", - "fieldsplit_vel_mg_levels_ksp_converged_maxits": "true", }) + cfg.update({f"fieldsplit_vel_{key}": value for key, value + in multigrid_options.gamg_bundle().settings.items()}) else: # Custom-FMG velocity block, wrapped in a short FGMRES — NOT `preonly`. # PCFieldSplit applies the Schur complement S = A11 - A10 A00^-1 A01 @@ -1054,20 +1074,23 @@ def _solve_rotated_iterative(solver, Ahat, bhat, Q, Qt, normal_rows, verbose=Fal A_vv, P_vv = vel_pc.getOperators() vel_pc.reset() vel_pc.setOperators(A_vv, P_vv) - custom_mg._configure_pcmg(vel_pc, custom_Pl) - # The Galerkin-coarsened ROTATED velocity block inherits every - # rigid-rotation nullspace mode of the constrained problem (a - # closed circle: one; a spherical shell: three) — the default - # redundant/LU coarse solve hits a zero pivot (SUBPC_ERROR, - # outer reason -11). SVD is nullspace-robust and the coarse - # level is small; same choice as the native spherical FMG setups. + # coarse="svd": the Galerkin-coarsened ROTATED velocity block + # inherits every rigid-rotation nullspace mode of the constrained + # problem (a closed circle: one; a spherical shell: three), and + # the default redundant/LU coarse solve hits a zero pivot there + # (SUBPC_ERROR, outer reason -11). Everything else in the bundle + # is identical to the native and standard custom-P routes. + custom_mg._configure_pcmg(vel_pc, custom_Pl, coarse="svd") + vel_pc.setUp() + # The bundle went into the GLOBAL options DB under the velocity + # sub-PC's prefix, which is derived from `pfx` and so is unique to + # this solve — drop it again now it is consumed, as the `finally` + # below does for the KSP's own keys. vopts = PETSc.Options() vpfx = vel_pc.getOptionsPrefix() or "" - vopts.setValue(vpfx + "mg_coarse_pc_type", "svd") - vopts.delValue(vpfx + "mg_coarse_redundant_pc_type") - vel_pc.setFromOptions() - vel_pc.setUp() - vopts.delValue(vpfx + "mg_coarse_pc_type") + for key in multigrid_options.geometric_mg_bundle( + coarse="svd").settings: + vopts.delValue(vpfx + key) # Constant-pressure nullspace on the Schur COMPLEMENT (enclosed # domains): the IS-built fieldsplit does not propagate the coupled # nullspace to the inner Schur solve, which is otherwise singular diff --git a/tests/parallel/test_1064_rotated_freeslip_parallel.py b/tests/parallel/test_1064_rotated_freeslip_parallel.py index 61c630ef..1ebe04ed 100644 --- a/tests/parallel/test_1064_rotated_freeslip_parallel.py +++ b/tests/parallel/test_1064_rotated_freeslip_parallel.py @@ -236,16 +236,23 @@ def fit(mask): ) -def _annulus_fmg_diagnostics(): +def _annulus_fmg_diagnostics(mesh_owned=False): """Annulus radial free-slip whose velocity block is CUSTOM GEOMETRIC FMG on a - nested hierarchy (coarse annulus -> refine -> refine), via set_custom_fmg — no - GAMG, no direct solve. Returns (velocity L2, leakage L2 on Lower, on Upper).""" + nested hierarchy (coarse annulus -> refine -> refine) — no GAMG, no direct + solve. Returns (velocity L2, leakage L2 on Lower, on Upper). + + ``mesh_owned`` selects how the hierarchy reaches the solver: ``False`` is an + explicit ``set_custom_fmg`` registration, ``True`` is the MESH-owned coarse + tail an ``adapt()`` child carries. Both must produce the same solve (#467).""" RI, RO = 0.5, 1.0 m0 = uw.meshing.Annulus(radiusInner=RI, radiusOuter=RO, cellSize=0.2, qdegree=3) dm1 = m0.dm.refine() dm2 = dm1.refine() coarse = [m0, _wrap(dm1, m0)] fine = _wrap(dm2, m0) + if mesh_owned: + fine._custom_mg_coarse_meshes = coarse + fine._custom_mg_builder = "barycentric" x, y = fine.X r = sympy.sqrt(x**2 + y**2) @@ -264,10 +271,14 @@ def _annulus_fmg_diagnostics(): s.saddle_preconditioner = 1.0 s.petsc_use_pressure_nullspace = True s.petsc_options["snes_type"] = "ksponly" - custom_mg.set_custom_fmg(s, coarse, builder="barycentric", field_id=0) + if not mesh_owned: + custom_mg.set_custom_fmg(s, coarse, builder="barycentric", field_id=0) s.solve() assert s._rotated_freeslip_info["ksp_reason"] > 0, "custom-FMG rotated solve did not converge" + assert s._rotated_freeslip_info["velocity_pc"] == "custom-FMG", ( + "rotated velocity block fell back to " + f"{s._rotated_freeslip_info['velocity_pc']} instead of geometric MG") L2 = float(np.sqrt(uw.maths.Integral(fine, v.sym.dot(v.sym)).evaluate())) vr = v.sym[0] * x / r + v.sym[1] * y / r leak_lo = float(np.sqrt(uw.maths.BdIntegral( @@ -424,6 +435,25 @@ def test_rotated_freeslip_annulus_fmg_partition_independent(): f"{leak_up_ref} vs {leak_up}") +def test_rotated_freeslip_mesh_owned_fmg_pickup(): + """#467, in parallel: a coarse tail owned by the MESH — what ``adapt()`` + leaves on a refinement child — must drive geometric MG under rotated + free-slip, on every partition. The rotated path used to consult only an + explicit ``set_custom_fmg`` registration and fell back to GAMG silently. + + The transfers are built cross-partition, so this is not implied by the serial + pickup test (``tests/test_1021_mg_option_bundle.py``); the tail is attached by + hand because what is under test is whether the rotated dispatch consults it, + not how ``adapt()`` produces it.""" + L2, leak_lo, leak_up = _annulus_fmg_diagnostics(mesh_owned=True) + L2_ref, leak_lo_ref, leak_up_ref = GOLDEN_ANNULUS_FMG + assert np.isclose(L2, L2_ref, rtol=1e-6, atol=0), ( + f"mesh-owned FMG annulus velocity L2 differs serial vs np={uw.mpi.size}: " + f"{L2_ref} vs {L2}") + assert np.isclose(leak_lo, leak_lo_ref, rtol=1e-4, atol=0) + assert np.isclose(leak_up, leak_up_ref, rtol=1e-4, atol=0) + + def test_rotated_freeslip_spherical3d_partition_independent(): """3D spherical shell (free-slip inner+outer, all three rotation nullspace modes): the parallel solve reproduces the serial velocity L2, converges, and diff --git a/tests/test_1021_mg_option_bundle.py b/tests/test_1021_mg_option_bundle.py new file mode 100644 index 00000000..53b8983d --- /dev/null +++ b/tests/test_1021_mg_option_bundle.py @@ -0,0 +1,240 @@ +"""The geometric-multigrid option bundle has ONE owner, and every route reads it. + +Three routes reach a multigrid velocity block — native (DMPlex refinement), +custom-P on the standard solve path, and custom-P through the rotated free-slip +path — and they are the same preconditioner, not alternatives: custom-P is +mandatory wherever native cannot go (rotated BCs, ``adapt()`` children). They +drifted apart because the bundle was written in two places (#468): the custom-P +route ran richardson smoothing at an iteration count nobody set, inherited from +whatever had last written that options prefix (3 left over from the GAMG bundle +on the standard path, PETSc's own default of 2 on the rotated path). + +These tests read the configuration back off the LIVE PETSc objects after setup, +not out of the options database, so they check what actually runs. + +The one legitimate per-route difference is the coarse solve: the +Galerkin-coarsened *rotated* velocity block inherits the rigid-rotation null +space, so redundant/LU hits a zero pivot there and it uses SVD (the #306 fix). +That difference is asserted, not tolerated — if it disappears, the rotated coarse +solve has silently reverted. + +Also here: the rotated path must pick up a MESH-OWNED hierarchy (the coarse tail +an ``adapt()`` child carries), which it previously ignored, silently solving on +GAMG (#467). +""" +import pytest +import sympy + +import underworld3 as uw +from underworld3.utilities import custom_mg, multigrid_options, rotated_bc + +pytestmark = [pytest.mark.level_2, pytest.mark.tier_b] + +RES = 0.15 +R_IN, R_OUT = 0.5, 1.0 + + +# --------------------------------------------------------------------------- # +# The owner itself — no solve needed +# --------------------------------------------------------------------------- # +def test_every_bundle_sets_the_smoother_iteration_count(): + """The drift mechanism: a bundle that omits ``mg_levels_ksp_max_it`` inherits + it from whatever last wrote that prefix. Every bundle must SET it.""" + bundles = [multigrid_options.gamg_bundle()] + bundles += [multigrid_options.geometric_mg_bundle(coarse=c) + for c in multigrid_options.GEOMETRIC_MG_COARSE_SOLVERS] + for bundle in bundles: + assert "mg_levels_ksp_max_it" in bundle.settings + assert "pc_type" in bundle.settings + + +def test_bundles_clear_each_others_keys(): + """Bundles share an options prefix, so switching between them must leave no + key behind — every key a bundle does not set, a sibling does, and it must + appear in that bundle's stale list.""" + bundles = [multigrid_options.gamg_bundle()] + bundles += [multigrid_options.geometric_mg_bundle(coarse=c) + for c in multigrid_options.GEOMETRIC_MG_COARSE_SOLVERS] + owned = set().union(*(set(b.settings) for b in bundles)) + for bundle in bundles: + assert set(bundle.stale) == owned - set(bundle.settings) + + +# --------------------------------------------------------------------------- # +# The three live routes +# --------------------------------------------------------------------------- # +def _stokes(mesh, tag, rotated): + """A buoyancy-driven annulus Stokes solve: fixed inner boundary, and an outer + boundary that is either an ordinary essential BC or rotated free-slip.""" + x, y = mesh.X + r = sympy.sqrt(x**2 + y**2) + rhat = sympy.Matrix([[x / r, y / r]]) + v = uw.discretisation.MeshVariable(f"V{tag}", mesh, mesh.dim, degree=2, + continuous=True) + p = uw.discretisation.MeshVariable(f"P{tag}", mesh, 1, degree=1, continuous=True) + s = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + s.constitutive_model = uw.constitutive_models.ViscousFlowModel + s.constitutive_model.Parameters.shear_viscosity_0 = 1.0 + blob = sympy.exp(-(((x - 0.75) ** 2 + y**2) / 0.05)) + s.bodyforce = sympy.Matrix([[50.0 * blob * x / r, 50.0 * blob * y / r]]) + s.add_essential_bc((0.0, 0.0), "Lower") + if rotated: + s.add_rotated_freeslip_bc(0.0, "Upper", normal=rhat) + else: + s.add_essential_bc((sympy.oo, 0.0), "Upper") + s.petsc_use_pressure_nullspace = True + s.tolerance = 1.0e-8 + return s + + +def _mg_config(vel_pc): + """Read the multigrid configuration off a live velocity sub-PC. The finest + smoother is representative; level 0 is the coarse solve.""" + assert vel_pc.getType() == "mg" + nlev = vel_pc.getMGLevels() + smoother = vel_pc.getMGSmoother(nlev - 1) + coarse = vel_pc.getMGCoarseSolve() + return { + "levels": nlev, + "mg_type": str(vel_pc.getMGType()), + "smoother_ksp": smoother.getType(), + "smoother_pc": smoother.getPC().getType(), + "smoother_max_it": smoother.getTolerances()[3], + "smoother_norm": str(smoother.getNormType()), + "coarse_ksp": coarse.getType(), + "coarse_pc": coarse.getPC().getType(), + } + + +def _annulus(cell_size): + return uw.meshing.Annulus(radiusInner=R_IN, radiusOuter=R_OUT, + cellSize=cell_size, qdegree=3) + + +def _native_config(): + """Route A: native geometric FMG on a refined DMPlex hierarchy.""" + s = _stokes(uw.meshing.Annulus(radiusInner=R_IN, radiusOuter=R_OUT, + cellSize=2 * RES, qdegree=3, refinement=1), + "nat", rotated=False) + s.solve() + return _mg_config(s.snes.getKSP().getPC().getFieldSplitSubKSP()[0].getPC()) + + +def _custom_standard_config(): + """Route B: custom-P FMG reached through the standard solve path.""" + s = _stokes(_annulus(RES), "std", rotated=False) + custom_mg.set_custom_fmg(s, [_annulus(2 * RES)], field_id=0) + s.solve() + return _mg_config(s.snes.getKSP().getPC().getFieldSplitSubKSP()[0].getPC()) + + +def _custom_rotated_config(monkeypatch): + """Route C: custom-P FMG reached through the rotated free-slip path. + + The rotated KSP is self-contained and is destroyed at the end of the solve, + so the configuration is captured from inside the solve — there is no live PC + left to interrogate afterwards. + """ + s = _stokes(_annulus(RES), "rot", rotated=True) + custom_mg.set_custom_fmg(s, [_annulus(2 * RES)], field_id=0) + + real_solve = rotated_bc._solve_rotated_iterative + captured = {} + + def capture(solver, Ahat, bhat, Q, Qt, normal_rows, **kw): + result = real_solve(solver, Ahat, bhat, Q, Qt, normal_rows, **kw) + captured.setdefault("config", _mg_config( + result[2]["pc"].getFieldSplitSubKSP()[0].getPC())) + return result + + monkeypatch.setattr(rotated_bc, "_solve_rotated_iterative", capture) + s.solve() + assert s._rotated_freeslip_info["velocity_pc"] == "custom-FMG" + return captured["config"] + + +def test_all_three_routes_share_one_bundle(monkeypatch): + """Native, standard custom-P and rotated custom-P must smooth identically — + same Krylov smoother, same preconditioner, same iteration count, same cycle. + Only the coarse solve legitimately differs.""" + native = _native_config() + standard = _custom_standard_config() + rotated = _custom_rotated_config(monkeypatch) + + shared = ("mg_type", "smoother_ksp", "smoother_pc", "smoother_max_it", + "smoother_norm", "coarse_ksp") + for key in shared: + assert native[key] == standard[key] == rotated[key], ( + f"route drift on {key!r}: native={native[key]!r} " + f"standard={standard[key]!r} rotated={rotated[key]!r}") + + # The measured bundle, not just self-consistency: gmres/sor at four + # iterations with no norm computation (see multigrid_options for why). + assert native["smoother_ksp"] == "gmres" + assert native["smoother_pc"] == "sor" + assert native["smoother_max_it"] == 4 + + # The ONE deliberate per-route difference (#306): the rotated coarse operator + # inherits the rigid-rotation null space, where redundant/LU hits a zero pivot. + assert native["coarse_pc"] == "redundant" + assert standard["coarse_pc"] == "redundant" + assert rotated["coarse_pc"] == "svd" + + +def test_rotated_picks_up_a_mesh_owned_hierarchy(): + """#467: a coarse tail owned by the MESH — what ``adapt()`` leaves on a + refinement child — must drive geometric MG under rotated free-slip. It used + to be ignored, and the solve fell back to GAMG silently. + + The tail is attached by hand rather than by running ``adapt()``: it is the + same attribute the child carries, and what is under test is whether the + rotated dispatch consults it. + """ + mesh = _annulus(RES) + mesh._custom_mg_coarse_meshes = [_annulus(2 * RES)] + mesh._custom_mg_builder = "barycentric" + + s = _stokes(mesh, "own", rotated=True) + s.solve() + + assert s._rotated_freeslip_info["velocity_pc"] == "custom-FMG" + assert s._rotated_freeslip_info["velocity_pc_type"] == "mg" + + +def test_rotated_fmg_survives_repeated_newton_increments(): + """The rotated path applies its bundle under a per-solve options prefix and + then drops the keys again, so the global database stays bounded under + time-stepping. That is only safe if nothing re-reads them: a later + ``setFromOptions`` on the velocity sub-PC would find ``pc_type`` gone and + silently abandon the multigrid. Several Newton increments, each re-solving + on a refreshed operator, must all keep ``pc=mg``.""" + mesh = _annulus(RES) + s = _stokes(mesh, "nl", rotated=True) + + # power-law viscosity: a genuinely nonlinear solve, several increments + x, y = mesh.X + v = s.Unknowns.u + grad = sympy.Matrix([[v.sym[0].diff(x), v.sym[0].diff(y)], + [v.sym[1].diff(x), v.sym[1].diff(y)]]) + edot = 0.5 * (grad + grad.T) + eII = sympy.sqrt(0.5 * (edot[0, 0] ** 2 + edot[1, 1] ** 2) + + edot[0, 1] ** 2 + 1.0e-12) + s.constitutive_model.Parameters.shear_viscosity_0 = eII ** (1.0 / 3.0 - 1.0) + s.consistent_jacobian = True + s.tolerance = 1.0e-7 + custom_mg.set_custom_fmg(s, [_annulus(2 * RES)], field_id=0) + s.solve() + + info = s._rotated_freeslip_info + assert len(info["ksp_its"]) > 1, "not a multi-increment solve — test is vacuous" + assert info["velocity_pc"] == "custom-FMG" + assert info["velocity_pc_type"] == "mg" + + +def test_rotated_without_a_hierarchy_falls_back_to_gamg(): + """The negative control for the test above: no hierarchy anywhere still + reports GAMG, so ``custom-FMG`` above is a real pickup and not a label that + is always set.""" + s = _stokes(_annulus(RES), "none", rotated=True) + s.solve() + assert s._rotated_freeslip_info["velocity_pc"] == "GAMG" From 2550d9a43a603b9d725dbbf07e9094d604b53920 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Wed, 29 Jul 2026 19:02:34 +1000 Subject: [PATCH 2/5] Correct the smoother's recorded mechanism: the velocity block stays symmetric under the consistent tangent MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The comment this PR moved into utilities/multigrid_options.py carried a claim we inherited without checking, and which the Layer 3 section of design/nonlinear-solver-homotopy-warmstart.md states outright: that richardson degrades because the consistent-Newton tangent makes the velocity block NON-SYMMETRIC. It does not. For an isotropic eta(eps_II) the Newton term is 2 (eta'/eps_II) (edot(u):edot(du)) (edot(u):edot(v)) a rank-one outer product a (x) a with a = edot(u) — symmetric under swapping du and v. eta depends on grad v only THROUGH its symmetric part, so both factors project onto the same tensor and nothing is left to break the symmetry. Measured as ||A - A^T||_F / ||A||_F on the assembled velocity block (velocity_block_symmetry.py in the rotated_pmat_audit study), against a linear calibration of 5.1e-17, with a control on the size of the Newton contribution so a null result cannot be vacuous: rheology asym Newton term linear (calibration) 5.1e-17 — shear-thinning, isotropic 5.3e-17 2.1e-01 power-law, isotropic 5.3e-17 4.0e-01 Drucker-Prager, pressure-dependent yield 5.1e-17 5.1e-02 transverse isotropic, PICARD tangent 7.2e-02 — transverse isotropic, Newton tangent 8.9e-02 2.0e-01 The consistent tangent contributed 5-40% of ||A|| in every nonlinear case and the block stayed symmetric anyway. Non-symmetry appears only under transverse isotropy — and there it appears under the PICARD tangent too, where a frozen-coefficient form int edot(du):C:edot(v) is symmetric by construction. That is a defect signature, not a property of the tangent: linear TI with anisotropy fully active is symmetric at 6.3e-17. Reported on #457. The Spiegelman-notch contraction measurement (0.75 richardson vs 0.56 gmres at four iterations, margin growing with depth) is unaffected — it is a measurement, not a derivation. The mechanism that fits it AND fits the cost we measured on well-conditioned problems is operator CONDITIONING: at eta contrast 1e26 the SOR-preconditioned spectrum is spread far enough that a stationary iteration stalls where a Krylov smoother adapts its polynomial. Also records the cost, which was previously unmeasured on this route: nested annulus, LINEAR (symmetric) velocity block, eta contrast 1e6, outer KSP timed in isolation, gmres/4 against richardson/3 — iterations x1.45 / x1.60 / x1.20 at 2 / 3 / 4 levels, wall clock x0.86 / x0.77 / x0.55. gmres wins iterations everywhere and loses wall clock everywhere. It stays the default on both routes because the failure it avoids is worse than the cost it carries, and it carries that cost exactly where the problems are easy; the number is recorded so the guardrail can be removed deliberately. No behavioural change — comments, docstrings and docs/developer/subsystems/solvers.md. A TODO(DESIGN) marks the design doc's Layer 3 claim rather than rewriting it. Underworld development team with AI support from Claude Code --- docs/developer/subsystems/solvers.md | 73 +++++++++++++++++-- .../utilities/multigrid_options.py | 49 +++++++++++-- 2 files changed, 109 insertions(+), 13 deletions(-) diff --git a/docs/developer/subsystems/solvers.md b/docs/developer/subsystems/solvers.md index 3f492df1..ba2e086b 100644 --- a/docs/developer/subsystems/solvers.md +++ b/docs/developer/subsystems/solvers.md @@ -371,12 +371,73 @@ Two consequences worth stating outright: last writer left. - **The measured smoother is `gmres`+`sor` at 4 iterations.** Chebyshev needs eigenvalue estimates of the smoothed operator, which are fragile on the - indefinite / variable-viscosity velocity block. Richardson is stationary and - degrades on the non-symmetric operator the consistent-Newton tangent produces: - measured on the Spiegelman notch (Drucker–Prager, η contrast 1e26, nested - 4-level hierarchy) the per-cycle contraction is ρ = 0.75 richardson against - 0.56 gmres at the *same* four iterations, and the gmres margin **grows with - depth** (5% at 3 levels, 25% at 4). + indefinite / variable-viscosity velocity block. Against richardson, measured on + the Spiegelman notch (Drucker–Prager, η contrast 1e26, nested 4-level + hierarchy): per-cycle contraction ρ = 0.75 richardson against 0.56 gmres at the + *same* four iterations, the margin growing with depth (5% at 3 levels, 25% at 4). + +### Why richardson loses — the recorded reason is wrong + +Both this bundle's comment and the Layer 3 section of +`design/nonlinear-solver-homotopy-warmstart.md` said richardson degrades because +the consistent-Newton tangent makes the velocity block **non-symmetric**. That +mechanism does not survive measurement. + +For an isotropic η(ε̇_II) the Newton term is + +$$2\,\frac{\eta'}{\dot\varepsilon_{II}}\big(\dot\varepsilon(u):\dot\varepsilon(\delta u)\big)\big(\dot\varepsilon(u):\dot\varepsilon(v)\big)$$ + +— a rank-one outer product $a\otimes a$ with $a = \dot\varepsilon(u)$, which is +symmetric. η depends on ∇v only *through* its symmetric part, so both factors +project onto the same tensor and there is nothing left to break the symmetry. +Measured as ‖A−Aᵀ‖_F/‖A‖_F on the assembled velocity block, against a linear +calibration of 5e-17: + +| rheology | ‖A−Aᵀ‖/‖A‖ | Newton term ‖A_N−A_P‖/‖A‖ | +|---|---|---| +| linear (calibration) | 5.1e-17 | — | +| shear-thinning, isotropic | 5.3e-17 | 2.1e-01 | +| power-law, isotropic | 5.3e-17 | 4.0e-01 | +| Drucker–Prager, pressure-dependent yield | 5.1e-17 | 5.1e-02 | +| transverse isotropic, Picard tangent | **7.2e-02** | — | +| transverse isotropic, Newton tangent | **8.9e-02** | 2.0e-01 | + +The third column is the control: if the consistent tangent had contributed +nothing, the symmetry reading would be vacuous. It contributed 5–40% of ‖A‖ in +every case, and the block stayed symmetric anyway. + +Non-symmetry *does* appear under transverse isotropy — but it appears under the +**Picard** tangent too, where a frozen-coefficient form ∫ ε̇(δu):C:ε̇(v) is +symmetric by construction. That is a defect signature rather than a property of +the tangent (see #457, and the linear TI rows, which are symmetric at 6.3e-17 +even with anisotropy fully active). + +The mechanism that does fit everything measured is **operator conditioning**: at +η contrast 1e26 the SOR-preconditioned spectrum is spread far enough that a +stationary iteration stalls where a Krylov smoother adapts its polynomial. It +predicts the notch result, and it predicts the cost below on well-conditioned +problems. + +### The smoother is not free + +On a *well-conditioned* problem gmres/4 is slower than richardson/3 despite +converging in fewer iterations. Nested annulus, linear (symmetric) velocity +block, η contrast 1e6, outer KSP timed in isolation: + +| depth | iterations | wall clock | +|---|---|---| +| 2 levels | ×1.45 | **×0.86** | +| 3 levels | ×1.60 | **×0.77** | +| 4 levels | ×1.20 | **×0.55** | + +The extra sweep costs ~3%; gmres itself costs ~40% per cycle. Custom-P builds its +Galerkin coarse operators from barycentric transfers, denser than the native +nested ones, so the smoother is a larger share of the cycle on that route. + +It stays the default on both routes because the failure it avoids (a stationary +smoother stalling on a high-contrast operator) is worse than the cost it carries, +and it carries that cost exactly where the problems are easy. Removing a guardrail +is the caller's decision; the number is recorded so the decision can be made. ### The one legitimate per-route difference diff --git a/src/underworld3/utilities/multigrid_options.py b/src/underworld3/utilities/multigrid_options.py index 2245a371..cd2e9d27 100644 --- a/src/underworld3/utilities/multigrid_options.py +++ b/src/underworld3/utilities/multigrid_options.py @@ -97,15 +97,32 @@ def _geometric_mg_settings(coarse): # two-level cycle is a coarse-grid correction, not a V-cycle, and is not # worth special-casing). Chebyshev needs eigenvalue estimates of the # smoothed operator, which are fragile on the indefinite / - # variable-viscosity velocity block and diverge. Richardson is stationary - # and degrades on the NON-SYMMETRIC operator produced by the - # consistent-Newton tangent. Measured on the Spiegelman notch - # (Drucker-Prager, eta contrast 1e26) over a nested 4-level hierarchy: + # variable-viscosity velocity block and diverge. + # + # The measurement, which is what this setting rests on: Spiegelman notch + # (Drucker-Prager, eta contrast 1e26) over a nested 4-level hierarchy, # contraction per V-cycle rho = 0.75 (richardson) vs 0.56 (gmres) at the - # SAME four smoother iterations -- and the gmres margin GROWS with depth - # (5% at 3 levels, 25% at 4), because deeper cycles apply the smoother on - # more coarse operators. Four iterations, not more: per unit work gmres/4 + # SAME four smoother iterations, the margin growing with depth (5% at 3 + # levels, 25% at 4). Four iterations, not more: per unit work gmres/4 # (rho^(1/4) = 0.87) beats gmres/8 (0.91). + # + # TODO(DESIGN): the MECHANISM on record for that result is wrong. Both + # this comment and design/nonlinear-solver-homotopy-warmstart.md (Layer 3) + # said richardson degrades because the consistent-Newton tangent makes the + # velocity block NON-SYMMETRIC. Measured (velocity_block_symmetry.py in the + # rotated_pmat_audit study), that block stays symmetric to machine precision + # -- 5e-17, against a linear calibration of 5e-17 -- under shear-thinning, + # power-law AND Drucker-Prager with pressure-dependent yield, with a control + # confirming the Newton term contributed 5-40% of ||A||. For an isotropic + # eta(eps_II) the Newton term is a rank-one outer product a (x) a with + # a = eps_dot(u), which is symmetric. Non-symmetry does appear under + # transverse isotropy, but it appears under the PICARD tangent there too + # (7e-2), where a frozen-coefficient form is symmetric by construction -- + # i.e. a defect signature, see #457. The likely real mechanism is operator + # CONDITIONING: at eta contrast 1e26 the SOR-preconditioned spectrum is + # spread far enough that a stationary iteration stalls where a Krylov + # smoother adapts its polynomial. That also predicts the cost measured on + # well-conditioned problems (see the note on geometric_mg_bundle). "mg_levels_ksp_type": "gmres", "mg_levels_pc_type": "sor", # SET the count, never inherit it. PCMG's own default is 2 and the GAMG @@ -183,6 +200,24 @@ def geometric_mg_bundle(coarse="redundant"): ------- MGBundle Settings and stale keys; see :meth:`MGBundle.apply`. + + Notes + ----- + The ``gmres``/4 smoother is not free on a well-conditioned problem. Measured on + a nested annulus with a linear (symmetric) velocity block at :math:`\eta` + contrast 1e6, outer KSP timed in isolation, ``gmres``/4 against + ``richardson``/3: it wins on iterations at every depth (×1.45, ×1.60, ×1.20 at + 2, 3, 4 levels) and **loses on wall clock at every depth** (×0.86, ×0.77, + ×0.55) — the extra sweep costs ~3% and ``gmres`` itself ~40% per cycle. Custom-P + builds its Galerkin coarse operators from barycentric transfers, denser than + the native nested ones, so the smoother is a larger share of the cycle here. + + It is the default anyway, on both routes, because the failure it avoids is + worse than the cost it carries: a stationary smoother stalls where a Krylov one + adapts, and it does so exactly on the badly-conditioned high-contrast problems + this code exists for (the Spiegelman-notch measurement above). Removing a + guardrail is the caller's decision to take; the number is quoted so that + decision can be made. """ return _bundle(_geometric_mg_settings(coarse)) From 5a390efe12e3bee730326f23b36e21b7e39f560a Mon Sep 17 00:00:00 2001 From: lmoresi Date: Thu, 30 Jul 2026 11:12:40 +1000 Subject: [PATCH 3/5] The bundle owner must not overwrite what the user set (audit F-1) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The multigrid bundle was applied wholesale on every rebuild. The only escape was the `_pc_user_override` latch, and that latch keys on `pc_type` ALONE — so of the bundle's ten keys, setting any of nine was silently discarded, and setting the tenth rescued all ten. Measured, asking for chebyshev / max_it 6 / coarse svd on the Stokes velocity block: what the user set result the three bundle keys only OVERWRITTEN (gmres/4/redundant) the three keys + fieldsplit_velocity_pc_type=mg honoured the three keys + preconditioner="fmg" OVERWRITTEN Note the inversion in the last row: explicit `preconditioner="fmg"` skips the latch entirely, because that mode "always applies". The more clearly a user stated their intent, the less control they had. And the workaround for the first row — "also set pc_type to the value it already has" — is undiscoverable. Making this PR the single owner of the bundle would have entrenched all of it, which is why the fix belongs here. `MGBundle.apply` now takes an `owned` record and leaves alone any key whose current value is not the one UW3 last wrote — in the stale-key clear too, since a user-set key is not ours to remove. Ownership is RECORDED, never inferred from the value: inference fails as soon as a second internal writer touches the same key, which is exactly how the `tolerance` and `strategy` setters defeated an earlier attempt at this (#477). Every internal writer of a bundle key therefore goes through the new `SolverBaseClass._push_managed_option`, so anything unrecorded is the user's by construction. The record is keyed by the GLOBALLY-QUALIFIED option name. The solver writes through a prefixed Options view (`Solver_N_`) while `custom_mg._configure_pcmg` reads the global database using the live PC's own full prefix; recording the unqualified name made every key look user-owned over there and the bundle silently stopped applying. That regression was caught by the defaults arm, not by any test — which is the whole argument for keeping a defaults arm. Tests: `test_user_set_bundle_keys_are_honoured` (parametrised over auto and explicit `fmg`) asserts the two keys the user set survive a rebuild AND that the key they left alone still carries the measured default — respecting one key must not abandon the bundle. `test_unset_bundle_keys_keep_the_managed_defaults` is its control: without it the first test could pass trivially by never applying the bundle, which is precisely how the first attempt failed. Validation: fmg_config_diff shows the three routes still agreeing on every row but the deliberate rotated `svd`. 68 passed across test_1014/1015/1016/1017/1018/1020 and test_0753/0835; `level_1 and tier_a` 526 passed with the two pre-existing #470 failures; parallel np2 and np4 green (6 / 9 / 2 passed). Found by the solver configuration audit, `docs/reviews/2026-07/solver-configuration-reachability-audit.md`. Underworld development team with AI support from Claude Code --- .../cython/petsc_generic_snes_solvers.pyx | 114 ++++++++++++------ src/underworld3/utilities/custom_mg.py | 14 ++- .../utilities/multigrid_options.py | 58 ++++++++- tests/test_1021_mg_option_bundle.py | 56 +++++++++ 4 files changed, 195 insertions(+), 47 deletions(-) diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index 21c3520b..ced17dd6 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -161,11 +161,37 @@ class SolverBaseClass(uw_object): # 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 + # Every multigrid option value UW3 itself has written on the managed block, + # keyed by full option name. This is what lets the bundle honour a + # user-set smoother while still managing the keys the user left alone: a + # present key whose value is not the one we recorded writing is theirs. + # Ownership is RECORDED, never inferred from the value — inference fails + # the moment a second internal writer touches the same key, which is how + # the `tolerance` and `strategy` setters defeated an earlier attempt (#477). + # Internal writers therefore go through _push_managed_option(). + self._managed_pc_options = {} # Custom multigrid prolongation hierarchy (see set_custom_mg / # utilities.custom_mg). None => standard FMG/GAMG path, unchanged. self._custom_mg = None + def _push_managed_option(self, key, value): + """Write a PETSc option UW3 owns, recording that we wrote it. + + Use this for any option a solver sets on its own behalf that the multigrid + bundles also write (``utilities.multigrid_options``). A plain + ``self.petsc_options[key] = value`` is indistinguishable from a user's own + write, and the bundle would then back off from a key nobody asked for. + """ + self.petsc_options[key] = value + # Key the record by the GLOBAL option name. `self.petsc_options` is a + # prefixed view (`Solver_N_`), but custom_mg._configure_pcmg reads the + # global database using the live PC's own full prefix — so an unqualified + # record makes every key look user-owned over there and the bundle + # silently stops applying. + self._managed_pc_options[self.petsc_options_prefix + key] = \ + multigrid_options.option_string(value) + @property def consistent_jacobian(self): r"""Jacobian tangent selection: ``False`` | ``True`` | ``"continuation"``. @@ -625,8 +651,18 @@ class SolverBaseClass(uw_object): # be configured from three places (#468). Coarse solve: the native # hierarchy is not rotated, so its coarse operator carries no inherited # null space and redundant+LU is right here. + # `_managed_pc_options` makes the bundle respect a key the USER set while + # still managing the ones they left alone. Before this, the bundle was + # applied wholesale on every rebuild and the only escape was the + # `_pc_user_override` latch above — which keys on `pc_type` ALONE, so a user + # who set (say) `mg_levels_ksp_max_it` had it silently discarded unless they + # also set `pc_type` to the value it already had. Explicit + # `preconditioner="fmg"` was worse: it skips the latch entirely, so the + # clearer the request the less control it carried. if want_fmg: - multigrid_options.geometric_mg_bundle().apply(opts, prefix) + multigrid_options.geometric_mg_bundle().apply( + PETSc.Options(), self.petsc_options_prefix + prefix, + owned=self._managed_pc_options) self._pc_managed_value = "mg" else: if self._preconditioner == "fmg" and n_levels <= 1 and uw.mpi.rank == 0: @@ -637,7 +673,9 @@ class SolverBaseClass(uw_object): f"mesh with refinement >= 1 to enable geometric multigrid.", stacklevel=2, ) - multigrid_options.gamg_bundle().apply(opts, prefix) + multigrid_options.gamg_bundle().apply( + PETSc.Options(), self.petsc_options_prefix + prefix, + owned=self._managed_pc_options) self._pc_managed_value = "gamg" def _enforce_galerkin_for_geometric_mg(self): @@ -2944,16 +2982,16 @@ class SNES_Scalar(SolverBaseClass): self.petsc_options["snes_type"] = "newtonls" self.petsc_options["ksp_type"] = "gmres" - self.petsc_options["pc_type"] = "gamg" - self.petsc_options["pc_gamg_type"] = "agg" - self.petsc_options["pc_gamg_repartition"] = True - self.petsc_options["pc_mg_type"] = "additive" - self.petsc_options["pc_gamg_agg_nsmooths"] = 2 - self.petsc_options["mg_levels_ksp_max_it"] = 3 - self.petsc_options["mg_levels_ksp_converged_maxits"] = None + self._push_managed_option("pc_type", "gamg") + self._push_managed_option("pc_gamg_type", "agg") + self._push_managed_option("pc_gamg_repartition", True) + self._push_managed_option("pc_mg_type", "additive") + self._push_managed_option("pc_gamg_agg_nsmooths", 2) + self._push_managed_option("mg_levels_ksp_max_it", 3) + self._push_managed_option("mg_levels_ksp_converged_maxits", None) self.petsc_options["snes_rtol"] = 1.0e-4 - self.petsc_options["mg_levels_ksp_max_it"] = 3 + self._push_managed_option("mg_levels_ksp_max_it", 3) if self.verbose == True: self.petsc_options["ksp_monitor"] = None @@ -3857,14 +3895,14 @@ class SNES_Vector(SolverBaseClass): self.petsc_options["snes_type"] = "newtonls" self.petsc_options["ksp_rtol"] = 1.0e-3 self.petsc_options["ksp_type"] = "gmres" - self.petsc_options["pc_type"] = "gamg" - self.petsc_options["pc_gamg_type"] = "agg" - self.petsc_options["pc_gamg_repartition"] = True - self.petsc_options["pc_mg_type"] = "additive" - self.petsc_options["pc_gamg_agg_nsmooths"] = 2 + self._push_managed_option("pc_type", "gamg") + self._push_managed_option("pc_gamg_type", "agg") + self._push_managed_option("pc_gamg_repartition", True) + self._push_managed_option("pc_mg_type", "additive") + self._push_managed_option("pc_gamg_agg_nsmooths", 2) self.petsc_options["snes_rtol"] = 1.0e-3 - self.petsc_options["mg_levels_ksp_max_it"] = 3 - self.petsc_options["mg_levels_ksp_converged_maxits"] = None + self._push_managed_option("mg_levels_ksp_max_it", 3) + self._push_managed_option("mg_levels_ksp_converged_maxits", None) if self.verbose == True: self.petsc_options["ksp_monitor"] = None @@ -4869,14 +4907,14 @@ class SNES_MultiComponent(SolverBaseClass): self.petsc_options["snes_type"] = "newtonls" self.petsc_options["ksp_rtol"] = 1.0e-3 self.petsc_options["ksp_type"] = "gmres" - self.petsc_options["pc_type"] = "gamg" - self.petsc_options["pc_gamg_type"] = "agg" - self.petsc_options["pc_gamg_repartition"] = True - self.petsc_options["pc_mg_type"] = "additive" - self.petsc_options["pc_gamg_agg_nsmooths"] = 2 + self._push_managed_option("pc_type", "gamg") + self._push_managed_option("pc_gamg_type", "agg") + self._push_managed_option("pc_gamg_repartition", True) + self._push_managed_option("pc_mg_type", "additive") + self._push_managed_option("pc_gamg_agg_nsmooths", 2) self.petsc_options["snes_rtol"] = 1.0e-3 - self.petsc_options["mg_levels_ksp_max_it"] = 3 - self.petsc_options["mg_levels_ksp_converged_maxits"] = None + self._push_managed_option("mg_levels_ksp_max_it", 3) + self._push_managed_option("mg_levels_ksp_converged_maxits", None) if self.verbose == True: self.petsc_options["ksp_monitor"] = None @@ -5736,13 +5774,13 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): self.petsc_options[f"fieldsplit_{v_name}_ksp_type"] = "fgmres" self.petsc_options[f"fieldsplit_{v_name}_ksp_max_it"] = 200 self.petsc_options[f"fieldsplit_{v_name}_ksp_rtol"] = self._tolerance * 0.1 - self.petsc_options[f"fieldsplit_{v_name}_pc_type"] = "gamg" - self.petsc_options[f"fieldsplit_{v_name}_pc_gamg_type"] = "agg" - self.petsc_options[f"fieldsplit_{v_name}_pc_gamg_repartition"] = True - self.petsc_options[f"fieldsplit_{v_name}_pc_mg_type"] = "additive" - self.petsc_options[f"fieldsplit_{v_name}_pc_gamg_agg_nsmooths"] = 2 - self.petsc_options[f"fieldsplit_{v_name}_mg_levels_ksp_max_it"] = 3 - self.petsc_options[f"fieldsplit_{v_name}_mg_levels_ksp_converged_maxits"] = None + self._push_managed_option(f"fieldsplit_{v_name}_pc_type", "gamg") + self._push_managed_option(f"fieldsplit_{v_name}_pc_gamg_type", "agg") + self._push_managed_option(f"fieldsplit_{v_name}_pc_gamg_repartition", True) + self._push_managed_option(f"fieldsplit_{v_name}_pc_mg_type", "additive") + self._push_managed_option(f"fieldsplit_{v_name}_pc_gamg_agg_nsmooths", 2) + self._push_managed_option(f"fieldsplit_{v_name}_mg_levels_ksp_max_it", 3) + self._push_managed_option(f"fieldsplit_{v_name}_mg_levels_ksp_converged_maxits", None) # Create this dict self.fields = {} @@ -6338,16 +6376,16 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): # preconditioning and weakly-indefinite coarse operators; issue #147). self.petsc_options[f"fieldsplit_velocity_ksp_type"] = "fgmres" self.petsc_options[f"fieldsplit_velocity_ksp_max_it"] = 200 - self.petsc_options[f"fieldsplit_velocity_pc_type"] = "gamg" - self.petsc_options[f"fieldsplit_velocity_pc_gamg_type"] = "agg" - self.petsc_options[f"fieldsplit_velocity_pc_gamg_repartition"] = True + self._push_managed_option(f"fieldsplit_velocity_pc_type", "gamg") + self._push_managed_option(f"fieldsplit_velocity_pc_gamg_type", "agg") + self._push_managed_option(f"fieldsplit_velocity_pc_gamg_repartition", True) # NOTE: "kaskade" here diverges from the "additive" set by __init__'s # velocity block — a long-standing (possibly unintentional) difference. # Do not change without benchmarking (READ-23 kept the value as-is). - self.petsc_options[f"fieldsplit_velocity_pc_mg_type"] = "kaskade" - self.petsc_options[f"fieldsplit_velocity_pc_gamg_agg_nsmooths"] = 2 - self.petsc_options[f"fieldsplit_velocity_mg_levels_ksp_max_it"] = 3 - self.petsc_options[f"fieldsplit_velocity_mg_levels_ksp_converged_maxits"] = None + self._push_managed_option(f"fieldsplit_velocity_pc_mg_type", "kaskade") + self._push_managed_option(f"fieldsplit_velocity_pc_gamg_agg_nsmooths", 2) + self._push_managed_option(f"fieldsplit_velocity_mg_levels_ksp_max_it", 3) + self._push_managed_option(f"fieldsplit_velocity_mg_levels_ksp_converged_maxits", None) diff --git a/src/underworld3/utilities/custom_mg.py b/src/underworld3/utilities/custom_mg.py index 5d96d272..50c8504f 100644 --- a/src/underworld3/utilities/custom_mg.py +++ b/src/underworld3/utilities/custom_mg.py @@ -595,7 +595,7 @@ def _assert_no_zero_columns_serial(P_csr, level): f"operator would be singular.") -def _configure_pcmg(pc, Ps, coarse="redundant"): +def _configure_pcmg(pc, Ps, coarse="redundant", owned=None): """Reconfigure ``pc`` as a fresh PCMG (FMG F-cycle) driven by the supplied reduced->reduced prolongations ``Ps``, Galerkin RAP for coarse operators. @@ -603,7 +603,10 @@ def _configure_pcmg(pc, Ps, coarse="redundant"): the same bundle the native (DMPlex-refinement) route applies — so custom-P and native multigrid cannot be configured differently (#468). ``coarse`` selects the coarse-solve variant: ``"svd"`` on the rotated path, whose - Galerkin-coarsened operator inherits the rigid-rotation null space. + Galerkin-coarsened operator inherits the rigid-rotation null space. ``owned`` + is the solver's record of the option values UW3 has written, so a key the USER + set is left alone (see :meth:`multigrid_options.MGBundle.apply`); ``None`` writes + unconditionally, which is right for the rotated path's per-solve prefix. The bundle is written into the options DB under the PC's OWN options prefix (``pc.getOptionsPrefix()`` — e.g. ``Solver_N_`` for the scalar top-level PC, @@ -617,7 +620,7 @@ def _configure_pcmg(pc, Ps, coarse="redundant"): nlev = len(Ps) + 1 prefix = pc.getOptionsPrefix() or "" multigrid_options.geometric_mg_bundle(coarse=coarse).apply( - PETSc.Options(), prefix) + PETSc.Options(), prefix, owned=owned) pc.setType("mg") pc.setMGLevels(nlev) pc.setMGType(PETSc.PC.MGType.FULL) @@ -650,7 +653,8 @@ def _install_transfers(solver, Ps, verbose=False): ksp = solver.snes.getKSP() solver.snes.setUp() ksp.setDMActive(PETSc.KSP.DMActive.OPERATOR, False) - _configure_pcmg(ksp.getPC(), Ps) + _configure_pcmg(ksp.getPC(), Ps, + owned=solver._managed_pc_options) if verbose: from underworld3 import mpi mpi.pprint(f"[{solver.name}] custom FMG installed: {nlev} levels, " @@ -711,7 +715,7 @@ def _install_velocity_block_transfers(solver, Ps, verbose=False): vel_ksp.setDMActive(PETSc.KSP.DMActive.OPERATOR, False) vel_pc.reset() vel_pc.setOperators(A_vv, P_vv) - _configure_pcmg(vel_pc, Ps) + _configure_pcmg(vel_pc, Ps, owned=solver._managed_pc_options) vel_pc.setUp() # 4. re-attach the coupled Stokes nullspace (operator state was touched) diff --git a/src/underworld3/utilities/multigrid_options.py b/src/underworld3/utilities/multigrid_options.py index cd2e9d27..6ecfa33e 100644 --- a/src/underworld3/utilities/multigrid_options.py +++ b/src/underworld3/utilities/multigrid_options.py @@ -49,10 +49,32 @@ from typing import NamedTuple -__all__ = ["MGBundle", "geometric_mg_bundle", "gamg_bundle", +__all__ = ["MGBundle", "geometric_mg_bundle", "gamg_bundle", "option_string", "GEOMETRIC_MG_COARSE_SOLVERS"] +def option_string(value): + """How PETSc stores ``value`` in its options database, so a value UW3 wrote can + be compared against what is in there now. ``None`` is a bare flag (reads back as + an empty string); booleans lower-case.""" + if value is None: + return "" + if isinstance(value, bool): + return str(value).lower() + return str(value) + + +def _user_owns(opts, name, owned): + """Was ``name`` set by someone other than UW3's own option writers? + + True when the key is present and its current value is not the one UW3 last + recorded writing. ``owned is None`` disables the check (the caller owns the + prefix outright).""" + if owned is None or not opts.hasName(name): + return False + return opts.getString(name) != owned.get(name) + + #: Coarse-solve variants of the geometric bundle. ``"redundant"`` (redundant+LU) #: is the default and is np-safe; ``"svd"`` is required whenever the coarse #: operator inherits a null space — see :func:`geometric_mg_bundle`. @@ -71,17 +93,45 @@ class MGBundle(NamedTuple): settings: dict stale: tuple - def apply(self, opts, prefix=""): + def apply(self, opts, prefix="", owned=None): """Write this bundle into ``opts`` (a ``PETSc.Options``) under ``prefix``. ``prefix`` is the managed block's own prefix — ``""`` for a top-level PC, ``"fieldsplit_velocity_"`` for the Stokes velocity sub-block — and is applied on top of whatever prefix ``opts`` itself carries. + + Parameters + ---------- + owned : dict or None + Bookkeeping of the option values UW3 itself has written, keyed by full + option name. A key that is already present but whose current value is + **not** what UW3 last wrote was set by the user, and is left alone — + both here and in the stale-key clear, because it is not ours to remove. + Pass ``None`` to write unconditionally, which is right when the caller + owns the whole prefix (the rotated path builds a fresh, per-solve + prefix no user can address). + + Ownership is *recorded*, never inferred from the value. Guessing by + value fails as soon as a second internal writer touches the same key — + which is exactly how the ``tolerance`` and ``strategy`` setters defeated + an earlier attempt (#477). Every internal writer of a bundle key goes + through this bookkeeping (``SolverBaseClass._push_managed_option``), so + anything unrecorded is the user's by construction. """ for key, value in self.settings.items(): - opts.setValue(prefix + key, value) + name = prefix + key + if _user_owns(opts, name, owned): + continue + opts.setValue(name, value) + if owned is not None: + owned[name] = option_string(value) for key in self.stale: - opts.delValue(prefix + key) + name = prefix + key + if _user_owns(opts, name, owned): + continue + opts.delValue(name) + if owned is not None: + owned.pop(name, None) def _geometric_mg_settings(coarse): diff --git a/tests/test_1021_mg_option_bundle.py b/tests/test_1021_mg_option_bundle.py index 53b8983d..c99c09a3 100644 --- a/tests/test_1021_mg_option_bundle.py +++ b/tests/test_1021_mg_option_bundle.py @@ -201,6 +201,62 @@ def test_rotated_picks_up_a_mesh_owned_hierarchy(): assert s._rotated_freeslip_info["velocity_pc_type"] == "mg" +def _velocity_mg_config(s): + """(smoother ksp, smoother max_it, coarse pc) off the live velocity sub-PC.""" + pc = s.snes.getKSP().getPC() + pc.setUp() # MUST precede getFieldSplitSubKSP + vpc = pc.getFieldSplitSubKSP()[0].getPC() + assert vpc.getType() == "mg" and vpc.getMGLevels() > 1, ( + f"expected geometric MG on the velocity block, got {vpc.getType()}") + smoother = vpc.getMGSmoother(vpc.getMGLevels() - 1) + return (smoother.getType(), smoother.getTolerances()[3], + vpc.getMGCoarseSolve().getPC().getType()) + + +@pytest.mark.parametrize("preconditioner", [None, "fmg"]) +def test_user_set_bundle_keys_are_honoured(preconditioner): + """A user who sets a bundle key must get it, and must keep the managed values + for the keys they left alone — per key, not all-or-nothing. + + The bundle used to be applied wholesale on every rebuild. The only escape was + the ``_pc_user_override`` latch, which keys on ``pc_type`` ALONE, so setting + ``mg_levels_ksp_max_it`` was silently discarded unless the user also set + ``pc_type`` to the value it already had. Explicit ``preconditioner="fmg"`` was + worse — it skips that latch entirely, so the clearer the request the less + control it carried, which is why both modes are parametrised here. + """ + mesh = uw.meshing.Annulus(radiusInner=R_IN, radiusOuter=R_OUT, + cellSize=2 * RES, qdegree=3, refinement=1) + s = _stokes(mesh, f"u{'f' if preconditioner else 'a'}", rotated=False) + if preconditioner is not None: + s.preconditioner = preconditioner + # set TWO of the bundle's keys and leave the rest to the framework + s.petsc_options["fieldsplit_velocity_mg_levels_ksp_max_it"] = 6 + s.petsc_options["fieldsplit_velocity_mg_coarse_pc_type"] = "svd" + s.solve() + s.solve() # a rebuild must not clobber it + + smoother, max_it, coarse = _velocity_mg_config(s) + assert max_it == 6, "user-set smoother iteration count was overwritten" + assert coarse == "svd", "user-set coarse solver was overwritten" + # ...and the key the user did NOT set still carries the measured default + assert smoother == "gmres", ( + "respecting one bundle key must not abandon the rest of the bundle") + + +def test_unset_bundle_keys_keep_the_managed_defaults(): + """The control for the test above: with nothing set, the whole managed bundle + applies. Without this, 'user-set keys are honoured' could pass trivially by + never applying the bundle at all — which is exactly how the first attempt at + this failed (the routes silently reverted to GAMG).""" + mesh = uw.meshing.Annulus(radiusInner=R_IN, radiusOuter=R_OUT, + cellSize=2 * RES, qdegree=3, refinement=1) + s = _stokes(mesh, "dflt", rotated=False) + s.solve() + s.solve() + assert _velocity_mg_config(s) == ("gmres", 4, "redundant") + + def test_rotated_fmg_survives_repeated_newton_increments(): """The rotated path applies its bundle under a per-solve options prefix and then drops the keys again, so the global database stays bounded under From b4159f5125a0c54a695e3a6e70e7e0239bfeda76 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Thu, 30 Jul 2026 11:59:36 +1000 Subject: [PATCH 4/5] Fill solver.strategy from the measured smoother regimes; rename MGBundle to MGSettings MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `solver.strategy` has been a public property accepting "default", "robust" and "fast" since it was written, validating your input and then configuring all three identically — its docstring said "(Reserved)". A property that checks a value and then ignores it is the same defect class as #477 and #478: the failure is invisible because the solve still converges. This fills it, and it turns out the machinery this PR built is exactly what was missing underneath. TWO LAYERS, which is the design point: solver.strategy the named intent public, already existed MGSettings the option values utilities/multigrid_options.py solver.petsc_options the escape hatch public, outranks both (F-1) `geometric_mg_bundle(smoother=...)` gains the variant axis, from the two regimes measured in this PR's own work — neither dominates: "robust" = gmres/4 survives an operator a stationary smoother stalls on. Spiegelman notch (eta contrast 1e26, 4 levels) per-V-cycle contraction 0.56 vs richardson's 0.75, margin growing with depth; TI rotated annulus 11 velocity its -> 5, 1.74x. "fast" = richardson/3 cheaper per cycle where the operator is benign. Linear, SYMMETRIC annulus at eta contrast 1e6: beats robust on wall clock at every depth (x1.16, x1.30, x1.82 at 2, 3, 4 levels) while taking more iterations. "default" stays "robust", so nobody's results move. "fast" is the documented opt-out for the cost the adversarial review on this PR quoted and could not otherwise offer a way out of. ALSO REMOVES DEAD CODE THE MEASUREMENT EXPOSED. The `strategy` setter wrote the whole GAMG velocity bundle plus `pc_mg_type=kaskade`, the latter carrying a comment warning against changing it without benchmarking. Every one of those writes was DEAD: `_apply_preconditioner_options` runs later, at `_build`, and is the single writer of that block — it overwrites them with the geometric bundle on a refined mesh, or the GAMG bundle ("additive") without one. Measured in both orders and with `preconditioner="fmg"` set before and after: the live PC was mg/FULL every time, so `kaskade` never once took effect. The strategy's effect on the velocity block now runs through `_mg_smoother_variant`, which selects a bundle variant rather than racing the bundle writer. `MGBundle` -> `MGSettings` (internal, no call site outside the module and custom_mg): "MG" is what the module actually holds, and "Bundle" said nothing. Tests: strategy selects a real variant; "default" reproduces the framework default exactly (the control — without it, filling the axis could silently move results); a user-set option still beats the strategy while the strategy keeps the keys the user left alone (the two layers compose in the right order). Validation: three routes still agree bar the deliberate rotated svd; 79 passed across the MG, adaptivity and solve-report suites; `level_1 and tier_a` 526 passed with the two pre-existing #470 failures; parallel np2 and np4 green (6 / 9 / 2). Underworld development team with AI support from Claude Code --- docs/developer/subsystems/solvers.md | 43 ++++++--- .../cython/petsc_generic_snes_solvers.pyx | 85 +++++++++++----- src/underworld3/utilities/custom_mg.py | 2 +- .../utilities/multigrid_options.py | 96 ++++++++++++------- tests/test_1021_mg_option_bundle.py | 45 +++++++++ 5 files changed, 201 insertions(+), 70 deletions(-) diff --git a/docs/developer/subsystems/solvers.md b/docs/developer/subsystems/solvers.md index ba2e086b..eda4c33a 100644 --- a/docs/developer/subsystems/solvers.md +++ b/docs/developer/subsystems/solvers.md @@ -418,26 +418,47 @@ stationary iteration stalls where a Krylov smoother adapts its polynomial. It predicts the notch result, and it predicts the cost below on well-conditioned problems. -### The smoother is not free +### Two regimes, and the named intent that selects them -On a *well-conditioned* problem gmres/4 is slower than richardson/3 despite -converging in fewer iterations. Nested annulus, linear (symmetric) velocity -block, η contrast 1e6, outer KSP timed in isolation: +Neither smoother dominates, so there are two variants and `solver.strategy` +chooses between them. This is the layering: + +| layer | what it is | where | +|---|---|---| +| `solver.strategy` | the **named intent** — `"default"` / `"robust"` / `"fast"` | public property on the solver | +| `MGSettings` | the **option values** an intent resolves to | `utilities/multigrid_options.py` | +| `solver.petsc_options` | the **escape hatch** — any key you set here outranks both | public, unchanged | + +`"robust"` (gmres/4) survives an operator a stationary smoother stalls on: the +Spiegelman notch contraction above, and 11 → 5 velocity iterations on the +transversely isotropic rotated annulus. + +`"fast"` (richardson/3) is cheaper per cycle and quicker where the operator is +benign. Nested annulus, **linear (symmetric)** velocity block, η contrast 1e6, +outer KSP timed in isolation — `"fast"` against `"robust"`: | depth | iterations | wall clock | |---|---|---| -| 2 levels | ×1.45 | **×0.86** | -| 3 levels | ×1.60 | **×0.77** | -| 4 levels | ×1.20 | **×0.55** | +| 2 levels | ×0.69 | **×1.16** | +| 3 levels | ×0.63 | **×1.30** | +| 4 levels | ×0.83 | **×1.82** | The extra sweep costs ~3%; gmres itself costs ~40% per cycle. Custom-P builds its Galerkin coarse operators from barycentric transfers, denser than the native nested ones, so the smoother is a larger share of the cycle on that route. -It stays the default on both routes because the failure it avoids (a stationary -smoother stalling on a high-contrast operator) is worse than the cost it carries, -and it carries that cost exactly where the problems are easy. Removing a guardrail -is the caller's decision; the number is recorded so the decision can be made. +`"default"` is `"robust"`: the failure it avoids is worse than the cost it carries, +and it carries that cost exactly where the problem is easy. `"fast"` is the +documented opt-out. + +```{warning} +Do not fill a strategy by writing velocity-block options from the strategy setter. +`_apply_preconditioner_options` runs later, at `_build`, and is the single writer +of that block — anything written earlier is overwritten. That is how +`pc_mg_type=kaskade` sat in the `strategy` setter for a long time carrying a +comment warning against changing it, while never once taking effect (measured: the +live PC was `mg`/FULL in every ordering). A strategy selects a bundle *variant*. +``` ### The one legitimate per-route difference diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index ced17dd6..3abc62af 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -175,6 +175,16 @@ class SolverBaseClass(uw_object): # utilities.custom_mg). None => standard FMG/GAMG path, unchanged. self._custom_mg = None + @property + def _mg_smoother_variant(self): + """Which measured smoother regime this solver's strategy asks for. + + ``solver.strategy`` is the named intent ("I want speed" / "I want this to + converge"); the values live in ``utilities.multigrid_options``. Solvers with + no strategy axis get the robust default. See + :func:`multigrid_options.geometric_mg_bundle` for the measurements.""" + return "fast" if getattr(self, "_strategy", "default") == "fast" else "robust" + def _push_managed_option(self, key, value): """Write a PETSc option UW3 owns, recording that we wrote it. @@ -660,9 +670,10 @@ class SolverBaseClass(uw_object): # `preconditioner="fmg"` was worse: it skips the latch entirely, so the # clearer the request the less control it carried. if want_fmg: - multigrid_options.geometric_mg_bundle().apply( - PETSc.Options(), self.petsc_options_prefix + prefix, - owned=self._managed_pc_options) + multigrid_options.geometric_mg_bundle( + smoother=self._mg_smoother_variant).apply( + PETSc.Options(), self.petsc_options_prefix + prefix, + owned=self._managed_pc_options) self._pc_managed_value = "mg" else: if self._preconditioner == "fmg" and n_levels <= 1 and uw.mpi.rank == 0: @@ -6309,19 +6320,39 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): @property def strategy(self): """ - Solver strategy controlling preconditioner configuration. - - Currently supports: - - ``"default"``: Standard Schur complement fieldsplit with GAMG - - ``"robust"``: (Reserved) More robust but slower configuration - - ``"fast"``: (Reserved) Faster but less robust configuration - - Setting this property reconfigures the entire preconditioner stack. + What this solve should optimise for — the named intent over the + multigrid smoother's two measured regimes. + + - ``"default"``, ``"robust"``: ``gmres``/4 smoothing. Survives an operator a + stationary smoother stalls on: Spiegelman notch (:math:`\eta` contrast + 1e26, 4 levels) per-V-cycle contraction 0.56 against richardson's 0.75, the + margin growing with depth; transversely isotropic rotated annulus 11 + velocity iterations down to 5. + - ``"fast"``: ``richardson``/3 smoothing. Cheaper per cycle and quicker where + the operator is benign — on a linear, symmetric annulus at + :math:`\eta` contrast 1e6 it beats ``"robust"`` on wall clock at every + hierarchy depth tested (x1.16, x1.30, x1.82 at 2, 3, 4 levels) while taking + more iterations. It gives up the regime ``"robust"`` exists for, so it is + an opt-in. + + ``"default"`` is ``"robust"``: the failure it avoids is worse than the cost it + carries, and it carries that cost exactly where the problem is easy. + + Setting this property also resets the fieldsplit / Schur / pressure sub-solve + configuration to the framework defaults. It does **not** write the velocity + block's preconditioner directly — that is applied later, from + :mod:`underworld3.utilities.multigrid_options`, which is the single writer of + those options; this property selects which variant it applies. Values written + by hand into :attr:`petsc_options` are respected and not overwritten. Returns ------- str Current strategy name. + + See Also + -------- + preconditioner : which multigrid FAMILY to use (geometric, algebraic, auto). """ return self._strategy @@ -6332,14 +6363,19 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): f"Unknown solver strategy {value!r}: " "expected 'default', 'robust', or 'fast'." ) - # 'robust' and 'fast' are accepted names but currently configure the - # same option bundle as 'default' (their dedicated branches were empty - # placeholders and have been removed — Charter S5, READ-23). + # 'fast' and 'robust' now select a real smoother variant, via + # `_mg_smoother_variant` -> `multigrid_options.geometric_mg_bundle`. They were + # accepted-and-inert placeholders for a long time: validated on input, then + # configured identically to 'default'. A property that checks your value and + # then ignores it is the same defect class as #477 and #478 — the failure is + # invisible, because the solve still converges. # self.is_setup = False self._strategy = value - # All strategies: reset to preferred + # Common to every strategy: reset the fieldsplit / Schur / pressure + # sub-solve to the framework defaults. The strategy's effect on the VELOCITY + # BLOCK is carried by `_mg_smoother_variant`, not by writes from here. self.petsc_options["snes_ksp_ew"] = None self.petsc_options["snes_ksp_ew_version"] = 3 @@ -6376,16 +6412,15 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): # preconditioning and weakly-indefinite coarse operators; issue #147). self.petsc_options[f"fieldsplit_velocity_ksp_type"] = "fgmres" self.petsc_options[f"fieldsplit_velocity_ksp_max_it"] = 200 - self._push_managed_option(f"fieldsplit_velocity_pc_type", "gamg") - self._push_managed_option(f"fieldsplit_velocity_pc_gamg_type", "agg") - self._push_managed_option(f"fieldsplit_velocity_pc_gamg_repartition", True) - # NOTE: "kaskade" here diverges from the "additive" set by __init__'s - # velocity block — a long-standing (possibly unintentional) difference. - # Do not change without benchmarking (READ-23 kept the value as-is). - self._push_managed_option(f"fieldsplit_velocity_pc_mg_type", "kaskade") - self._push_managed_option(f"fieldsplit_velocity_pc_gamg_agg_nsmooths", 2) - self._push_managed_option(f"fieldsplit_velocity_mg_levels_ksp_max_it", 3) - self._push_managed_option(f"fieldsplit_velocity_mg_levels_ksp_converged_maxits", None) + # The velocity BLOCK's preconditioner is not set here. `strategy` used to + # write the whole GAMG bundle plus `pc_mg_type=kaskade`, and every one of + # those writes was DEAD: `_apply_preconditioner_options` runs later (at + # `_build`) and overwrites them with the geometric bundle on a refined mesh, + # or the GAMG bundle (`additive`) without one. Measured both orders — the + # live PC was `mg`/FULL every time, so `kaskade` never once took effect + # despite the comment warning against changing it. The strategy's effect on + # the velocity block now runs through `_mg_smoother_variant`, which selects + # a bundle variant instead of racing the bundle writer. diff --git a/src/underworld3/utilities/custom_mg.py b/src/underworld3/utilities/custom_mg.py index 50c8504f..588b0031 100644 --- a/src/underworld3/utilities/custom_mg.py +++ b/src/underworld3/utilities/custom_mg.py @@ -605,7 +605,7 @@ def _configure_pcmg(pc, Ps, coarse="redundant", owned=None): selects the coarse-solve variant: ``"svd"`` on the rotated path, whose Galerkin-coarsened operator inherits the rigid-rotation null space. ``owned`` is the solver's record of the option values UW3 has written, so a key the USER - set is left alone (see :meth:`multigrid_options.MGBundle.apply`); ``None`` writes + set is left alone (see :meth:`multigrid_options.MGSettings.apply`); ``None`` writes unconditionally, which is right for the rotated path's per-solve prefix. The bundle is written into the options DB under the PC's OWN options prefix diff --git a/src/underworld3/utilities/multigrid_options.py b/src/underworld3/utilities/multigrid_options.py index 6ecfa33e..4f2700a5 100644 --- a/src/underworld3/utilities/multigrid_options.py +++ b/src/underworld3/utilities/multigrid_options.py @@ -38,7 +38,9 @@ Apply the geometric bundle to a solver's managed block:: from underworld3.utilities import multigrid_options - multigrid_options.geometric_mg_bundle().apply(self.petsc_options, prefix) + multigrid_options.geometric_mg_bundle().apply( + PETSc.Options(), self.petsc_options_prefix + prefix, + owned=self._managed_pc_options) Read the settings without writing them (the rotated path stages its options in a dict so it can drop them again after ``setUp``):: @@ -49,8 +51,8 @@ from typing import NamedTuple -__all__ = ["MGBundle", "geometric_mg_bundle", "gamg_bundle", "option_string", - "GEOMETRIC_MG_COARSE_SOLVERS"] +__all__ = ["MGSettings", "geometric_mg_bundle", "gamg_bundle", "option_string", + "GEOMETRIC_MG_COARSE_SOLVERS", "GEOMETRIC_MG_SMOOTHERS"] def option_string(value): @@ -80,8 +82,14 @@ def _user_owns(opts, name, owned): #: operator inherits a null space — see :func:`geometric_mg_bundle`. GEOMETRIC_MG_COARSE_SOLVERS = ("redundant", "svd") +#: Smoother variants. These are the two measured regimes, not a taste setting: +#: ``"robust"`` survives a badly-conditioned operator where a stationary smoother +#: stalls, ``"fast"`` is cheaper per cycle where the operator is benign. Selected by +#: ``solver.strategy``; see :func:`geometric_mg_bundle` for the numbers. +GEOMETRIC_MG_SMOOTHERS = ("robust", "fast") -class MGBundle(NamedTuple): + +class MGSettings(NamedTuple): """A preconditioner option bundle: the keys to set, and the keys to clear. ``settings`` maps an option suffix (no prefix, no leading ``-``) to its @@ -134,8 +142,11 @@ def apply(self, opts, prefix="", owned=None): owned.pop(name, None) -def _geometric_mg_settings(coarse): - """The geometric-MG settings for one coarse-solve variant.""" +def _geometric_mg_settings(coarse, smoother="robust"): + """The geometric-MG settings for one coarse-solve and smoother variant. + + Both variants set the SAME keys — only values differ — so the derived stale-key + sets are variant-independent.""" settings = { "pc_type": "mg", "pc_mg_type": "full", # FMG (F-cycle) @@ -173,12 +184,7 @@ def _geometric_mg_settings(coarse): # spread far enough that a stationary iteration stalls where a Krylov # smoother adapts its polynomial. That also predicts the cost measured on # well-conditioned problems (see the note on geometric_mg_bundle). - "mg_levels_ksp_type": "gmres", "mg_levels_pc_type": "sor", - # SET the count, never inherit it. PCMG's own default is 2 and the GAMG - # bundle below leaves 3 under the same prefix, so a bundle that omits - # this key smooths differently depending on what ran before it (#468). - "mg_levels_ksp_max_it": 4, # Run EXACTLY max_it smoother iterations: no residual-norm computation and # no convergence test, so every V-cycle costs the same. A Krylov smoother # makes the cycle non-stationary, which is why the velocity block is @@ -186,6 +192,20 @@ def _geometric_mg_settings(coarse): "mg_levels_ksp_norm_type": "none", "mg_levels_ksp_converged_maxits": None, } + if smoother == "robust": + settings["mg_levels_ksp_type"] = "gmres" + settings["mg_levels_ksp_max_it"] = 4 + elif smoother == "fast": + # Cheaper per cycle and measurably faster where the operator is benign, + # at the cost of the regime "robust" exists for. See the docstring. + settings["mg_levels_ksp_type"] = "richardson" + settings["mg_levels_ksp_max_it"] = 3 + else: + raise ValueError( + f"smoother must be one of {GEOMETRIC_MG_SMOOTHERS} (got {smoother!r})") + # SET the iteration count, never inherit it. PCMG's own default is 2 and the + # GAMG bundle leaves 3 under the same prefix, so a bundle that omits this key + # smooths differently depending on what ran before it (#468). if coarse == "redundant": # 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 @@ -231,15 +251,19 @@ def _all_keys(): def _bundle(settings): - return MGBundle(settings=settings, + return MGSettings(settings=settings, stale=tuple(sorted(_all_keys() - set(settings)))) -def geometric_mg_bundle(coarse="redundant"): +def geometric_mg_bundle(coarse="redundant", smoother="robust"): """Geometric multigrid (FMG F-cycle) on an explicit hierarchy. Parameters ---------- + smoother : {"robust", "fast"} + Which of the two MEASURED regimes to configure. ``"robust"`` is + ``gmres``/4 and is the default; ``"fast"`` is ``richardson``/3. Chosen by + ``solver.strategy``, not usually here. coarse : {"redundant", "svd"} The coarse-level solve. ``"redundant"`` (redundant+LU) is the default. Use ``"svd"`` when the coarse operator inherits a null space — which the @@ -248,28 +272,34 @@ def geometric_mg_bundle(coarse="redundant"): Returns ------- - MGBundle - Settings and stale keys; see :meth:`MGBundle.apply`. + MGSettings + Settings and stale keys; see :meth:`MGSettings.apply`. Notes ----- - The ``gmres``/4 smoother is not free on a well-conditioned problem. Measured on - a nested annulus with a linear (symmetric) velocity block at :math:`\eta` - contrast 1e6, outer KSP timed in isolation, ``gmres``/4 against - ``richardson``/3: it wins on iterations at every depth (×1.45, ×1.60, ×1.20 at - 2, 3, 4 levels) and **loses on wall clock at every depth** (×0.86, ×0.77, - ×0.55) — the extra sweep costs ~3% and ``gmres`` itself ~40% per cycle. Custom-P - builds its Galerkin coarse operators from barycentric transfers, denser than - the native nested ones, so the smoother is a larger share of the cycle here. - - It is the default anyway, on both routes, because the failure it avoids is - worse than the cost it carries: a stationary smoother stalls where a Krylov one - adapts, and it does so exactly on the badly-conditioned high-contrast problems - this code exists for (the Spiegelman-notch measurement above). Removing a - guardrail is the caller's decision to take; the number is quoted so that - decision can be made. + The two smoother variants are two measured regimes, and neither dominates. + + ``"robust"`` (``gmres``/4) survives an operator a stationary smoother stalls on. + Spiegelman notch, :math:`\eta` contrast 1e26, nested 4-level hierarchy: + per-V-cycle contraction :math:`\rho` = 0.56 against richardson's 0.75 at the same + four iterations, the margin growing with depth. Transversely isotropic rotated + annulus: 11 velocity iterations down to 5, and 1.74x on the linear solve. + + ``"fast"`` (``richardson``/3) is cheaper per cycle and measurably quicker where + the operator is benign. Nested annulus, LINEAR (symmetric) velocity block at + :math:`\eta` contrast 1e6, outer KSP timed in isolation: robust wins on + iterations at every depth (x1.45, x1.60, x1.20 at 2, 3, 4 levels) and **loses on + wall clock at every depth** (x0.86, x0.77, x0.55). The extra sweep costs ~3% and + ``gmres`` itself ~40% per cycle; custom-P builds its Galerkin coarse operators + from barycentric transfers, denser than the native nested ones, so the smoother + is a larger share of the cycle on that route. + + ``"robust"`` is the default because the failure it avoids is worse than the cost + it carries, and it carries that cost exactly where the problem is easy. Choose + between them with ``solver.strategy``, which is the named-intent layer over + these values. """ - return _bundle(_geometric_mg_settings(coarse)) + return _bundle(_geometric_mg_settings(coarse, smoother)) def gamg_bundle(): @@ -277,7 +307,7 @@ def gamg_bundle(): Returns ------- - MGBundle - Settings and stale keys; see :meth:`MGBundle.apply`. + MGSettings + Settings and stale keys; see :meth:`MGSettings.apply`. """ return _bundle(_gamg_settings()) diff --git a/tests/test_1021_mg_option_bundle.py b/tests/test_1021_mg_option_bundle.py index c99c09a3..9934e252 100644 --- a/tests/test_1021_mg_option_bundle.py +++ b/tests/test_1021_mg_option_bundle.py @@ -257,6 +257,51 @@ def test_unset_bundle_keys_keep_the_managed_defaults(): assert _velocity_mg_config(s) == ("gmres", 4, "redundant") +def _strategy_mg_config(strategy): + mesh = uw.meshing.Annulus(radiusInner=R_IN, radiusOuter=R_OUT, + cellSize=2 * RES, qdegree=3, refinement=1) + s = _stokes(mesh, f"s{strategy or 'none'}"[:6], rotated=False) + if strategy is not None: + s.strategy = strategy + s.solve() + s.solve() + return _velocity_mg_config(s) + + +def test_strategy_selects_a_real_smoother_variant(): + """``solver.strategy`` must actually change the smoother. + + ``"fast"`` and ``"robust"`` were accepted-and-inert for a long time: validated + on input, then configured identically to ``"default"``. A property that checks + your value and then ignores it is the same defect class as #477/#478 — invisible, + because the solve still converges. They now select a measured variant: + ``"robust"`` is gmres/4, ``"fast"`` is richardson/3. + """ + assert _strategy_mg_config("fast")[:2] == ("richardson", 3) + assert _strategy_mg_config("robust")[:2] == ("gmres", 4) + + +def test_default_strategy_does_not_change_behaviour(): + """The control: ``"default"`` must reproduce the framework default exactly, so + filling the strategy axis moves nobody's results.""" + assert _strategy_mg_config("default") == _strategy_mg_config(None) + + +def test_a_user_key_still_beats_the_strategy(): + """The two layers compose in the right order: an explicit option the user wrote + outranks the strategy's choice of variant.""" + mesh = uw.meshing.Annulus(radiusInner=R_IN, radiusOuter=R_OUT, + cellSize=2 * RES, qdegree=3, refinement=1) + s = _stokes(mesh, "sxu", rotated=False) + s.strategy = "fast" # would ask for richardson/3 + s.petsc_options["fieldsplit_velocity_mg_levels_ksp_type"] = "chebyshev" + s.solve() + s.solve() + smoother, max_it, _ = _velocity_mg_config(s) + assert smoother == "chebyshev", "the user's smoother lost to the strategy" + assert max_it == 3, "the strategy should still own the keys the user left alone" + + def test_rotated_fmg_survives_repeated_newton_increments(): """The rotated path applies its bundle under a per-solve options prefix and then drops the keys again, so the global database stays bounded under From 13bd82ceb6ec05952d544b16bde2fb0dbc6a0a24 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Thu, 30 Jul 2026 12:27:12 +1000 Subject: [PATCH 5/5] solver.strategy reports the preconditioner it resolved to (#484 shape) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Asking "what am I actually running?" required knowing which nine PETSc option keys to look up, and `_pc_managed_value` / `_pc_user_override` were private. That is the gap the configuration audit filed as #484 — of twelve preconditioner fallbacks, ten left no readable state — and `solver.strategy` is the natural place to close it for this block, since the strategy is what the user set. The getter now returns the strategy NAME that also reports what it resolved to: >>> stokes.strategy 'default' — geometric multigrid (2 levels), full cycle, smoother gmresx4 + sor, coarse redundant/lu >>> stokes.strategy == "default" True `_StrategyName` subclasses `str` deliberately, so comparisons, formatting and serialisation are unchanged — nothing that used `solver.strategy` before behaves differently. Verified nothing in src/ or tests/ compares it to a string in the first place; they only set it. Three properties of the report matter more than its wording: - It NAMES a user override rather than absorbing it, so the summary cannot hide the difference between what the strategy asked for and what is running. - It REFUSES to report before it knows. The preconditioner resolves at the first solve; until then the report says so instead of presenting the constructor defaults as the answer. The first version of this did present them, which is precisely the stale-but-authoritative summary this reporting exists to prevent — hence the new `_pc_resolved` flag. - `solver.preconditioner_settings` is the same information as a dict, so a test can assert on it instead of parsing prose or inferring from timings. `multigrid_options.describe()` does the formatting, next to the values it describes, so a new bundle variant cannot be added without the report following it. Validation: 14 passed in test_1021; 49 passed across test_1014/1017/1018/1020/1058; `level_1 and tier_a` 526 passed with the two pre-existing #470 failures; parallel np2 and np4 15 passed each. Underworld development team with AI support from Claude Code --- docs/developer/subsystems/solvers.md | 30 ++++++ .../cython/petsc_generic_snes_solvers.pyx | 98 ++++++++++++++++++- .../utilities/multigrid_options.py | 38 +++++++ tests/test_1021_mg_option_bundle.py | 42 ++++++++ 4 files changed, 207 insertions(+), 1 deletion(-) diff --git a/docs/developer/subsystems/solvers.md b/docs/developer/subsystems/solvers.md index eda4c33a..e03c9309 100644 --- a/docs/developer/subsystems/solvers.md +++ b/docs/developer/subsystems/solvers.md @@ -451,6 +451,36 @@ nested ones, so the smoother is a larger share of the cycle on that route. and it carries that cost exactly where the problem is easy. `"fast"` is the documented opt-out. +### Asking what you actually got + +`solver.strategy` reports it. The value is still the strategy *name* — it compares +and formats as the plain string, so nothing that used it before changes — but +displaying it shows the preconditioner that name resolved to: + +```python +>>> stokes.strategy +'default' — geometric multigrid (2 levels), full cycle, smoother gmresx4 + sor, + coarse redundant/lu +>>> stokes.strategy == "default" +True +``` + +Three properties of that report matter more than its formatting: + +- **It names a user override rather than absorbing it.** A key you set yourself is + listed as such (`; overridden by the user: mg_levels_ksp_max_it=6`), so the + summary cannot hide the difference between what the strategy asked for and what + is running. +- **It refuses to report before it knows.** The preconditioner is resolved at the + first solve; until then the report says so instead of presenting the constructor + defaults as though they were the answer. A summary that looks authoritative and + is stale is worse than no summary. +- **`solver.preconditioner_settings`** is the same information as a dict, so a test + can assert on it rather than parsing prose or inferring from timings. + +This is the general shape #484 asks for across all the managed fallbacks: report the +resolved choice as readable state, and distinguish "could not" from "chose not to". + ```{warning} Do not fill a strategy by writing velocity-block options from the strategy setter. `_apply_preconditioner_options` runs later, at `_build`, and is the single writer diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index 3abc62af..715c78f1 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -14,6 +14,31 @@ import underworld3.timing as timing from underworld3.utilities._api_tools import uw_object from underworld3.utilities import multigrid_options + +class _StrategyName(str): + """A strategy name that also reports what it resolved to. + + Subclasses ``str`` deliberately: ``solver.strategy == "fast"``, string + formatting and serialisation all behave exactly as before, but displaying it — + in a REPL, a notebook, or a log line — shows the preconditioner it actually + configured. Asking "what am I running?" should not require knowing which nine + PETSc option keys to look up. + + Use :attr:`SolverBaseClass.preconditioner_settings` for the machine-readable + form. + """ + + def __new__(cls, name, summary): + obj = super().__new__(cls, name) + obj._summary = summary + return obj + + def __repr__(self): + return f"{str.__repr__(self)} — {self._summary}" + + def _repr_markdown_(self): + return f"**`{str(self)}`** — {self._summary}" + from underworld3.function import expression as public_expression expression = lambda *x, **X: public_expression(*x, _unique_name_generation=True, **X) @@ -170,11 +195,59 @@ class SolverBaseClass(uw_object): # the `tolerance` and `strategy` setters defeated an earlier attempt (#477). # Internal writers therefore go through _push_managed_option(). self._managed_pc_options = {} + # Has _apply_preconditioner_options actually made the resolution decision + # yet? Until it has, the options database holds only this solver's __init__ + # defaults, which is NOT what the next solve will run — reporting them as + # resolved would be exactly the stale-but-authoritative-looking summary this + # reporting exists to prevent. + self._pc_resolved = False # Custom multigrid prolongation hierarchy (see set_custom_mg / # utilities.custom_mg). None => standard FMG/GAMG path, unchanged. self._custom_mg = None + @property + def preconditioner_settings(self): + """The option values the managed preconditioner block is configured with. + + A read-only dict of the multigrid keys UW3 currently has in the options + database for this solver's managed block, so what actually got applied can + be asserted on instead of inferred from timings. Empty for a solver with no + managed block (``_pc_option_prefix is None``), or before the first build + resolves one. + + The values are what the strategy resolved to, *including* any key you set + yourself — those are respected (see :attr:`petsc_options`). Use + :attr:`strategy` for a readable summary of the same thing. + """ + prefix = self._pc_option_prefix + if prefix is None: + return {} + keys = set(multigrid_options.gamg_bundle().settings) + for coarse in multigrid_options.GEOMETRIC_MG_COARSE_SOLVERS: + keys |= set(multigrid_options.geometric_mg_bundle(coarse=coarse).settings) + out = {} + for key in sorted(keys): + name = prefix + key + if self.petsc_options.hasName(name): + out[key] = self.petsc_options.getString(name) + return out + + @property + def _user_overridden_pc_options(self): + """The managed-block keys the USER set, as (key, value) pairs. + + A key present in the options database whose value is not the one UW3 + recorded writing is theirs — the same test the bundle writer uses to decide + what to leave alone.""" + prefix = self._pc_option_prefix + if prefix is None: + return () + qualified = self.petsc_options_prefix + return tuple( + (key, value) for key, value in self.preconditioner_settings.items() + if self._managed_pc_options.get(qualified + prefix + key) != value) + @property def _mg_smoother_variant(self): """Which measured smoother regime this solver's strategy asks for. @@ -591,6 +664,9 @@ class SolverBaseClass(uw_object): prefix = self._pc_option_prefix if prefix is None: return + # Every path from here is a resolution decision, including "the user owns + # these options, leave them alone". + self._pc_resolved = True opts = self.petsc_options @@ -6350,11 +6426,31 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): str Current strategy name. + The value returned is the strategy name (it compares and formats as the + plain string) and additionally reports the preconditioner it resolved to + when displayed:: + + >>> stokes.strategy + 'default' — geometric multigrid (3 levels), full cycle, + smoother gmresx4 + sor, coarse redundant/lu + See Also -------- preconditioner : which multigrid FAMILY to use (geometric, algebraic, auto). + preconditioner_settings : the same information as a dict, for assertions. """ - return self._strategy + settings = self.preconditioner_settings + if not settings or not self._pc_resolved: + summary = "not resolved yet — configured at the first solve" + if settings: + summary += (f" (framework defaults in place: " + f"{multigrid_options.describe(settings)})") + else: + levels = len(getattr(self.mesh, "dm_hierarchy", []) or []) or None + summary = multigrid_options.describe( + settings, levels=levels, + overridden=self._user_overridden_pc_options) + return _StrategyName(self._strategy, summary) @strategy.setter def strategy(self, value): diff --git a/src/underworld3/utilities/multigrid_options.py b/src/underworld3/utilities/multigrid_options.py index 4f2700a5..2ba90a26 100644 --- a/src/underworld3/utilities/multigrid_options.py +++ b/src/underworld3/utilities/multigrid_options.py @@ -52,6 +52,7 @@ from typing import NamedTuple __all__ = ["MGSettings", "geometric_mg_bundle", "gamg_bundle", "option_string", + "describe", "GEOMETRIC_MG_COARSE_SOLVERS", "GEOMETRIC_MG_SMOOTHERS"] @@ -142,6 +143,43 @@ def apply(self, opts, prefix="", owned=None): owned.pop(name, None) +def describe(settings, levels=None, overridden=()): + """One line of plain English for a resolved bundle. + + Parameters + ---------- + settings : dict + The bundle's settings, as applied. + levels : int or None + Hierarchy depth, when known. + overridden : sequence of (key, value) + Keys the user has taken over, which the bundle left alone. + """ + pc = settings.get("pc_type", "?") + if pc == "mg": + what = "geometric multigrid" + if levels: + what += f" ({levels} levels)" + cycle = settings.get("pc_mg_type", "?") + smoother = (f"{settings.get('mg_levels_ksp_type', '?')}" + f"x{settings.get('mg_levels_ksp_max_it', '?')}" + f" + {settings.get('mg_levels_pc_type', '?')}") + coarse = settings.get("mg_coarse_pc_type", "?") + if coarse == "redundant": + coarse += f"/{settings.get('mg_coarse_redundant_pc_type', '?')}" + text = f"{what}, {cycle} cycle, smoother {smoother}, coarse {coarse}" + elif pc == "gamg": + text = (f"algebraic multigrid (gamg/{settings.get('pc_gamg_type', '?')}), " + f"{settings.get('pc_mg_type', '?')} cycle, " + f"smoother max_it {settings.get('mg_levels_ksp_max_it', '?')}") + else: + text = f"pc_type={pc}" + if overridden: + text += ("; overridden by the user: " + + ", ".join(f"{k}={v}" for k, v in overridden)) + return text + + def _geometric_mg_settings(coarse, smoother="robust"): """The geometric-MG settings for one coarse-solve and smoother variant. diff --git a/tests/test_1021_mg_option_bundle.py b/tests/test_1021_mg_option_bundle.py index 9934e252..7468388d 100644 --- a/tests/test_1021_mg_option_bundle.py +++ b/tests/test_1021_mg_option_bundle.py @@ -302,6 +302,48 @@ def test_a_user_key_still_beats_the_strategy(): assert max_it == 3, "the strategy should still own the keys the user left alone" +def test_strategy_reports_what_it_resolved_to(): + """``solver.strategy`` reports the preconditioner it configured, so "what am I + running?" does not require knowing which nine PETSc keys to look up (#484). + + It must still BE the string: comparisons, formatting and serialisation of a + strategy name are unchanged. + """ + mesh = uw.meshing.Annulus(radiusInner=R_IN, radiusOuter=R_OUT, + cellSize=2 * RES, qdegree=3, refinement=1) + s = _stokes(mesh, "rep", rotated=False) + + # before any solve: must NOT present the __init__ defaults as resolved + assert s.strategy == "default" + assert f"{s.strategy}" == "default" + assert "not resolved yet" in repr(s.strategy) + + s.solve() + assert s.strategy == "default" # still the plain string + summary = repr(s.strategy) + assert "geometric multigrid (2 levels)" in summary + assert "gmres" in summary and "sor" in summary + assert "not resolved yet" not in summary + + # and the machine-readable form agrees with it + assert s.preconditioner_settings["mg_levels_ksp_type"] == "gmres" + assert s.preconditioner_settings["mg_levels_ksp_max_it"] == "4" + + +def test_strategy_report_names_a_user_override(): + """A key the user took over must be called out, not silently folded into the + summary — otherwise the report reintroduces the ambiguity it exists to remove.""" + mesh = uw.meshing.Annulus(radiusInner=R_IN, radiusOuter=R_OUT, + cellSize=2 * RES, qdegree=3, refinement=1) + s = _stokes(mesh, "rov", rotated=False) + s.petsc_options["fieldsplit_velocity_mg_levels_ksp_max_it"] = 6 + s.solve() + summary = repr(s.strategy) + assert "overridden by the user" in summary + assert "mg_levels_ksp_max_it=6" in summary + assert s.preconditioner_settings["mg_levels_ksp_max_it"] == "6" + + def test_rotated_fmg_survives_repeated_newton_increments(): """The rotated path applies its bundle under a per-solve options prefix and then drops the keys again, so the global database stays bounded under