From ad9cead7d3f4bd84919447f17d9d84155ce67320 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Fri, 18 Sep 2026 16:24:37 +1000 Subject: [PATCH] Add a "bulletproof" smoother strategy: Schwarz, to take the solver off the suspect list A third variant on the geometric-MG smoother axis, selected by `solver.strategy = "bulletproof"`: gmres/8 preconditioned by additive Schwarz instead of SOR. Each subdomain solve inverts the local coupling directly, so the smoother does not rely on the operator being close to symmetric or elliptic. Its purpose is diagnostic as much as numerical. When a run misbehaves and the linear solve, the transport scheme and the mesh resolution are all candidates, rerunning under this variant removes the first from the list. Measured on creeping Oldroyd-B past a confined cylinder (Wi 0.4, beta 0.59, three levels), where the velocity block carries upper-convected stretching terms and a cell-scale stress layer. Same problem, one key changed: sor -> one step took 26947 s and the drag reached -4642 asm -> the same step took 115 s and stayed finite 90% of the inner velocity solves hit their iteration cap under sor, 77% under asm. It then produced the SAME wrong drag as every other smoother, which is how the real fault was identified as being outside the linear algebra. Cost on a benign operator is smaller than expected. SolCx (free-slip, viscosity jump), identical velocity error to four significant figures in every case: 32^2, contrast 1e6: default 15.5 s fast 11.3 s bulletproof 12.0 s 32^2, contrast 1e8: default 15.9 s fast 10.6 s bulletproof 12.6 s 48^2, 3 levels, 1e6: default 157 s fast 112 s bulletproof 136 s Faster than the current default and 10-20% behind "fast", so it is an opt-in and not a change of default. Every variant still writes the same option keys, which is what keeps the stale-key derivation variant-independent; this is why the variant takes PETSc's own Schwarz defaults rather than tuning overlap. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL --- .../cython/petsc_generic_snes_solvers.pyx | 30 +++++++++--- .../utilities/multigrid_options.py | 46 +++++++++++++++---- 2 files changed, 61 insertions(+), 15 deletions(-) diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index 8cbe26c54..1cf41317c 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -365,10 +365,14 @@ class SolverBaseClass(uw_object): """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" + converge" / "I need the solver off the list of suspects"); 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.""" + strategy = getattr(self, "_strategy", "default") + if strategy in ("fast", "bulletproof"): + return strategy + return "robust" def _push_managed_option(self, key, value): """Write a PETSc option UW3 owns, recording that we wrote it. @@ -7043,7 +7047,7 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): def strategy(self): """ What this solve should optimise for — the named intent over the - multigrid smoother's two measured regimes. + multigrid smoother's measured regimes. - ``"default"``, ``"robust"``: ``gmres``/4 smoothing. Survives an operator a stationary smoother stalls on: Spiegelman notch (:math:`\eta` contrast @@ -7056,6 +7060,18 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): 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. + - ``"bulletproof"``: ``gmres``/8 smoothing preconditioned by **additive + Schwarz**, whose subdomain solves invert the local coupling directly rather + than relying on the operator being close to symmetric or elliptic. Slower + per sweep, and it exists to take the linear solve OFF the list of suspects. + When a run misbehaves and the solver, the transport scheme and the mesh + resolution are all candidates, run it again under ``"bulletproof"``: if the + answer is unchanged, the fault is not in the linear algebra. Measured on + creeping Oldroyd-B past a confined cylinder (Wi 0.4, three levels), where + ``"robust"`` spent 26947 s on one step and reached a drag of -4642 while + ``"bulletproof"`` took 115 s and stayed finite — and produced the same + wrong drag as every other smoother, which is what identified the real fault + as the under-resolved elastic layer rather than the solve. ``"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. @@ -7100,10 +7116,10 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): @strategy.setter def strategy(self, value): - if value not in ("default", "robust", "fast"): + if value not in ("default", "robust", "fast", "bulletproof"): raise ValueError( f"Unknown solver strategy {value!r}: " - "expected 'default', 'robust', or 'fast'." + "expected 'default', 'robust', 'fast', or 'bulletproof'." ) # 'fast' and 'robust' now select a real smoother variant, via # `_mg_smoother_variant` -> `multigrid_options.geometric_mg_bundle`. They were diff --git a/src/underworld3/utilities/multigrid_options.py b/src/underworld3/utilities/multigrid_options.py index 3651cccce..d2fbc4c64 100644 --- a/src/underworld3/utilities/multigrid_options.py +++ b/src/underworld3/utilities/multigrid_options.py @@ -83,11 +83,13 @@ 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: +#: Smoother variants. These are 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") +#: stalls, ``"fast"`` is cheaper per cycle where the operator is benign, and +#: ``"bulletproof"`` exists to take the linear solve OFF the list of suspects when a +#: run misbehaves and the cause could be the solver, the transport scheme or the +#: resolution. Selected by ``solver.strategy``; see :func:`geometric_mg_bundle`. +GEOMETRIC_MG_SMOOTHERS = ("robust", "fast", "bulletproof") class MGSettings(NamedTuple): @@ -183,8 +185,10 @@ def describe(settings, levels=None, overridden=()): 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.""" + Every variant sets the SAME keys — only values differ — so the derived stale-key + sets are variant-independent. A variant that needs a key the others do not have + must add it to all of them (with a default), or the invariant breaks: this is why + ``"bulletproof"`` takes PETSc's own ASM defaults rather than tuning overlap.""" settings = { # The KSP this preconditioner serves must be FLEXIBLE, for both # variants. The "robust" smoother is a Krylov solve, so the @@ -253,6 +257,32 @@ def _geometric_mg_settings(coarse, smoother="robust"): # 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 + elif smoother == "bulletproof": + # Additive Schwarz under the Krylov smoother: each subdomain solve inverts + # the LOCAL coupling directly, so the smoother does not rely on the operator + # being close to symmetric or elliptic. SOR does, and has no purchase once it + # is neither. + # + # The measurement: creeping Oldroyd-B past a confined cylinder (Wi 0.4, + # beta 0.59, three levels), where the viscoelastic velocity block carries + # upper-convected stretching terms and a cell-scale stress layer. Same + # problem, one key changed: + # + # sor -> step 15 took 26947 s and the drag reached -4642 + # asm -> step 15 took 115 s and the drag stayed finite + # + # 90% of the inner velocity solves hit their iteration cap under sor against + # 77% under asm. It is SLOWER per sweep on a benign operator, which is why + # this is a third named variant and not a change to "robust". + # + # Its purpose is diagnostic as much as numerical. When a run misbehaves and + # the linear solve, the transport scheme and the mesh resolution are all + # candidates, this variant removes the first from the list: the cylinder run + # above survived under asm and still produced the same wrong drag at the same + # step, which is how we learned the fault was not in the solver. + settings["mg_levels_ksp_type"] = "gmres" + settings["mg_levels_ksp_max_it"] = 8 + settings["mg_levels_pc_type"] = "asm" else: raise ValueError( f"smoother must be one of {GEOMETRIC_MG_SMOOTHERS} (got {smoother!r})") @@ -355,7 +385,7 @@ def geometric_mg_bundle(coarse="redundant", smoother="robust"): Parameters ---------- - smoother : {"robust", "fast"} + smoother : {"robust", "fast", "bulletproof"} 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. @@ -372,7 +402,7 @@ def geometric_mg_bundle(coarse="redundant", smoother="robust"): Notes ----- - The two smoother variants are two measured regimes, and neither dominates. + The smoother variants are measured regimes, and none dominates. ``"robust"`` (``gmres``/4) survives an operator a stationary smoother stalls on. Spiegelman notch, :math:`\eta` contrast 1e26, nested 4-level hierarchy: