Add a "bulletproof" smoother strategy: Schwarz, to take the solver off the suspect list

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

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: