From bfe13aba2f9cb8b6a0bf22966fe3eda38d58f70f Mon Sep 17 00:00:00 2001 From: Louis Moresi Date: Wed, 25 Mar 2026 16:26:06 +1100 Subject: [PATCH 1/3] Add nullspace_rotations property to mesh classes Mesh factories now set _nullspace_rotations with symbolic velocity fields for rigid-body rotation null modes: - Annulus, AnnulusWithSpokes, AnnulusInternalBoundary, DiscInternalBoundaries: 1 mode (z-rotation) - SphericalShell, SphericalShellInternalBoundary, CubedSphere, SegmentedSphericalShell, SegmentedSphericalBall: 3 modes (x, y, z rotation) - Box, QuarterAnnulus, SegmentofAnnulus, SegmentofSphere: 0 modes (walls break symmetry, default [] from base class) Each mode is a SymPy Matrix velocity field in Cartesian coordinates. The solver can project these onto its FE space for PETSc MatSetNullSpace. Underworld development team with AI support from Claude Code --- .../discretisation/discretisation_mesh.py | 30 +++++++++++++++++++ src/underworld3/meshing/annulus.py | 16 ++++++++++ src/underworld3/meshing/segmented.py | 16 ++++++++++ src/underworld3/meshing/spherical.py | 24 +++++++++++++++ 4 files changed, 86 insertions(+) diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py index 4b11dbd6e..b0985e3b8 100644 --- a/src/underworld3/discretisation/discretisation_mesh.py +++ b/src/underworld3/discretisation/discretisation_mesh.py @@ -664,6 +664,11 @@ def mesh_update_callback(array, change_context): self._lvec = None self.petsc_fe = None + # Rigid-body rotation null modes for this geometry. + # Mesh factories override this for closed surfaces (annulus, sphere). + # Each entry is a SymPy Matrix velocity field in mesh coordinates. + self._nullspace_rotations = [] + self.degree = degree self.qdegree = qdegree @@ -1531,6 +1536,31 @@ def t(self): """ return self._t + @property + def nullspace_rotations(self): + """Symbolic velocity fields for rigid-body rotation null modes. + + Returns a list of SymPy Matrix expressions in mesh Cartesian + coordinates. Empty for meshes with no rotation nullspace (boxes, + wedge segments with walls). Set by mesh factory functions for + closed surfaces (annulus, spherical shell, etc.). + + Each entry represents a rigid rotation: v = omega x r. + + Returns + ------- + list of sympy.Matrix + Velocity fields for each independent rotation mode. + + Examples + -------- + >>> annulus = uw.meshing.Annulus(...) + >>> annulus.nullspace_rotations # [Matrix([-y, x])] + >>> shell = uw.meshing.SphericalShell(...) + >>> shell.nullspace_rotations # 3 rotation matrices + """ + return self._nullspace_rotations + @property def r(self) -> Tuple[sympy.vector.BaseScalar]: r"""Tuple of coordinate scalars :math:`(x, y)` or :math:`(x, y, z)`. diff --git a/src/underworld3/meshing/annulus.py b/src/underworld3/meshing/annulus.py index 1a3631d67..cf8d73e91 100644 --- a/src/underworld3/meshing/annulus.py +++ b/src/underworld3/meshing/annulus.py @@ -538,6 +538,10 @@ class boundary_normals(Enum): new_mesh.boundary_normals = boundary_normals + # Full annulus: rigid rotation about z-axis + x, y = new_mesh.X + new_mesh._nullspace_rotations = [sympy.Matrix([-y, x])] + return new_mesh @@ -1110,6 +1114,10 @@ class boundary_normals(Enum): new_mesh.boundary_normals = boundary_normals + # Full annulus with spokes: rigid rotation about z-axis + x, y = new_mesh.X + new_mesh._nullspace_rotations = [sympy.Matrix([-y, x])] + return new_mesh @@ -1392,6 +1400,10 @@ class boundary_normals(Enum): new_mesh.boundary_normals = boundary_normals + # Full annulus with internal boundary: rigid rotation about z-axis + x, y = new_mesh.X + new_mesh._nullspace_rotations = [sympy.Matrix([-y, x])] + return new_mesh @@ -1671,4 +1683,8 @@ class boundary_normals(Enum): new_mesh.boundary_normals = boundary_normals + # Full disc with internal boundaries: rigid rotation about z-axis + x, y = new_mesh.X + new_mesh._nullspace_rotations = [sympy.Matrix([-y, x])] + return new_mesh diff --git a/src/underworld3/meshing/segmented.py b/src/underworld3/meshing/segmented.py index 4c555191f..b5d264415 100644 --- a/src/underworld3/meshing/segmented.py +++ b/src/underworld3/meshing/segmented.py @@ -676,6 +676,14 @@ class boundary_normals(Enum): new_mesh.boundary_normals = boundary_normals + # Full segmented spherical shell: 3 rigid rotation modes + x, y, z = new_mesh.X + new_mesh._nullspace_rotations = [ + sympy.Matrix([0, -z, y]), + sympy.Matrix([z, 0, -x]), + sympy.Matrix([-y, x, 0]), + ] + return new_mesh @@ -1088,4 +1096,12 @@ class boundary_normals(Enum): new_mesh.boundary_normals = boundary_normals + # Solid sphere: 3 rigid rotation modes + x, y, z = new_mesh.X + new_mesh._nullspace_rotations = [ + sympy.Matrix([0, -z, y]), + sympy.Matrix([z, 0, -x]), + sympy.Matrix([-y, x, 0]), + ] + return new_mesh diff --git a/src/underworld3/meshing/spherical.py b/src/underworld3/meshing/spherical.py index 732705d86..7d74116db 100644 --- a/src/underworld3/meshing/spherical.py +++ b/src/underworld3/meshing/spherical.py @@ -269,6 +269,14 @@ class boundary_normals(Enum): Upper = 12 Centre = 1 + # Full spherical shell: 3 rigid rotation modes + x, y, z = new_mesh.X + new_mesh._nullspace_rotations = [ + sympy.Matrix([0, -z, y]), # rotation about x + sympy.Matrix([z, 0, -x]), # rotation about y + sympy.Matrix([-y, x, 0]), # rotation about z + ] + return new_mesh @@ -475,6 +483,14 @@ class boundary_normals(Enum): Upper = 13 Centre = 1 + # Full spherical shell with internal boundary: 3 rigid rotation modes + x, y, z = new_mesh.X + new_mesh._nullspace_rotations = [ + sympy.Matrix([0, -z, y]), + sympy.Matrix([z, 0, -x]), + sympy.Matrix([-y, x, 0]), + ] + return new_mesh @@ -1011,4 +1027,12 @@ class boundary_normals(Enum): new_mesh.boundary_normals = boundary_normals + # Full cubed sphere: 3 rigid rotation modes + x, y, z = new_mesh.X + new_mesh._nullspace_rotations = [ + sympy.Matrix([0, -z, y]), + sympy.Matrix([z, 0, -x]), + sympy.Matrix([-y, x, 0]), + ] + return new_mesh From b43f8941b56188e6e610efb6bcf394d2e157797a Mon Sep 17 00:00:00 2001 From: Louis Moresi Date: Wed, 25 Mar 2026 17:20:15 +1100 Subject: [PATCH 2/3] Add petsc_use_nullspace convenience property to Stokes solver Setting stokes.petsc_use_nullspace = True enables: 1. Constant-pressure nullspace (petsc_use_pressure_nullspace) 2. Velocity rotation modes from mesh.nullspace_rotations This auto-populates petsc_velocity_nullspace_basis from the mesh geometry, so users don't need to construct rotation modes manually. For annulus: 1 pressure + 1 rotation = 2 modes For spherical shell: 1 pressure + 3 rotations = 4 modes For box: 1 pressure + 0 rotations = 1 mode Existing PR #91 tests pass unchanged. Underworld development team with AI support from Claude Code --- .../cython/petsc_generic_snes_solvers.pyx | 26 +++++++++++++++++++ 1 file changed, 26 insertions(+) diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index 62c411630..6e0ef4a62 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -3553,6 +3553,32 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): self.is_setup = False self._saddle_preconditioner = function + @property + def petsc_use_nullspace(self): + """Enable full nullspace handling: constant pressure + mesh rotation modes. + + Convenience property that enables the pressure nullspace and + auto-populates velocity nullspace modes from ``mesh.nullspace_rotations``. + + For finer control, use ``petsc_use_pressure_nullspace`` and + ``petsc_velocity_nullspace_basis`` separately. + """ + return (self._petsc_use_pressure_nullspace + and len(self._petsc_velocity_nullspace_basis) > 0) + + @petsc_use_nullspace.setter + def petsc_use_nullspace(self, value): + value = bool(value) + self._petsc_use_pressure_nullspace = value + if value and hasattr(self.mesh, 'nullspace_rotations'): + modes = self.mesh.nullspace_rotations + if modes: + self.petsc_velocity_nullspace_basis = modes + elif not value: + self._petsc_velocity_nullspace_basis = () + self._reset_stokes_nullspace() + self.is_setup = False + @property def petsc_use_pressure_nullspace(self): """ From 367219937b12a9cb844e773e3128f7d8a4f398e5 Mon Sep 17 00:00:00 2001 From: Louis Moresi Date: Thu, 2 Apr 2026 22:44:58 +1100 Subject: [PATCH 3/3] Fix petsc_use_nullspace getter and setter logic - Getter: use 'or' instead of 'and' so pressure-only nullspace (no rotation modes) reports True correctly - Setter: always sync velocity basis to mesh.nullspace_rotations when enabling (even if empty) to prevent stale modes from a previous mesh carrying over Addresses Copilot review comments on PR #105. Underworld development team with AI support from Claude Code --- .../cython/petsc_generic_snes_solvers.pyx | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index 6e0ef4a62..8b8735956 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -3564,17 +3564,20 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): ``petsc_velocity_nullspace_basis`` separately. """ return (self._petsc_use_pressure_nullspace - and len(self._petsc_velocity_nullspace_basis) > 0) + or len(self._petsc_velocity_nullspace_basis) > 0) @petsc_use_nullspace.setter def petsc_use_nullspace(self, value): value = bool(value) self._petsc_use_pressure_nullspace = value - if value and hasattr(self.mesh, 'nullspace_rotations'): - modes = self.mesh.nullspace_rotations - if modes: - self.petsc_velocity_nullspace_basis = modes - elif not value: + if value: + # Sync velocity nullspace basis to mesh rotations (even if empty) + # to avoid stale modes carrying over between meshes/runs + if hasattr(self.mesh, 'nullspace_rotations'): + self.petsc_velocity_nullspace_basis = self.mesh.nullspace_rotations or () + else: + self._petsc_velocity_nullspace_basis = () + else: self._petsc_velocity_nullspace_basis = () self._reset_stokes_nullspace() self.is_setup = False