Skip to content

Add SNES_Scalar.petsc_use_constant_nullspace for singular Poisson#201

Open
lmoresi wants to merge 1 commit into
developmentfrom
feature/snes-constant-nullspace
Open

Add SNES_Scalar.petsc_use_constant_nullspace for singular Poisson#201
lmoresi wants to merge 1 commit into
developmentfrom
feature/snes-constant-nullspace

Conversation

@lmoresi
Copy link
Copy Markdown
Member

@lmoresi lmoresi commented May 21, 2026

Summary

Adds an opt-in petsc_use_constant_nullspace flag on SNES_Scalar that attaches a constant MatNullSpace to the Jacobian. Needed for scalar problems whose linear operator has a constant kernel and no Dirichlet BC to break it: pure-Neumann domains, closed manifolds (sphere, torus, periodic boxes).

Without it, PETSc either diverges on the singular system or returns an arbitrary solution within the constant null space, dependent on initial guess and solver path. With it, PETSc projects each KSP right-hand side onto the orthogonal complement of the nullspace and returns the minimum-norm (near-zero-mean) solution.

Default False, so every existing scalar solver is unaffected.

Implementation

src/underworld3/cython/petsc_generic_snes_solvers.pyx:

  • _petsc_use_constant_nullspace = False field in SNES_Scalar.__init__
  • petsc_use_constant_nullspace property + setter (setter invalidates is_setup so the rebuild attaches the nullspace)
  • _attach_constant_nullspace() helper: calls snes.setUp() to materialise the Jacobian, then sets MatNullSpace().create(constant=True) on operator + preconditioner + their transposes
  • Single call site in SNES_Scalar._setup_solver

Mirrors the existing Stokes saddle-point nullspace pattern (_attach_stokes_nullspace) but specialised for the scalar case (no field decomposition, just the one Jacobian Mat).

Example

poisson = uw.systems.Poisson(box_mesh, u_Field=T)
poisson.constitutive_model = uw.constitutive_models.DiffusionModel
poisson.constitutive_model.Parameters.diffusivity = 1.0
poisson.petsc_use_constant_nullspace = True   # opt-in
poisson.f = x - 0.5    # zero-mean source — equation is solvable
poisson.solve()
# T has near-zero discrete mean and matches the analytic answer
# u(x,y) = -x^3/6 + x^2/4 + C  with C ≈ 0

Why this matters / context

Surfaced while testing the SphericalManifold work (cd-1 mesh, separate WIP branch). Surface Poisson -Δ_S T = z on a closed sphere fails DIVERGED_LINE_SEARCH at iteration 0 without the nullspace declaration; with the flag it solves cleanly and recovers T = z/2 + const to discretisation error with monotone refinement convergence. The same flag works on any pure-Neumann scalar problem, manifold or not — flagging this as a small standalone change keeps the API addition independent of the larger manifold-mesh work.

Test plan

  • tests/test_1055_snes_scalar_constant_nullspace.py — 2 tests:
    • test_constant_nullspace_canonical_solution: pure-Neumann Poisson on a 2-D box; verifies convergence on both paths, shape-correctness matches analytic to FE accuracy on both paths, and the flag pulls the discrete mean toward zero.
    • test_constant_nullspace_property_setter: property round-trip, setup invalidation.
  • Regression: test_0001_meshes.py, test_0506_tensor_evaluate.py, test_1010_stokesCart.py (25/25 Stokes tests) all pass.

Underworld development team with AI support from Claude Code

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
Copilot AI review requested due to automatic review settings May 21, 2026 08:12
Copy link
Copy Markdown
Contributor

Copilot AI left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Pull request overview

This PR adds an opt-in PETSc constant-nullspace attachment for SNES_Scalar to improve robustness and solution selection for singular scalar operators (e.g., pure-Neumann Poisson problems and closed-manifold Laplacians) by informing PETSc of the constant kernel via MatNullSpace(constant=True).

Changes:

  • Introduces SNES_Scalar.petsc_use_constant_nullspace (default False) and wires it into solver setup via a new _attach_constant_nullspace() helper.
  • Adds a regression test covering the new flag on a pure-Neumann Poisson problem and validating the property setter behavior.

Reviewed changes

Copilot reviewed 2 out of 2 changed files in this pull request and generated 4 comments.

File Description
src/underworld3/cython/petsc_generic_snes_solvers.pyx Adds the petsc_use_constant_nullspace property and attaches a constant MatNullSpace to the scalar Jacobian.
tests/test_1055_snes_scalar_constant_nullspace.py New regression tests for solver convergence/shape and the flag/property behavior.

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

Comment on lines +2117 to +2119
if self._petsc_use_constant_nullspace:
self._attach_constant_nullspace()

Comment on lines +2142 to +2144
nullspace = PETSc.NullSpace().create(
constant=True, vectors=(), comm=self.dm.comm,
)
Comment on lines +89 to +98
# 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}"
)
Comment on lines +1676 to +1680
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.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants