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..e03c9309 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,193 @@ 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. 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. + +### Two regimes, and the named intent that selects them + +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 | ×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. + +`"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. + +### 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 +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 + +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 +524,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..715c78f1 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -12,6 +12,32 @@ 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 + + +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) @@ -160,11 +186,95 @@ 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 = {} + # 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. + + ``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. + + 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"``. @@ -554,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 @@ -618,46 +731,25 @@ 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. + # `_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: - # 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( + 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: @@ -668,19 +760,9 @@ 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( + PETSc.Options(), self.petsc_options_prefix + prefix, + owned=self._managed_pc_options) self._pc_managed_value = "gamg" def _enforce_galerkin_for_geometric_mg(self): @@ -2987,16 +3069,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 @@ -3900,14 +3982,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 @@ -4912,14 +4994,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 @@ -5779,13 +5861,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 = {} @@ -6314,21 +6396,61 @@ 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. + + 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): @@ -6337,14 +6459,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 @@ -6381,16 +6508,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.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 - # 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 + # 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 d4de0371..588b0031 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,32 @@ 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", owned=None): """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. ``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.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 (``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, owned=owned) pc.setType("mg") pc.setMGLevels(nlev) pc.setMGType(PETSc.PC.MGType.FULL) @@ -649,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, " @@ -710,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) @@ -963,25 +968,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 +1020,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 +1060,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..2ba90a26 --- /dev/null +++ b/src/underworld3/utilities/multigrid_options.py @@ -0,0 +1,351 @@ +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( + 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``):: + + cfg.update({f"fieldsplit_vel_{k}": v + for k, v in multigrid_options.gamg_bundle().settings.items()}) +""" + +from typing import NamedTuple + +__all__ = ["MGSettings", "geometric_mg_bundle", "gamg_bundle", "option_string", + "describe", + "GEOMETRIC_MG_COARSE_SOLVERS", "GEOMETRIC_MG_SMOOTHERS"] + + +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`. +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 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 + 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="", 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(): + 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: + name = prefix + key + if _user_owns(opts, name, owned): + continue + opts.delValue(name) + if owned is not 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. + + 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) + # 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. + # + # 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, 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_pc_type": "sor", + # 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 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 + # 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 MGSettings(settings=settings, + stale=tuple(sorted(_all_keys() - set(settings)))) + + +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 + Galerkin-coarsened *rotated* velocity block always does, because the + rigid rotations survive the constraint. + + Returns + ------- + MGSettings + Settings and stale keys; see :meth:`MGSettings.apply`. + + Notes + ----- + 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, smoother)) + + +def gamg_bundle(): + """Algebraic multigrid — the fallback when no geometric hierarchy exists. + + Returns + ------- + MGSettings + Settings and stale keys; see :meth:`MGSettings.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..7468388d --- /dev/null +++ b/tests/test_1021_mg_option_bundle.py @@ -0,0 +1,383 @@ +"""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 _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 _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_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 + 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"