fix(mg): lock out native geometric FMG for single-field solvers (#276)

src/underworld3/cython/petsc_generic_snes_solvers.pyx

@@ -473,6 +473,32 @@ class SolverBaseClass(uw_object):
         else:  # "gamg" — explicit, always applied
             want_fmg = False
 
+        # Native geometric FMG relies on DMCreateInjection between the refined
+        # DMPlex levels. For a **single-field** (scalar / vector) discretisation
+        # this is fragile: whether PETSc can build the injection depends on the
+        # geometry×element-degree×refinement combination, and it fails at solve
+        # time with err62 "Could not locate matching functional for injection"
+        # on the common curved-shell cases (issue #276) as well as some flat
+        # high-degree ones. The Stokes velocity sub-block (prefix
+        # "fieldsplit_velocity_") is the validated, robust native-FMG path and is
+        # unaffected. So never auto-route a single-field solver to native FMG —
+        # fall back to GAMG. Geometric MG on a scalar/vector solver is available,
+        # robustly, via ``underworld3.utilities.custom_mg.set_custom_fmg`` (own
+        # barycentric/RBF prolongation + Galerkin coarse operators; no injection).
+        if want_fmg and prefix == "":
+            if self._preconditioner == "fmg" and uw.mpi.rank == 0:
+                import warnings
+                warnings.warn(
+                    f"[{self.name}] preconditioner='fmg' is not supported on a "
+                    f"single-field (scalar/vector) solver: native geometric FMG "
+                    f"needs DMCreateInjection, which PETSc cannot reliably build "
+                    f"on a refined DMPlex (issue #276). Falling back to GAMG. For "
+                    f"geometric MG on this solver use "
+                    f"underworld3.utilities.custom_mg.set_custom_fmg().",
+                    stacklevel=2,
+                )
+            want_fmg = False
+
         if want_fmg:
             # Geometric Full Multigrid on the refinement hierarchy. Galerkin
             # (RAP) coarse operators are required because UW3 does not install
@@ -500,7 +526,7 @@ class SolverBaseClass(uw_object):
                 opts.delValue(f"{prefix}{key}")
             self._pc_managed_value = "mg"
         else:
-            if self._preconditioner == "fmg" and uw.mpi.rank == 0:
+            if self._preconditioner == "fmg" and n_levels <= 1 and uw.mpi.rank == 0:
                 import warnings
                 warnings.warn(
                     f"[{self.name}] preconditioner='fmg' requested but the mesh "

tests/test_1014_stokes_multigrid.py

@@ -91,15 +91,20 @@ def test_invalid_preconditioner_raises():
         stokes.preconditioner = "wibble"
 
 
-def test_scalar_poisson_auto_geometric_mg():
+def test_scalar_poisson_auto_falls_back_to_gamg():
+    # #276: native geometric FMG is locked out for single-field (scalar/vector)
+    # solvers — DMCreateInjection is not reliably constructible on a refined
+    # DMPlex for a single field (fails on curved shells and some high-degree flat
+    # cases). So a scalar solver on a refined hierarchy falls back to GAMG rather
+    # than crashing; robust geometric MG for scalars is via custom_mg.set_custom_fmg.
     poisson = uw.systems.Poisson(mesh_refined)
     poisson.constitutive_model = uw.constitutive_models.DiffusionModel
     poisson.constitutive_model.Parameters.diffusivity = 1
     poisson.f = 0.0
     poisson.add_dirichlet_bc(0.0, "Bottom")
     poisson.add_dirichlet_bc(1.0, "Top")
     poisson.solve()
-    assert poisson.petsc_options.getString("pc_type") == "mg"
+    assert poisson.petsc_options.getString("pc_type") == "gamg"
     assert poisson.snes.getConvergedReason() > 0
 
 

tests/test_1020_fmg_single_field_lockout.py

@@ -0,0 +1,94 @@
+#!/usr/bin/env python3
+"""Regression: native geometric FMG is locked out for single-field solvers (#276).
+
+Native geometric FMG relies on ``DMCreateInjection`` between refined DMPlex
+levels. PETSc can build that for the Stokes velocity sub-block but NOT for a
+single-field (scalar/vector) discretisation on a refined DMPlex — it fails at
+solve time with err62 ("Could not locate matching functional for injection").
+
+So ``preconditioner="fmg"``/``"auto"`` on a scalar/vector solver must fall back
+to GAMG (and solve), never route to native FMG and crash. Geometric MG on such
+a solver remains available via ``utilities.custom_mg.set_custom_fmg`` (covered by
+test_1016).
+"""
+import sympy
+import pytest
+
+import underworld3 as uw
+
+pytestmark = [pytest.mark.level_2, pytest.mark.tier_a]
+
+
+@pytest.fixture
+def annulus_hierarchy():
+    # refinement>=1 => a multi-level dm_hierarchy (the FMG trigger)
+    return uw.meshing.Annulus(
+        radiusOuter=1.0, radiusInner=0.5, cellSize=1.0 / 6.0, refinement=2, qdegree=3
+    )
+
+
+def _poisson(mesh):
+    T = uw.discretisation.MeshVariable("T", mesh, 1, degree=3)
+    p = uw.systems.Poisson(mesh, T)
+    p.constitutive_model = uw.constitutive_models.DiffusionModel
+    p.constitutive_model.Parameters.diffusivity = 1.0
+    p.add_dirichlet_bc(1.0, "Lower")
+    p.add_dirichlet_bc(0.0, "Upper")
+    return p
+
+
+def test_scalar_fmg_falls_back_and_solves(annulus_hierarchy):
+    """preconditioner='fmg' on a scalar solver falls back to GAMG and solves
+    (previously err62 DMCreateInjection)."""
+    assert len(annulus_hierarchy.dm_hierarchy) > 1  # hierarchy present
+    p = _poisson(annulus_hierarchy)
+    p.preconditioner = "fmg"
+    p.solve()  # must not raise err62
+    assert p.snes.getConvergedReason() > 0
+    # resolved to GAMG, NOT native geometric mg
+    assert p.petsc_options.getString("pc_type") == "gamg"
+
+
+def test_scalar_auto_uses_gamg_not_native_fmg(annulus_hierarchy):
+    """preconditioner='auto' on a scalar solver with a hierarchy must NOT select
+    native FMG (it would crash) — it silently uses GAMG."""
+    p = _poisson(annulus_hierarchy)
+    # 'auto' is the default; solve and confirm no native-mg pc_type
+    p.solve()
+    assert p.snes.getConvergedReason() > 0
+    assert p.petsc_options.getString("pc_type") != "mg"
+
+
+def test_cartesian_box_scalar_also_falls_back():
+    """The lockout is single-field-wide, not geometry-specific: native single-
+    field FMG is fragile across geometry×degree×refinement (it errors 62 even on
+    some flat high-degree boxes), so a scalar solver on a flat box also falls
+    back to GAMG and solves rather than risk the crash."""
+    box = uw.meshing.UnstructuredSimplexBox(
+        minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=0.4, refinement=1, qdegree=3
+    )
+    assert len(box.dm_hierarchy) > 1
+    T = uw.discretisation.MeshVariable("Tbox", box, 1, degree=3)
+    p = uw.systems.Poisson(box, T)
+    p.constitutive_model = uw.constitutive_models.DiffusionModel
+    p.constitutive_model.Parameters.diffusivity = 1.0
+    p.add_dirichlet_bc(0.0, "Bottom")
+    p.add_dirichlet_bc(1.0, "Top")
+    p.preconditioner = "fmg"
+    p.solve()  # must not raise err62
+    assert p.snes.getConvergedReason() > 0
+    assert p.petsc_options.getString("pc_type") == "gamg"
+
+
+def test_stokes_velocity_fmg_still_selected(annulus_hierarchy):
+    """The lockout is single-field only: the Stokes velocity sub-block still gets
+    native geometric FMG on the same hierarchy."""
+    v = uw.discretisation.MeshVariable("v", annulus_hierarchy, annulus_hierarchy.dim, degree=2)
+    pp = uw.discretisation.MeshVariable("p", annulus_hierarchy, 1, degree=1)
+    stokes = uw.systems.Stokes(annulus_hierarchy, velocityField=v, pressureField=pp)
+    stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel
+    stokes.constitutive_model.Parameters.viscosity = 1.0
+    stokes.preconditioner = "fmg"
+    stokes._build(False, False, None)
+    # the velocity sub-block keeps native geometric multigrid
+    assert stokes.petsc_options.getString("fieldsplit_velocity_pc_type") == "mg"