From 1469d92483789ffcab48de25f15977fd3b811c4d Mon Sep 17 00:00:00 2001 From: Tyagi Date: Tue, 24 Mar 2026 17:56:50 +1100 Subject: [PATCH 1/2] Add PETSc pressure nullspace support for Stokes --- .../cython/petsc_generic_snes_solvers.pyx | 107 ++++++++++++++++++ tests/test_1013_stokes_pressure_nullspace.py | 63 +++++++++++ 2 files changed, 170 insertions(+) create mode 100644 tests/test_1013_stokes_pressure_nullspace.py diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index 1f9cebea2..d3c1dd06b 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -447,6 +447,10 @@ class SolverBaseClass(uw_object): self.dm = None # Should be able to avoid nuking this if we # can insert new functions in template (surface integrals problematic in # the current implementation ) + if hasattr(self, "_pressure_nullspace"): + self._pressure_nullspace = None + if hasattr(self, "_pressure_nullspace_basis"): + self._pressure_nullspace_basis = None # This is a workaround for some problem in the PETSc machinery # where we need a surface integral term somewhere on every process @@ -2829,6 +2833,9 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): self.boundary_conditions = False # self._constitutive_model = None self._saddle_preconditioner = None + self._petsc_use_pressure_nullspace = False + self._pressure_nullspace = None + self._pressure_nullspace_basis = None # Construct strainrate tensor for future usage. # Grab gradients, and let's switch out to sympy.Matrix notation @@ -3088,6 +3095,100 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): self.is_setup = False self._saddle_preconditioner = function + @property + def petsc_use_pressure_nullspace(self): + """ + Enable PETSc handling of the constant-pressure nullspace. + + When enabled, the solver attaches a nullspace basis with zero + velocity entries and constant pressure entries to the coupled + Stokes matrices before solve. + """ + return self._petsc_use_pressure_nullspace + + @petsc_use_pressure_nullspace.setter + def petsc_use_pressure_nullspace(self, value): + self._petsc_use_pressure_nullspace = bool(value) + self._pressure_nullspace = None + self._pressure_nullspace_basis = None + self.is_setup = False + + def _pressure_dirichlet_bcs(self): + """Return essential boundary conditions applied to the pressure field.""" + + pressure_field_ids = {1} + try: + pressure_field_ids.add(self.p.field_id) + except Exception: + pass + + return [bc for bc in self.essential_bcs if bc.f_id in pressure_field_ids] + + def _build_pressure_nullspace(self): + """Create a constant-pressure nullspace basis for the coupled Stokes DM.""" + + template_vec = self.dm.getGlobalVec() + try: + null_vec = template_vec.duplicate() + finally: + self.dm.restoreGlobalVec(template_vec) + + null_vec.set(0.0) + + pressure_is = self._subdict["pressure"][0] + pressure_subvec = null_vec.getSubVector(pressure_is) + pressure_subvec.set(1.0) + null_vec.restoreSubVector(pressure_is, pressure_subvec) + + self._pressure_nullspace_basis = null_vec + self._pressure_nullspace = PETSc.NullSpace().create( + constant=False, + vectors=(null_vec,), + comm=self.dm.comm, + ) + + return self._pressure_nullspace + + def _attach_pressure_nullspace(self): + """Attach the configured pressure nullspace to the Stokes matrices.""" + + if not self._petsc_use_pressure_nullspace: + return + + pressure_bcs = self._pressure_dirichlet_bcs() + if pressure_bcs: + boundaries = ", ".join(sorted({bc.boundary for bc in pressure_bcs})) + raise ValueError( + "petsc_use_pressure_nullspace=True requires the pressure field to be " + f"free of Dirichlet boundary conditions. Found pressure Dirichlet BCs on: {boundaries}" + ) + + if "pressure" not in self._subdict: + raise RuntimeError("Pressure field decomposition is unavailable; cannot attach nullspace.") + + self.snes.setUp() + + jacobian = self.snes.getJacobian() + operator_matrix = jacobian[0] + preconditioner_matrix = jacobian[1] if len(jacobian) > 1 else None + + nullspace = self._pressure_nullspace + if nullspace is None: + nullspace = self._build_pressure_nullspace() + + operator_matrix.setNullSpace(nullspace) + operator_matrix.setTransposeNullSpace(nullspace) + + if preconditioner_matrix is not None: + preconditioner_matrix.setNullSpace(nullspace) + preconditioner_matrix.setTransposeNullSpace(nullspace) + + if self.verbose and uw.mpi.rank == 0: + print( + f"Stokes Saddle Pt ({self.name}): attached constant-pressure nullspace", + flush=True, + ) + ## F0, F1 should be f0 and F1, (pf0 for Saddles can be added here) ## don't add new ones uf0, uF1 are redundant @@ -3744,6 +3845,8 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): for index,name in enumerate(names): self._subdict[name] = (isets[index],dms[index]) + self._attach_pressure_nullspace() + self.is_setup = True self.constitutive_model._solver_is_setup = True @@ -3846,6 +3949,7 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): self.petsc_options.setValue("snes_max_it", 0) self.snes.setType("nrichardson") self.snes.setFromOptions() + self._attach_pressure_nullspace() self.snes.solve(None, gvec) # with self.mesh.access(): @@ -3871,6 +3975,7 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): self.snes.atol = self.atol self.snes.setType("nrichardson") self.snes.setFromOptions() + self._attach_pressure_nullspace() self.snes.solve(None, gvec) self._warn_on_divergence(phase="picard") @@ -3880,6 +3985,7 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): self.snes.atol = self.atol self.petsc_options.setValue("snes_max_it", snes_max_it) self.snes.setFromOptions() + self._attach_pressure_nullspace() self.snes.solve(None, gvec) else: @@ -3889,6 +3995,7 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): self.snes.atol = self.atol self.petsc_options.setValue("snes_max_it", snes_max_it) self.snes.setFromOptions() + self._attach_pressure_nullspace() self.snes.solve(None, gvec) cdef DM dm = self.dm diff --git a/tests/test_1013_stokes_pressure_nullspace.py b/tests/test_1013_stokes_pressure_nullspace.py new file mode 100644 index 000000000..0be4f41fd --- /dev/null +++ b/tests/test_1013_stokes_pressure_nullspace.py @@ -0,0 +1,63 @@ +import pytest + +pytestmark = pytest.mark.level_3 + +import sympy +import underworld3 as uw + + +def test_stokes_pressure_nullspace_solves_without_pressure_bc(): + mesh = uw.meshing.StructuredQuadBox(elementRes=(3, 3)) + x, y = mesh.X + + u = uw.discretisation.MeshVariable( + "u_nullspace", + mesh, + mesh.dim, + vtype=uw.VarType.VECTOR, + degree=2, + ) + p = uw.discretisation.MeshVariable( + "p_nullspace", + mesh, + 1, + vtype=uw.VarType.SCALAR, + degree=1, + continuous=True, + ) + + stokes = uw.systems.Stokes(mesh, velocityField=u, pressureField=p) + stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel + stokes.constitutive_model.Parameters.shear_viscosity_0 = 1.0 + stokes.bodyforce = sympy.Matrix([0.0, x]) + stokes.petsc_use_pressure_nullspace = True + + stokes.tolerance = 1.0e-4 + stokes.petsc_options["snes_type"] = "ksponly" + stokes.petsc_options["ksp_type"] = "fgmres" + stokes.petsc_options["ksp_rtol"] = 1.0e-4 + stokes.petsc_options["ksp_atol"] = 0.0 + + stokes.petsc_options.setValue("fieldsplit_velocity_pc_mg_type", "kaskade") + stokes.petsc_options.setValue("fieldsplit_velocity_pc_mg_cycle_type", "w") + stokes.petsc_options["fieldsplit_velocity_mg_coarse_pc_type"] = "svd" + stokes.petsc_options["fieldsplit_velocity_ksp_type"] = "fcg" + stokes.petsc_options["fieldsplit_velocity_mg_levels_ksp_type"] = "chebyshev" + stokes.petsc_options["fieldsplit_velocity_mg_levels_ksp_max_it"] = 5 + stokes.petsc_options["fieldsplit_velocity_mg_levels_ksp_converged_maxits"] = None + stokes.petsc_options.setValue("fieldsplit_pressure_pc_type", "mg") + stokes.petsc_options.setValue("fieldsplit_pressure_pc_mg_type", "multiplicative") + stokes.petsc_options.setValue("fieldsplit_pressure_pc_mg_cycle_type", "v") + + stokes.add_dirichlet_bc((0.0, 0.0), "Bottom") + stokes.add_dirichlet_bc((0.0, None), "Top") + stokes.add_dirichlet_bc((0.0, None), "Left") + stokes.add_condition(conds=(0.0, None), label="Right", f_id=0, c_type="dirichlet") + + stokes.solve() + + assert stokes.snes.getConvergedReason() > 0 + + jacobian = stokes.snes.getJacobian() + nullspace = jacobian[0].getNullSpace() + assert nullspace is not None From 5f6d09309382b4f1ea909cd42f63626d05ca9ed4 Mon Sep 17 00:00:00 2001 From: Tyagi Date: Wed, 25 Mar 2026 01:24:26 +1100 Subject: [PATCH 2/2] Generalize Stokes nullspace support for shell modes --- .../cython/petsc_generic_snes_solvers.pyx | 243 +++++++++++++++--- tests/test_1013_stokes_pressure_nullspace.py | 15 ++ tests/test_1014_stokes_shell_nullspace.py | 129 ++++++++++ 3 files changed, 350 insertions(+), 37 deletions(-) create mode 100644 tests/test_1014_stokes_shell_nullspace.py diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index d3c1dd06b..34cd3a328 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -1,5 +1,6 @@ from xmlrpc.client import Boolean +import numpy as np import sympy from sympy import sympify @@ -447,10 +448,10 @@ class SolverBaseClass(uw_object): self.dm = None # Should be able to avoid nuking this if we # can insert new functions in template (surface integrals problematic in # the current implementation ) - if hasattr(self, "_pressure_nullspace"): - self._pressure_nullspace = None - if hasattr(self, "_pressure_nullspace_basis"): - self._pressure_nullspace_basis = None + if hasattr(self, "_stokes_nullspace"): + self._stokes_nullspace = None + if hasattr(self, "_stokes_nullspace_basis"): + self._stokes_nullspace_basis = () # This is a workaround for some problem in the PETSc machinery # where we need a surface integral term somewhere on every process @@ -2834,8 +2835,9 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): # self._constitutive_model = None self._saddle_preconditioner = None self._petsc_use_pressure_nullspace = False - self._pressure_nullspace = None - self._pressure_nullspace_basis = None + self._petsc_velocity_nullspace_basis = () + self._stokes_nullspace = None + self._stokes_nullspace_basis = () # Construct strainrate tensor for future usage. # Grab gradients, and let's switch out to sympy.Matrix notation @@ -3100,32 +3102,128 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): """ Enable PETSc handling of the constant-pressure nullspace. - When enabled, the solver attaches a nullspace basis with zero - velocity entries and constant pressure entries to the coupled - Stokes matrices before solve. + When enabled, the solver attaches the constant-pressure mode to + the coupled Stokes nullspace basis before solve. Additional + user-supplied velocity nullspace modes can be configured through + ``petsc_velocity_nullspace_basis``. + + For free-slip shell problems this is typically used together with + the rigid-body rotation modes documented on + ``petsc_velocity_nullspace_basis``. + + Examples + -------- + 2-D annulus with pressure gauge only: + >>> stokes.petsc_use_pressure_nullspace = True + + 2-D annulus with pressure plus rigid rotation: + >>> x, y = mesh.X + >>> stokes.petsc_use_pressure_nullspace = True + >>> stokes.petsc_velocity_nullspace_basis = [sympy.Matrix([-y, x])] + + 3-D spherical shell with pressure plus the three rigid rotations: + >>> x, y, z = mesh.X + >>> stokes.petsc_use_pressure_nullspace = True + >>> stokes.petsc_velocity_nullspace_basis = [ + ... sympy.Matrix([0, -z, y]), + ... sympy.Matrix([z, 0, -x]), + ... sympy.Matrix([-y, x, 0]), + ... ] """ return self._petsc_use_pressure_nullspace @petsc_use_pressure_nullspace.setter def petsc_use_pressure_nullspace(self, value): self._petsc_use_pressure_nullspace = bool(value) - self._pressure_nullspace = None - self._pressure_nullspace_basis = None + self._reset_stokes_nullspace() self.is_setup = False + @property + def petsc_velocity_nullspace_basis(self): + """ + Optional exact velocity nullspace modes for the coupled Stokes solve. + + Each entry must be a vector-valued SymPy expression defined in the + mesh coordinate system and representing an exact null mode of the + configured Stokes operator. Typical examples are rigid-body rotation + modes for annulus or spherical-shell free-slip problems. + + For centered shell geometries with exact free-slip / no-penetration + boundary conditions, the rigid-body rotation modes are: + + - 2-D annulus: one mode, ``(-y, x)``, equivalent to ``r e_theta`` + - 3-D spherical shell: three modes, + ``(0, -z, y)``, ``(z, 0, -x)``, and ``(-y, x, 0)`` + + These are the velocity fields generated by rigid rotations + ``u = omega x x``. They are tangent to concentric circles / spheres + and have zero strain rate, so they are exact velocity null modes for + the free-slip shell Stokes operator. + + They are not exact null modes when the boundary conditions select a + specific tangential velocity, for example: + + - essential velocity boundary conditions + - penalty boundary conditions on the full velocity error + ``u - u_analytic`` + + To remove shell nullspaces in Stokes, set the pressure mode and then + provide the exact rotation basis: + + Examples + -------- + 2-D annulus: + >>> x, y = mesh.X + >>> stokes.petsc_use_pressure_nullspace = True + >>> stokes.petsc_velocity_nullspace_basis = [sympy.Matrix([-y, x])] + + 3-D spherical shell: + >>> x, y, z = mesh.X + >>> stokes.petsc_use_pressure_nullspace = True + >>> stokes.petsc_velocity_nullspace_basis = [ + ... sympy.Matrix([0, -z, y]), + ... sympy.Matrix([z, 0, -x]), + ... sympy.Matrix([-y, x, 0]), + ... ] + """ + return self._petsc_velocity_nullspace_basis + + @petsc_velocity_nullspace_basis.setter + def petsc_velocity_nullspace_basis(self, modes): + if modes is None: + modes = () + + velocity_modes = [] + for mode in modes: + matrix_mode = sympy.Matrix(mode) + if matrix_mode.shape == (1, self.mesh.dim): + matrix_mode = matrix_mode.T + if matrix_mode.shape != (self.mesh.dim, 1): + raise ValueError( + "Each petsc_velocity_nullspace_basis mode must have shape " + f"({self.mesh.dim}, 1) or (1, {self.mesh.dim}); got {matrix_mode.shape}." + ) + velocity_modes.append(matrix_mode) + + self._petsc_velocity_nullspace_basis = tuple(velocity_modes) + self._reset_stokes_nullspace() + self.is_setup = False + + def _reset_stokes_nullspace(self): + self._stokes_nullspace = None + self._stokes_nullspace_basis = () + def _pressure_dirichlet_bcs(self): """Return essential boundary conditions applied to the pressure field.""" - pressure_field_ids = {1} - try: - pressure_field_ids.add(self.p.field_id) - except Exception: - pass + pressure_field_id = getattr(getattr(self, "p", None), "field_id", None) + if pressure_field_id is None: + raise RuntimeError("Pressure field is unavailable; cannot inspect pressure Dirichlet BCs.") - return [bc for bc in self.essential_bcs if bc.f_id in pressure_field_ids] + return [bc for bc in self.essential_bcs if bc.f_id == pressure_field_id] - def _build_pressure_nullspace(self): - """Create a constant-pressure nullspace basis for the coupled Stokes DM.""" + def _build_pressure_nullspace_vector(self): + """Create the constant-pressure basis vector for the coupled Stokes DM.""" template_vec = self.dm.getGlobalVec() try: @@ -3140,31 +3238,96 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): pressure_subvec.set(1.0) null_vec.restoreSubVector(pressure_is, pressure_subvec) - self._pressure_nullspace_basis = null_vec - self._pressure_nullspace = PETSc.NullSpace().create( + return null_vec + + def _build_velocity_nullspace_vector(self, mode): + """Create a velocity nullspace basis vector from a user-supplied mode.""" + + mode_values = np.asarray(uw.function.evaluate(mode, self.u.coords_nd), dtype=np.float64) + if mode_values.ndim == 1: + mode_values = mode_values.reshape(-1, 1) + elif mode_values.ndim > 2: + mode_values = mode_values.reshape(mode_values.shape[0], -1) + + if mode_values.shape != (self.u.coords_nd.shape[0], self.u.num_components): + raise ValueError( + "Velocity nullspace mode evaluation must return an array with shape " + f"({self.u.coords_nd.shape[0]}, {self.u.num_components}); got {mode_values.shape}." + ) + + template_vec = self.dm.getGlobalVec() + try: + null_vec = template_vec.duplicate() + finally: + self.dm.restoreGlobalVec(template_vec) + + null_vec.set(0.0) + + velocity_is, velocity_subdm = self._subdict["velocity"] + velocity_basis = self.u.vec.duplicate() + try: + velocity_basis.array[:] = mode_values.reshape(velocity_basis.array.shape) + velocity_subvec = null_vec.getSubVector(velocity_is) + velocity_subdm.localToGlobal(velocity_basis, velocity_subvec, addv=False) + null_vec.restoreSubVector(velocity_is, velocity_subvec) + finally: + velocity_basis.destroy() + + return null_vec + + def _build_stokes_nullspace(self): + """Create the configured coupled Stokes nullspace basis.""" + + basis_vectors = [] + + if self._petsc_use_pressure_nullspace: + basis_vectors.append(self._build_pressure_nullspace_vector()) + + for mode in self._petsc_velocity_nullspace_basis: + basis_vectors.append(self._build_velocity_nullspace_vector(mode)) + + if not basis_vectors: + return None + + orthonormal_basis = [] + for basis_vec in basis_vectors: + for orth_vec in orthonormal_basis: + basis_vec.axpy(-orth_vec.dot(basis_vec), orth_vec) + + basis_norm = basis_vec.norm() + if np.isclose(basis_norm, 0.0): + raise ValueError( + "Configured Stokes nullspace basis contains a dependent or zero mode." + ) + + basis_vec.scale(1.0 / basis_norm) + orthonormal_basis.append(basis_vec) + + self._stokes_nullspace_basis = tuple(orthonormal_basis) + self._stokes_nullspace = PETSc.NullSpace().create( constant=False, - vectors=(null_vec,), + vectors=self._stokes_nullspace_basis, comm=self.dm.comm, ) - return self._pressure_nullspace + return self._stokes_nullspace - def _attach_pressure_nullspace(self): - """Attach the configured pressure nullspace to the Stokes matrices.""" + def _attach_stokes_nullspace(self): + """Attach the configured coupled Stokes nullspace to the solver matrices.""" - if not self._petsc_use_pressure_nullspace: + if not self._petsc_use_pressure_nullspace and not self._petsc_velocity_nullspace_basis: return pressure_bcs = self._pressure_dirichlet_bcs() if pressure_bcs: boundaries = ", ".join(sorted({bc.boundary for bc in pressure_bcs})) raise ValueError( - "petsc_use_pressure_nullspace=True requires the pressure field to be " + "PETSc Stokes nullspace support requires the pressure field to be " f"free of Dirichlet boundary conditions. Found pressure Dirichlet BCs on: {boundaries}" ) - if "pressure" not in self._subdict: - raise RuntimeError("Pressure field decomposition is unavailable; cannot attach nullspace.") + if "pressure" not in self._subdict or "velocity" not in self._subdict: + raise RuntimeError("Velocity/pressure field decomposition is unavailable; cannot attach nullspace.") self.snes.setUp() @@ -3172,9 +3335,12 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): operator_matrix = jacobian[0] preconditioner_matrix = jacobian[1] if len(jacobian) > 1 else None - nullspace = self._pressure_nullspace + nullspace = self._stokes_nullspace if nullspace is None: - nullspace = self._build_pressure_nullspace() + nullspace = self._build_stokes_nullspace() + + if nullspace is None: + return operator_matrix.setNullSpace(nullspace) operator_matrix.setTransposeNullSpace(nullspace) @@ -3185,7 +3351,8 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): if self.verbose and uw.mpi.rank == 0: print( - f"Stokes Saddle Pt ({self.name}): attached constant-pressure nullspace", + f"Stokes Saddle Pt ({self.name}): attached Stokes nullspace with " + f"{len(self._stokes_nullspace_basis)} basis mode(s)", flush=True, ) @@ -3845,7 +4012,7 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): for index,name in enumerate(names): self._subdict[name] = (isets[index],dms[index]) - self._attach_pressure_nullspace() + self._attach_stokes_nullspace() self.is_setup = True self.constitutive_model._solver_is_setup = True @@ -3949,7 +4116,9 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): self.petsc_options.setValue("snes_max_it", 0) self.snes.setType("nrichardson") self.snes.setFromOptions() - self._attach_pressure_nullspace() + # PETSc may rebuild operator state after setFromOptions(), so reattach + # the configured Stokes nullspace before each solve path. + self._attach_stokes_nullspace() self.snes.solve(None, gvec) # with self.mesh.access(): @@ -3975,7 +4144,7 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): self.snes.atol = self.atol self.snes.setType("nrichardson") self.snes.setFromOptions() - self._attach_pressure_nullspace() + self._attach_stokes_nullspace() self.snes.solve(None, gvec) self._warn_on_divergence(phase="picard") @@ -3985,7 +4154,7 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): self.snes.atol = self.atol self.petsc_options.setValue("snes_max_it", snes_max_it) self.snes.setFromOptions() - self._attach_pressure_nullspace() + self._attach_stokes_nullspace() self.snes.solve(None, gvec) else: @@ -3995,7 +4164,7 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): self.snes.atol = self.atol self.petsc_options.setValue("snes_max_it", snes_max_it) self.snes.setFromOptions() - self._attach_pressure_nullspace() + self._attach_stokes_nullspace() self.snes.solve(None, gvec) cdef DM dm = self.dm diff --git a/tests/test_1013_stokes_pressure_nullspace.py b/tests/test_1013_stokes_pressure_nullspace.py index 0be4f41fd..ad5cbb86c 100644 --- a/tests/test_1013_stokes_pressure_nullspace.py +++ b/tests/test_1013_stokes_pressure_nullspace.py @@ -57,7 +57,22 @@ def test_stokes_pressure_nullspace_solves_without_pressure_bc(): stokes.solve() assert stokes.snes.getConvergedReason() > 0 + assert len(stokes._stokes_nullspace_basis) == 1 jacobian = stokes.snes.getJacobian() nullspace = jacobian[0].getNullSpace() assert nullspace is not None + + basis_vec = stokes._stokes_nullspace_basis[0] + velocity_is = stokes._subdict["velocity"][0] + pressure_is = stokes._subdict["pressure"][0] + + velocity_subvec = basis_vec.getSubVector(velocity_is) + pressure_subvec = basis_vec.getSubVector(pressure_is) + + try: + assert velocity_subvec.norm() == pytest.approx(0.0, abs=1.0e-12) + assert pressure_subvec.norm() > 0.0 + finally: + basis_vec.restoreSubVector(velocity_is, velocity_subvec) + basis_vec.restoreSubVector(pressure_is, pressure_subvec) diff --git a/tests/test_1014_stokes_shell_nullspace.py b/tests/test_1014_stokes_shell_nullspace.py new file mode 100644 index 000000000..51c20325e --- /dev/null +++ b/tests/test_1014_stokes_shell_nullspace.py @@ -0,0 +1,129 @@ +import pytest + +pytestmark = pytest.mark.level_3 + +import sympy +import underworld3 as uw + + +def _configure_shell_stokes(mesh): + u = uw.discretisation.MeshVariable( + "u_shell_nullspace", + mesh, + mesh.dim, + vtype=uw.VarType.VECTOR, + degree=2, + ) + p = uw.discretisation.MeshVariable( + "p_shell_nullspace", + mesh, + 1, + vtype=uw.VarType.SCALAR, + degree=1, + continuous=True, + ) + + stokes = uw.systems.Stokes(mesh, velocityField=u, pressureField=p) + stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel + stokes.constitutive_model.Parameters.shear_viscosity_0 = 1.0 + stokes.bodyforce = sympy.Matrix([0.0] * mesh.dim) + + stokes.tolerance = 1.0e-4 + stokes.petsc_options["snes_type"] = "ksponly" + stokes.petsc_options["ksp_type"] = "fgmres" + stokes.petsc_options["ksp_rtol"] = 1.0e-4 + stokes.petsc_options["ksp_atol"] = 0.0 + + stokes.petsc_options.setValue("fieldsplit_velocity_pc_mg_type", "kaskade") + stokes.petsc_options.setValue("fieldsplit_velocity_pc_mg_cycle_type", "w") + stokes.petsc_options["fieldsplit_velocity_mg_coarse_pc_type"] = "svd" + stokes.petsc_options["fieldsplit_velocity_ksp_type"] = "fcg" + stokes.petsc_options["fieldsplit_velocity_mg_levels_ksp_type"] = "chebyshev" + stokes.petsc_options["fieldsplit_velocity_mg_levels_ksp_max_it"] = 5 + stokes.petsc_options["fieldsplit_velocity_mg_levels_ksp_converged_maxits"] = None + stokes.petsc_options.setValue("fieldsplit_pressure_pc_type", "mg") + stokes.petsc_options.setValue("fieldsplit_pressure_pc_mg_type", "multiplicative") + stokes.petsc_options.setValue("fieldsplit_pressure_pc_mg_cycle_type", "v") + + gamma = mesh.Gamma + stokes.add_natural_bc(1.0e4 * gamma.dot(u) * gamma, "Upper") + stokes.add_natural_bc(1.0e4 * gamma.dot(u) * gamma, "Lower") + + stokes.petsc_use_pressure_nullspace = True + + return stokes + + +@pytest.mark.parametrize( + ("mesh", "rotation_modes", "expected_mode_count"), + [ + ( + uw.meshing.Annulus(radiusOuter=1.0, radiusInner=0.6, cellSize=0.25, qdegree=2), + [sympy.Matrix([-sympy.Symbol("y"), sympy.Symbol("x")])], + 2, + ), + ( + uw.meshing.SphericalShell(radiusOuter=1.0, radiusInner=0.6, cellSize=0.5, qdegree=2), + [ + sympy.Matrix([0, -sympy.Symbol("z"), sympy.Symbol("y")]), + sympy.Matrix([sympy.Symbol("z"), 0, -sympy.Symbol("x")]), + sympy.Matrix([-sympy.Symbol("y"), sympy.Symbol("x"), 0]), + ], + 4, + ), + ], +) +def test_stokes_shell_rotation_nullspace(mesh, rotation_modes, expected_mode_count): + if mesh.dim == 2: + x, y = mesh.X + coordinate_subs = { + sympy.Symbol("x"): x, + sympy.Symbol("y"): y, + } + else: + x, y, z = mesh.X + coordinate_subs = { + sympy.Symbol("x"): x, + sympy.Symbol("y"): y, + sympy.Symbol("z"): z, + } + + resolved_modes = [mode.subs(coordinate_subs) for mode in rotation_modes] + + stokes = _configure_shell_stokes(mesh) + stokes.petsc_velocity_nullspace_basis = resolved_modes + stokes.solve() + + assert stokes.snes.getConvergedReason() > 0 + assert len(stokes._stokes_nullspace_basis) == expected_mode_count + + jacobian = stokes.snes.getJacobian() + nullspace = jacobian[0].getNullSpace() + + assert nullspace is not None + + velocity_is = stokes._subdict["velocity"][0] + pressure_is = stokes._subdict["pressure"][0] + pressure_modes = 0 + velocity_modes = 0 + + for basis_vec in stokes._stokes_nullspace_basis: + velocity_subvec = basis_vec.getSubVector(velocity_is) + pressure_subvec = basis_vec.getSubVector(pressure_is) + + try: + velocity_norm = velocity_subvec.norm() + pressure_norm = pressure_subvec.norm() + finally: + basis_vec.restoreSubVector(velocity_is, velocity_subvec) + basis_vec.restoreSubVector(pressure_is, pressure_subvec) + + if pressure_norm > 1.0e-10: + pressure_modes += 1 + assert velocity_norm == pytest.approx(0.0, abs=1.0e-12) + else: + velocity_modes += 1 + assert velocity_norm > 0.0 + + assert pressure_modes == 1 + assert velocity_modes == expected_mode_count - 1