From 3ca944949f54ec306474ac8fb782f56cc1dcfb89 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Sat, 13 Jun 2026 09:34:59 +1000 Subject: [PATCH] tests: cover SNES_Scalar.constant_nullspace (pure-Neumann Poisson) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Development's constant_nullspace (added in 074a660 for the Monge-Ampere smoother's pure-Neumann solve) had no test coverage. This adds a tier_a/level_1 regression test salvaged from the now-closed #201 (which proposed a duplicate petsc_use_constant_nullspace flag): - a pure-Neumann Poisson (-Δu = x - 1/2, no Dirichlet BCs) is singular up to a constant; with constant_nullspace=True PETSc returns the minimum-norm (discrete zero-mean) solution. The test checks convergence, that the non-constant part matches the analytic solution to FE accuracy, and that the mean drops vs the no-nullspace solve; - property-setter round-trip and default. Verified against current development: constant_nullspace=True gives mean 3.7e-16 vs -0.224 drift when off. Underworld development team with AI support from Claude Code --- tests/test_1056_constant_nullspace.py | 125 ++++++++++++++++++++++++++ 1 file changed, 125 insertions(+) create mode 100644 tests/test_1056_constant_nullspace.py diff --git a/tests/test_1056_constant_nullspace.py b/tests/test_1056_constant_nullspace.py new file mode 100644 index 00000000..f1414c97 --- /dev/null +++ b/tests/test_1056_constant_nullspace.py @@ -0,0 +1,125 @@ +"""Regression test for ``SNES_Scalar.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 ``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 + +pytestmark = [pytest.mark.tier_a, pytest.mark.level_1] + + +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.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 ``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 and default.""" + 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.constant_nullspace is False, "Default should be False" + + poisson.constant_nullspace = True + assert poisson.constant_nullspace is True + + poisson.constant_nullspace = False + assert poisson.constant_nullspace is False + + +if __name__ == "__main__": + test_constant_nullspace_property_setter() + test_constant_nullspace_canonical_solution() + print("PASSED")