From d99e87230f20305d6cffd72afa07deb197e58d99 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Thu, 21 May 2026 18:12:03 +1000 Subject: [PATCH] Add SNES_Scalar.petsc_use_constant_nullspace for singular Poisson MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Scalar problems with no Dirichlet boundary conditions have a constant kernel — Poisson on a fully-Neumann box, on a closed manifold (sphere, torus), on a periodic domain. The linear operator is singular and the solver either diverges or returns an arbitrary (initial-guess-dependent) solution within the constant null space. The new opt-in flag attaches MatNullSpace(constant=True) to the Jacobian operator (and its transpose, and the preconditioner) at setup. PETSc projects each KSP RHS onto the orthogonal complement of the nullspace before solving and returns the minimum-norm (near-zero-mean) solution. Default False, so every existing scalar solver is unaffected. Verified on a pure-Neumann 2D box: pushes the discrete mean of the solution from -0.27 (arbitrary) to -0.04 (canonical); non-constant part of the solution is unchanged. Test at tests/test_1055_snes_scalar_constant_nullspace.py. Underworld development team with AI support from Claude Code --- .../cython/petsc_generic_snes_solvers.pyx | 72 ++++++++++ ...est_1055_snes_scalar_constant_nullspace.py | 129 ++++++++++++++++++ 2 files changed, 201 insertions(+) create mode 100644 tests/test_1055_snes_scalar_constant_nullspace.py diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index 6bb6b71b..1b7026e6 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/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, diff --git a/tests/test_1055_snes_scalar_constant_nullspace.py b/tests/test_1055_snes_scalar_constant_nullspace.py new file mode 100644 index 00000000..52362db1 --- /dev/null +++ b/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")