From b297a3e357b60f591570af2b74164f9507ba14ab Mon Sep 17 00:00:00 2001 From: lmoresi Date: Sun, 9 Aug 2026 10:02:00 +1000 Subject: [PATCH] The geometric-MG bundle owns its outer-KSP pairing: fgmres, applied to the live KSP (#514) Since #471 the bundle's default smoother is gmres/4. A Krylov smoother makes the preconditioner vary between applications, and a plain-gmres outer's recurrence assumes it fixed: measured on a 3-D adapt child (PR #510 CI, then locally), the KSP reported CONVERGED_RTOL at 3 iterations with the true residual stalled at 1.3e-6 - a hundred times the gate. The same failure was found and fixed locally twice before, on the Stokes velocity sub-KSP and on the free-surface solver, which is exactly the piecemeal drift #471 exists to end; this puts the pairing in the one owner. Three parts. ksp_type=fgmres joins the geometric bundle for BOTH smoother variants - richardson does not need flexibility, but the stale-key derivation requires every variant to own the same keys, and fgmres judges convergence on the true residual, which is never wrong. The three solver-default writers of ksp_type go through _push_managed_option, so the bundle may upgrade the DEFAULT but never a user's explicit choice (the ownership latch). And the top-level install path applies the resolved type to the LIVE KSP: the database write alone is provably inert there, because the KSP consumed its options long before injection - the new live test caught exactly that in its first run. The fieldsplit velocity sub-KSP needs no live application (its setFromOptions runs at the parent's PCSetUp, after the write) and the rotated path already writes fgmres unconditionally on its per-solve prefix. On hierarchies where one FMG application converges the cycle, the two pairings are indistinguishable (measured: 3.3e-8 either way on the 2-D fixture); the drift needs a hierarchy the smoother varies on, which is the 3-D one-level-per-doubling shape arriving with PR #510 - its test_0842 is the discriminating gate, and passes on the composed tree. Closes #514. Underworld development team with AI support from Claude Code --- .../cython/petsc_generic_snes_solvers.pyx | 6 +- src/underworld3/utilities/custom_mg.py | 18 ++++- .../utilities/multigrid_options.py | 15 ++++ tests/test_1021_mg_option_bundle.py | 71 +++++++++++++++++++ 4 files changed, 104 insertions(+), 6 deletions(-) diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index 5d590b2f..5babe0e2 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -3076,7 +3076,7 @@ class SNES_Scalar(SolverBaseClass): self._pc_option_prefix = "" self.petsc_options["snes_type"] = "newtonls" - self.petsc_options["ksp_type"] = "gmres" + self._push_managed_option("ksp_type", "gmres") self._push_managed_option("pc_type", "gamg") self._push_managed_option("pc_gamg_type", "agg") self._push_managed_option("pc_gamg_repartition", True) @@ -3989,7 +3989,7 @@ class SNES_Vector(SolverBaseClass): # Here we can set some defaults for this set of KSP / SNES solvers self.petsc_options["snes_type"] = "newtonls" self.petsc_options["ksp_rtol"] = 1.0e-3 - self.petsc_options["ksp_type"] = "gmres" + self._push_managed_option("ksp_type", "gmres") self._push_managed_option("pc_type", "gamg") self._push_managed_option("pc_gamg_type", "agg") self._push_managed_option("pc_gamg_repartition", True) @@ -5001,7 +5001,7 @@ class SNES_MultiComponent(SolverBaseClass): # PETSc options — mirror SNES_Vector defaults self.petsc_options["snes_type"] = "newtonls" self.petsc_options["ksp_rtol"] = 1.0e-3 - self.petsc_options["ksp_type"] = "gmres" + self._push_managed_option("ksp_type", "gmres") self._push_managed_option("pc_type", "gamg") self._push_managed_option("pc_gamg_type", "agg") self._push_managed_option("pc_gamg_repartition", True) diff --git a/src/underworld3/utilities/custom_mg.py b/src/underworld3/utilities/custom_mg.py index 542c8791..9719e6e0 100644 --- a/src/underworld3/utilities/custom_mg.py +++ b/src/underworld3/utilities/custom_mg.py @@ -595,7 +595,8 @@ def _assert_no_zero_columns_serial(P_csr, level): f"operator would be singular.") -def _configure_pcmg(pc, Ps, coarse="redundant", smoother="robust", owned=None): +def _configure_pcmg(pc, Ps, coarse="redundant", smoother="robust", owned=None, + ksp=None): """Reconfigure ``pc`` as a fresh PCMG (FMG F-cycle) driven by the supplied reduced->reduced prolongations ``Ps``, Galerkin RAP for coarse operators. @@ -623,8 +624,19 @@ def _configure_pcmg(pc, Ps, coarse="redundant", smoother="robust", owned=None): fixed).""" nlev = len(Ps) + 1 prefix = pc.getOptionsPrefix() or "" + opts = PETSc.Options() multigrid_options.geometric_mg_bundle(coarse=coarse, smoother=smoother).apply( - PETSc.Options(), prefix, owned=owned) + opts, prefix, owned=owned) + # ``ksp_type`` is in the bundle (#514: a Krylov smoother makes this PC vary + # between applications, so its KSP must judge convergence flexibly), but on + # the top-level path the KSP consumed its options long before this + # injection runs, so a database write alone never takes effect. Apply the + # RESOLVED value — the bundle's, or the user's own where the ownership + # latch left it alone — to the live object. The fieldsplit velocity + # sub-KSP does not need this: its ``setFromOptions`` runs at the parent's + # ``PCSetUp``, after the write. + if ksp is not None: + ksp.setType(opts.getString(prefix + "ksp_type", ksp.getType())) pc.setType("mg") pc.setMGLevels(nlev) pc.setMGType(PETSc.PC.MGType.FULL) @@ -659,7 +671,7 @@ def _install_transfers(solver, Ps, verbose=False): ksp.setDMActive(PETSc.KSP.DMActive.OPERATOR, False) _configure_pcmg(ksp.getPC(), Ps, smoother=solver._mg_smoother_variant, - owned=solver._managed_pc_options) + owned=solver._managed_pc_options, ksp=ksp) if verbose: from underworld3 import mpi mpi.pprint(f"[{solver.name}] custom FMG installed: {nlev} levels, " diff --git a/src/underworld3/utilities/multigrid_options.py b/src/underworld3/utilities/multigrid_options.py index 2ba90a26..10a5c491 100644 --- a/src/underworld3/utilities/multigrid_options.py +++ b/src/underworld3/utilities/multigrid_options.py @@ -186,6 +186,21 @@ def _geometric_mg_settings(coarse, smoother="robust"): Both variants set the SAME keys — only values differ — so the derived stale-key sets are variant-independent.""" settings = { + # The KSP this preconditioner serves must be FLEXIBLE, for both + # variants. The "robust" smoother is a Krylov solve, so the + # preconditioner VARIES between applications and a plain gmres outer's + # recurrence has no guarantee: measured on a 3-D adapt child (#514), + # the preconditioned residual dropped nine orders and the KSP reported + # CONVERGED_RTOL at 3 iterations while the TRUE residual stalled at + # 1.3e-6. The same failure was found and fixed locally twice before — + # the Stokes velocity sub-KSP (see the smoother note below) and the + # free-surface solver (systems/solvers.py, "false-converges ... TRUE + # residual blew up") — which is exactly the drift-by-parallel-fixes + # this module exists to end. "fast" (richardson) is stationary and + # does not need flexibility, but fgmres is sound there too, and both + # variants setting the SAME keys is what keeps the stale-key + # derivation variant-independent. + "ksp_type": "fgmres", "pc_type": "mg", "pc_mg_type": "full", # FMG (F-cycle) # Galerkin (RAP) coarse operators are REQUIRED: UW3 installs no diff --git a/tests/test_1021_mg_option_bundle.py b/tests/test_1021_mg_option_bundle.py index d3122b6a..1b60f39e 100644 --- a/tests/test_1021_mg_option_bundle.py +++ b/tests/test_1021_mg_option_bundle.py @@ -48,6 +48,77 @@ def test_every_bundle_sets_the_smoother_iteration_count(): assert "pc_type" in bundle.settings +def test_the_geometric_bundle_pairs_krylov_smoothing_with_a_flexible_outer(): + """#514: a Krylov smoother makes the preconditioner VARY between + applications, and a plain-gmres outer's recurrence then has no guarantee. + Measured on a 3-D adapt child: the KSP reported CONVERGED_RTOL at 3 + iterations while the TRUE residual stalled 100x above the tolerance. + + The pairing lives in the bundle so it cannot be fixed piecemeal a THIRD + time — the Stokes velocity sub-KSP and the free-surface solver each + already fixed exactly this failure locally, which is the drift this + module exists to end. Both smoother variants set the key, because the + stale-key derivation requires every variant to own the same key set. + """ + for smoother in multigrid_options.GEOMETRIC_MG_SMOOTHERS: + bundle = multigrid_options.geometric_mg_bundle(smoother=smoother) + assert bundle.settings["ksp_type"] == "fgmres", ( + f"{smoother!r}: a geometric-MG solve must judge convergence " + f"flexibly; got {bundle.settings.get('ksp_type')}") + # ... and falling back to GAMG must not leave the flexible outer behind + # as an unexplained leftover on the prefix. + assert "ksp_type" in multigrid_options.gamg_bundle().stale + + +def test_an_adapt_child_scalar_solve_converges_truly(): + """The live form of the pairing test, on the route that caught #514. + + A scalar solve on an adapt child auto-injects the geometric hierarchy, so + its outer KSP must come out flexible ON THE LIVE OBJECT — the database + write alone is provably not enough, because the top-level KSP consumes its + options long before the injection runs (that inert write was the first + form of the fix, and this assert is what caught it). + + The KSP-type assert is the regression guard; the error bound is sanity + only. Measured on this fixture the two pairings are indistinguishable + (3.3e-8 either way, one FMG application) because the cycle converges + before the recurrence can drift — the hierarchy shape that separates them + is the 3-D subsampled one, guarded by test_0842's 1e-8 gate. + """ + import numpy as np + + mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=0.2, + regular=False, qdegree=2, refinement=1) + + def metric(points): + r = np.linalg.norm(np.asarray(points) - 0.5, axis=1) + return 1.0 / np.where(r < 0.25, 0.05, 0.2) ** 2 + + child = mesh.adapt(metric, max_levels=1) + + u = uw.discretisation.MeshVariable("u_flex", child, 1, degree=1) + poisson = uw.systems.Poisson(child, u_Field=u) + poisson.constitutive_model = uw.constitutive_models.DiffusionModel + poisson.constitutive_model.Parameters.diffusivity = 1.0 + poisson.f = 0.0 + poisson.add_dirichlet_bc(0.0, "Bottom") + poisson.add_dirichlet_bc(1.0, "Top") + poisson.petsc_options["ksp_rtol"] = 1e-8 + poisson.solve() + + ksp = poisson.snes.getKSP() + assert ksp.getPC().getType() == "mg", "the hierarchy was not injected" + assert ksp.getType() == "fgmres", ( + "a Krylov-smoothed geometric-MG solve is running under a plain gmres " + "outer; #514 has regressed") + err = np.linalg.norm(poisson.Unknowns.u.data[:, 0] + - poisson.Unknowns.u.coords[:, 1]) + nrm = np.linalg.norm(poisson.Unknowns.u.coords[:, 1]) + 1e-30 + assert err / nrm < 1e-6, ( + f"true relative error {err / nrm:.3e} on an exactly-linear solution") + + 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