Add PETSc pressure nullspace support to Stokes

[gthyagi]

Summary

This PR adds PETSc constant-pressure nullspace support to the Underworld Stokes saddle-point solver.

What changed

• Added a new petsc_use_pressure_nullspace property on SNES_Stokes_SaddlePt.
• When enabled, the solver builds a coupled nullspace vector with zero velocity entries and constant pressure entries.
• The nullspace is attached to the Stokes Jacobian and transpose Jacobian before each solve path, including after SNES reconfiguration.
• The cached nullspace objects are cleared when the solver DM is rebuilt.
• The solver now raises a clear error if pressure Dirichlet BCs are present while pressure-nullspace handling is enabled.
• Added a regression test that solves a Stokes problem without any pressure BC and checks that the Jacobian has a nullspace attached.

Why this is important

At the moment, users often make singular Stokes systems solvable by imposing a pressure Dirichlet condition, sometimes on an entire boundary. That removes the pressure gauge freedom, but it is stronger than the mathematical requirement and can change the solution.

This PR makes it possible to follow the standard PETSc approach instead:

• solve the singular Stokes system with the constant-pressure nullspace registered algebraically, and
• choose the pressure gauge afterwards by subtracting a mean or applying another normalization.

This matters for annulus and spherical-shell benchmarks, especially cases like the Thieulot benchmark k=0 branch, where pressure is only defined up to a constant and a boundary-wide pressure Dirichlet condition is not the right tool.

It also makes Underworld align better with PETSc/ASPECT-style pressure handling for incompressible Stokes problems.

PETSc precedent

This follows the same idea used in PETSc Stokes examples, where PETSc builds a pressure nullspace and attaches it to the operator instead of forcing pressure pointwise on a boundary:

• ex62.c CreatePressureNullSpace()
• ex69.c CreatePressureNullSpace()
• MatSetNullSpace manual page

Validation

• pixi run -e amr-dev pytest tests/test_1013_stokes_pressure_nullspace.py -q

This test passes locally and verifies that a Stokes solve without pressure Dirichlet BCs converges with the nullspace attached.

[lmoresi]

Review: Good feature, a few suggestions

This is a useful addition — the PETSc nullspace approach is the right way to handle pressure gauge freedom, especially for annulus and spherical-shell benchmarks where pressure Dirichlet BCs are not physically meaningful.

Design question: rigid-body modes

The current implementation handles the constant-pressure nullspace (1 mode). For spherical shell / annulus models, the system can also have rigid-body rotation null modes (net rotation for free-slip shells). Does the architecture support extending the nullspace basis to include these? The PETSc.NullSpace().create(vectors=(...)) call accepts multiple basis vectors, so it should be straightforward to add rotation modes later — worth noting in the docstring.

Suggestions

1. Branch naming: The project uses descriptive kebab-case branch names like feature/pressure-nullspace or bugfix/... — see docs/developer/guides/branching-strategy.md for the convention. The codex/ prefix doesn't follow this pattern.

2. Test marker: The test uses @pytest.mark.level_3 but it's a 3x3 mesh solve that runs in seconds. level_1 would be more appropriate (see docs/developer/TESTING-RELIABILITY-SYSTEM.md).

3. _pressure_dirichlet_bcs field ID: Hardcoding {1} as the pressure field ID is fragile. Consider using self.p.field_id exclusively (with a fallback error if not available), rather than assuming 1 as a default.

4. Multiple _attach_pressure_nullspace calls: The call appears in 4 places in the solve paths. This is fine (cheap operation), but a comment explaining why it's needed after each setFromOptions() would help future readers.

Summary

The implementation follows PETSc best practices (ex62.c, ex69.c patterns). The validation check for conflicting pressure BCs is a nice touch. Good to merge after the minor items above.

[gthyagi]

Addressing the review items directly:

1. Branch naming

• I renamed the working branch on the fork to feature/stokes-nullspace.
• GitHub does not let me change the head branch of an existing PR, so #90 itself still points to codex/petsc-pressure-nullspace.
• To reopen this review thread, I restored codex/petsc-pressure-nullspace as an alias that points to the same commit as feature/stokes-nullspace.
• So the naming convention is implemented on the branch itself, but this specific PR cannot be retargeted to the renamed head branch.

2. Rigid-body modes design question

• Implemented.
• The Stokes nullspace support is now generalized beyond the constant-pressure mode.
• A new petsc_velocity_nullspace_basis property allows multiple exact velocity nullspace modes to be registered alongside the pressure mode.
• The docstrings now document the annulus rigid rotation (-y, x) and the three spherical-shell rigid rotations (0, -z, y), (z, 0, -x), and (-y, x, 0), plus how to set them on Stokes.

3. @pytest.mark.level_3

• Not changed intentionally.
• This test remains level_3 because it executes a real Stokes solve. The mesh is small and runtime is short, but by the project testing definitions it is still a solver / physics test rather than a lightweight setup-only test.

4. _pressure_dirichlet_bcs field id

• Implemented.
• The old hardcoded {1} fallback is removed.
• _pressure_dirichlet_bcs() now uses self.p.field_id exclusively and raises a clear error if the pressure field is unavailable.

5. Repeated nullspace attach calls after setFromOptions()

• Implemented.
• The repeated reattach calls remain, and there is now an inline comment explaining that PETSc may rebuild operator state after setFromOptions(), so the nullspace must be reattached before each solve path.
• The comment is placed once before the first repeated call rather than duplicated at every identical call site.

6. PETSc-style validation against conflicting pressure BCs

• Implemented.
• The solver now checks for pressure Dirichlet BCs before attaching the Stokes nullspace and raises a clear error if they are present.

7. Extra validation added

• Implemented.
• The original pressure-nullspace regression test was kept and strengthened.
• A new shell test covers pressure plus rigid-rotation nullspace bases for a 2-D annulus and a 3-D spherical shell.

The current fork branches feature/stokes-nullspace and codex/petsc-pressure-nullspace both point to the same latest commit so this reopened PR now includes the recent implementation updates as well.

[gthyagi]

Closing this PR in favor of #91, which carries the same implementation on the correctly named head branch : #91 . I also copied the reviewer-response summary from this thread into #91 so the follow-up context is preserved there.

[gthyagi]

Closing this PR in favor of #91, which carries the same implementation on the correctly named head branch feature/stokes-nullspace: #91 . I also copied the reviewer-response summary from this thread into #91 so the follow-up context is preserved there.