Add SNES_Scalar.petsc_use_constant_nullspace for singular Poisson

src/underworld3/cython/petsc_generic_snes_solvers.pyx

@@ -1646,6 +1646,15 @@ class SNES_Scalar(SolverBaseClass):
         self.boundary_conditions = False
         # self._constitutive_model = None
 
+        # Constant-nullspace handling. When True, attach a MatNullSpace
+        # with the constant mode to the Jacobian operator before each
+        # KSP solve. Needed for scalar problems on a closed manifold or
+        # a fully-Neumann domain where the operator (e.g. Laplacian)
+        # has a 1-D constant kernel and the linear system is otherwise
+        # singular. Default False — every solver with any Dirichlet BC
+        # or reaction term is non-singular and shouldn't pay the cost.
+        self._petsc_use_constant_nullspace = False
+
         self.verbose = verbose
 
         self._rebuild_after_mesh_update = self._build  # Maybe just reboot the dm
@@ -1656,6 +1665,30 @@ class SNES_Scalar(SolverBaseClass):
 
         self.is_setup = False
 
+    @property
+    def petsc_use_constant_nullspace(self):
+        """Whether to attach a constant MatNullSpace to the Jacobian.
+
+        Set to ``True`` for scalar problems whose linear operator has a
+        constant kernel and no Dirichlet boundary conditions to break
+        it — fully-Neumann domains or closed manifolds (e.g. a Poisson
+        problem on a periodic box or on a spherical surface). PETSc
+        projects the right-hand side onto the orthogonal complement of
+        the nullspace and selects the minimum-norm solution from the
+        null-affine family, so the system becomes uniquely solvable up
+        to that nullspace.
+
+        Default ``False`` — any solver with at least one Dirichlet BC
+        or a reaction term is non-singular and shouldn't pay the
+        projection cost.
+        """
+        return self._petsc_use_constant_nullspace
+
+    @petsc_use_constant_nullspace.setter
+    def petsc_use_constant_nullspace(self, value):
+        self._petsc_use_constant_nullspace = bool(value)
+        self.is_setup = False
+
     @property
     def tolerance(self):
         """
@@ -2081,9 +2114,48 @@ class SNES_Scalar(SolverBaseClass):
 
             UW_DMPlexSetSNESLocalFEM(cdm.dm, PETSC_FALSE, NULL)
 
+        if self._petsc_use_constant_nullspace:
+            self._attach_constant_nullspace()
+
         self.is_setup = True
         self.constitutive_model._solver_is_setup = True
 
+    def _attach_constant_nullspace(self):
+        """Attach a constant ``MatNullSpace`` to the Jacobian.
+
+        Calls ``snes.setUp()`` to materialise the Jacobian template,
+        then sets a constant-mode nullspace on the operator and
+        preconditioner matrices (and their transposes — the projector
+        is symmetric). PETSc projects each KSP right-hand side onto
+        the orthogonal complement of the nullspace before solving and
+        returns the minimum-norm solution within the affine null
+        space.
+
+        Used by scalar Poisson on a fully-Neumann domain or a closed
+        manifold. See :attr:`petsc_use_constant_nullspace`.
+        """
+        self.snes.setUp()
+        jacobian = self.snes.getJacobian()
+        operator_matrix = jacobian[0]
+        preconditioner_matrix = jacobian[1] if len(jacobian) > 1 else None
+
+        nullspace = PETSc.NullSpace().create(
+            constant=True, vectors=(), comm=self.dm.comm,
+        )
+
+        operator_matrix.setNullSpace(nullspace)
+        operator_matrix.setTransposeNullSpace(nullspace)
+        if (preconditioner_matrix is not None
+                and preconditioner_matrix.handle != operator_matrix.handle):
+            preconditioner_matrix.setNullSpace(nullspace)
+            preconditioner_matrix.setTransposeNullSpace(nullspace)
+
+        if self.verbose and uw.mpi.rank == 0:
+            print(
+                f"SNES_Scalar ({self.name}): attached constant nullspace",
+                flush=True,
+            )
+
     @timing.routine_timer_decorator
     def solve(self,
               zero_init_guess: bool =True,

tests/test_1055_snes_scalar_constant_nullspace.py

@@ -0,0 +1,129 @@
+"""Regression test for ``SNES_Scalar.petsc_use_constant_nullspace``.
+
+A scalar Poisson problem with no Dirichlet boundary conditions has a
+constant kernel (any constant solves :math:`-\\Delta u = 0`). The linear
+system is singular and the solution is determined only up to an additive
+constant.
+
+The new ``petsc_use_constant_nullspace`` flag tells PETSc about the
+constant nullspace via ``MatNullSpace().create(constant=True)``. PETSc
+projects each KSP right-hand side onto the orthogonal complement of the
+nullspace before solving and returns the minimum-norm solution within
+the affine null space — i.e. the solution with (discrete) zero mean.
+
+We verify:
+  1. The shape of the solution (non-constant part) is the same with or
+     without the flag — PETSc's projection lives entirely in the
+     constant kernel, so the non-constant component is unaffected.
+  2. With the flag, the discrete mean of the solution is small
+     (minimum-norm pinning).
+  3. The recovered field agrees with the analytic solution
+     :math:`u(x, y) = -x^3/6 + x^2/4 + C` to FE-discretisation
+     accuracy.
+"""
+import numpy as np
+import pytest
+import sympy
+
+import underworld3 as uw
+from underworld3.systems import Poisson
+
+
+def _build_neumann_poisson(mesh, use_nullspace):
+    """Build a pure-Neumann Poisson with f = x - 1/2 (zero mean)."""
+    T = uw.discretisation.MeshVariable(
+        f"T_ns_{use_nullspace}", mesh, num_components=1, degree=2,
+    )
+    poisson = Poisson(mesh, u_Field=T)
+    poisson.constitutive_model = uw.constitutive_models.DiffusionModel
+    poisson.constitutive_model.Parameters.diffusivity = 1.0
+    poisson.petsc_use_constant_nullspace = use_nullspace
+    x, _ = mesh.X
+    poisson.f = x - sympy.Rational(1, 2)
+    return poisson, T
+
+
+def _shape_error(T, mesh):
+    """Discretisation error in the non-constant part of the solution."""
+    coords = np.asarray(T.coords)
+    expected = -coords[:, 0] ** 3 / 6.0 + coords[:, 0] ** 2 / 4.0
+    diff = T.data[:, 0] - expected
+    diff_centered = diff - diff.mean()
+    return np.abs(diff_centered).max()
+
+
+def test_constant_nullspace_canonical_solution():
+    """With ``petsc_use_constant_nullspace=True`` the solver returns the
+    minimum-norm (near-zero-mean) solution."""
+    mesh = uw.meshing.UnstructuredSimplexBox(
+        minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=0.1,
+    )
+
+    poisson_off, T_off = _build_neumann_poisson(mesh, use_nullspace=False)
+    poisson_off.solve()
+    assert poisson_off.snes.getConvergedReason() > 0, (
+        "Pure-Neumann Poisson without nullspace declaration failed to "
+        f"converge: SNES reason {poisson_off.snes.getConvergedReason()}"
+    )
+
+    poisson_on, T_on = _build_neumann_poisson(mesh, use_nullspace=True)
+    poisson_on.solve()
+    assert poisson_on.snes.getConvergedReason() > 0, (
+        "Pure-Neumann Poisson with nullspace declaration failed to "
+        f"converge: SNES reason {poisson_on.snes.getConvergedReason()}"
+    )
+
+    # Both should recover the analytic shape to FE accuracy.
+    err_off = _shape_error(T_off, mesh)
+    err_on  = _shape_error(T_on,  mesh)
+    assert err_off < 1.0e-4, f"shape error off = {err_off:.3e}"
+    assert err_on  < 1.0e-4, f"shape error on  = {err_on:.3e}"
+
+    # The flag should leave the non-constant part of the solution
+    # essentially unchanged — PETSc's projection lives in the kernel.
+    assert abs(err_off - err_on) < 1.0e-5, (
+        f"shape error differs more than expected between paths: "
+        f"off={err_off:.3e}, on={err_on:.3e}"
+    )
+
+    # With the flag, the mean of the solution should be much closer to
+    # zero than without (canonical minimum-norm pinning). The exact
+    # value depends on the discrete inner product against the constant
+    # mode, but the flag should move it toward zero by a clear margin.
+    mean_off = T_off.data[:, 0].mean()
+    mean_on  = T_on.data[:, 0].mean()
+    assert abs(mean_on) < abs(mean_off), (
+        f"with nullspace declaration, |mean(T)| should drop; "
+        f"got |mean_off|={abs(mean_off):.3e}, |mean_on|={abs(mean_on):.3e}"
+    )
+
+
+def test_constant_nullspace_property_setter():
+    """Property setter round-trip + ``is_setup`` invalidation."""
+    mesh = uw.meshing.UnstructuredSimplexBox(
+        minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=0.2,
+    )
+    T = uw.discretisation.MeshVariable(
+        "T_setter", mesh, num_components=1, degree=1,
+    )
+    poisson = Poisson(mesh, u_Field=T)
+
+    assert poisson.petsc_use_constant_nullspace is False, (
+        "Default should be False"
+    )
+
+    poisson.petsc_use_constant_nullspace = True
+    assert poisson.petsc_use_constant_nullspace is True
+    assert poisson.is_setup is False, (
+        "Setting the flag should invalidate the setup so the solver "
+        "rebuilds with the nullspace attached."
+    )
+
+    poisson.petsc_use_constant_nullspace = False
+    assert poisson.petsc_use_constant_nullspace is False
+
+
+if __name__ == "__main__":
+    test_constant_nullspace_property_setter()
+    test_constant_nullspace_canonical_solution()
+    print("PASSED")