diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index 8cbe26c5..1cf41317 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 3651cccc..d2fbc4c6 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: