From e956fccdf728d739b2bf853408291af8deda6953 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Sun, 16 Aug 2026 22:20:14 +1000 Subject: [PATCH 1/2] Refuse a rotated constraint and a block constraint on one solver; drop dead label lookups Two independent fixes in the same boundary-condition code. #464. add_rotated_freeslip_bc, add_fault_bc and add_constraint_bc could all be called on the same solver and nothing checked. The rotated driver builds its own index-set fieldsplit over velocity and pressure and build_rotation addresses them by field number, so a block constraint's multiplier DOFs are in neither index set and the preconditioner is built over a strict subset of the operator's rows. The two impose the same wall-normal condition by different means, so _reject_mixed_constraint_mechanisms refuses the combination at whichever call comes second, naming both. test_0062 covers both orderings and, as controls, that several rotated boundaries and several block constraints are each still allowed. With the guard neutered the solver accepts one rotated boundary condition and one multiplier together and says nothing, which is what the issue reports. #506. bc_is = bc_label.getStratumIS(value) was assigned at four sites in petsc_generic_snes_solvers.pyx and read at none; bc_label at those sites existed only to produce it. A stratum IS beside a collective DS registration invites a rank-local skip, which is why the issue asks for it to go. Closes #464, #506. Underworld development team with AI support from Claude Code --- .../cython/petsc_generic_snes_solvers.pyx | 61 ++++++++++++++--- src/underworld3/systems/solvers.py | 2 + ...est_0062_rotated_constraint_exclusivity.py | 65 +++++++++++++++++++ 3 files changed, 119 insertions(+), 9 deletions(-) create mode 100644 tests/test_0062_rotated_constraint_exclusivity.py diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index 3ff11af0c..4d3b3386d 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -3545,8 +3545,6 @@ class SNES_Scalar(SolverBaseClass): value = mesh.boundaries[bc.boundary].value ind = value - bc_label = self.dm.getLabel(boundary) - bc_is = bc_label.getStratumIS(value) self.natural_bcs[index] = self.natural_bcs[index]._replace(boundary_label_val=value) # use type 5 bc for `DM_BC_ESSENTIAL_FIELD` enum @@ -3685,9 +3683,6 @@ class SNES_Scalar(SolverBaseClass): boundary = bc.boundary value = mesh.boundaries[bc.boundary].value - bc_label = mesh.dm.getLabel(boundary) - bc_is = bc_label.getStratumIS(value) - if bc.fn_f is not None: bd_F0 = sympy.Array(bc.fn_f) @@ -4490,8 +4485,6 @@ class SNES_Vector(SolverBaseClass): value = mesh.boundaries[bc.boundary].value ind = value - bc_label = self.dm.getLabel(boundary) - bc_is = bc_label.getStratumIS(value) self.natural_bcs[index] = self.natural_bcs[index]._replace(boundary_label_val=value) # use type 5 bc for `DM_BC_ESSENTIAL_FIELD` enum @@ -5297,8 +5290,6 @@ class SNES_MultiComponent(SolverBaseClass): value = mesh.boundaries[bc.boundary].value ind = value - bc_label = self.dm.getLabel(boundary) - bc_is = bc_label.getStratumIS(value) self.natural_bcs[index] = self.natural_bcs[index]._replace(boundary_label_val=value) bc_type = 6 @@ -6108,6 +6099,54 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): # this attrib records if we need to re-setup self.is_setup = False + def _reject_mixed_constraint_mechanisms(self, adding): + """Refuse a rotated constraint and a block constraint on one solver (#464). + + The rotated driver builds its own index-set fieldsplit over exactly two + fields (``rotated_bc._solve_rotated_iterative``), and ``build_rotation`` + addresses velocity and pressure by field number. A block constraint + registers a multiplier field of its own, and those DOFs are in neither + index set — the preconditioner would then be built over a strict subset + of the operator's rows, with nothing said about it. + + The two are alternative ways to impose the same wall-normal condition, + so asking for both is a configuration error rather than a case to + support. Supporting it would need a third split for the multipliers and + a ``build_rotation`` that knows about them. + + Parameters + ---------- + adding : str + Name of the method being called, so the message can say which + mechanism is already in place and which one was refused. + """ + + rotated = list(getattr(self, "_rotated_freeslip_bcs", None) or []) + list( + getattr(self, "_fault_contact_faults", None) or [] + ) + multipliers = list(getattr(self, "_multipliers", None) or []) + + if adding == "add_constraint_bc": + if not rotated: + return + present, refused = len(rotated), "a block constraint" + present_kind = "rotated (free-slip or fault contact)" + else: + if not multipliers: + return + present, refused = len(multipliers), "a rotated constraint" + present_kind = "block-constraint multiplier" + + raise RuntimeError( + f"{adding}(): this solver already carries {present} " + f"{present_kind} boundary condition(s), so it cannot also take " + f"{refused}. The rotated solve splits velocity and pressure by " + f"field number and a block constraint adds a multiplier field " + f"outside that split, so the preconditioner would cover only part " + f"of the operator. Both impose the same wall-normal condition — " + f"use one of them (issue #464)." + ) + def add_rotated_freeslip_bc(self, conds=None, boundary=None, normal=None): r"""Add STRONG free-slip (:math:`\mathbf{u}\cdot\hat{\mathbf n}=0`) by rotating the boundary velocity DOFs into a per-node (normal, tangential) frame and @@ -6163,6 +6202,8 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): DeprecationWarning: the string becomes ``boundary`` and a second positional argument, if present, becomes ``normal``. """ + self._reject_mixed_constraint_mechanisms("add_rotated_freeslip_bc") + if isinstance(conds, str): # legacy boundary-first call: (boundary[, normal]) if boundary is not None: @@ -6241,6 +6282,8 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): """ from underworld3.utilities import fault_contact + self._reject_mixed_constraint_mechanisms("add_fault_bc") + if not isinstance(boundary, str): raise TypeError( f"add_fault_bc() requires the fault's boundary name string; " diff --git a/src/underworld3/systems/solvers.py b/src/underworld3/systems/solvers.py index 77e1aecdc..98d0c1b7f 100644 --- a/src/underworld3/systems/solvers.py +++ b/src/underworld3/systems/solvers.py @@ -2659,6 +2659,8 @@ def add_constraint_bc(self, conds=None, boundary=None, normal=None, screening=No DeprecationWarning; see ``SolverBaseClass._value_first_bc_args``. """ + self._reject_mixed_constraint_mechanisms("add_constraint_bc") + conds, boundary = self._value_first_bc_args( "add_constraint_bc", conds, boundary, alias=g) g = conds if conds is not None else 0.0 diff --git a/tests/test_0062_rotated_constraint_exclusivity.py b/tests/test_0062_rotated_constraint_exclusivity.py new file mode 100644 index 000000000..bd713313c --- /dev/null +++ b/tests/test_0062_rotated_constraint_exclusivity.py @@ -0,0 +1,65 @@ +"""A solver takes rotated constraints or block constraints, not both (#464). + +The rotated driver builds its own index-set fieldsplit over velocity and +pressure and addresses them by field number. A block constraint registers a +multiplier field, whose DOFs are then in neither index set, so the +preconditioner covers a strict subset of the operator's rows and says nothing +about it. Both mechanisms impose the same wall-normal condition. + +The two controls matter as much as the two refusals: a guard that also refused +each mechanism on its own would pass a test that only checked that mixing +raises. Neutering `_reject_mixed_constraint_mechanisms` and repeating the first +case shows what the guard is for — the solver accepts one rotated boundary +condition and one multiplier together, and says nothing. +""" + +import pytest + +import underworld3 as uw + +pytestmark = [pytest.mark.level_1, pytest.mark.tier_a] + + +@pytest.fixture +def stokes(): + """A constrained Stokes solver, which is the class that has both APIs.""" + + mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=0.5, qdegree=2 + ) + v = uw.discretisation.MeshVariable("Uc", mesh, mesh.dim, degree=2) + p = uw.discretisation.MeshVariable("Pc", mesh, 1, degree=1) + + return uw.systems.Stokes_Constrained(mesh, velocityField=v, pressureField=p) + + +def test_block_constraint_after_rotated_freeslip_is_refused(stokes): + stokes.add_rotated_freeslip_bc(0.0, "Left") + + with pytest.raises(RuntimeError, match=r"add_constraint_bc.*#464"): + stokes.add_constraint_bc(0.0, "Bottom") + + +def test_rotated_freeslip_after_block_constraint_is_refused(stokes): + stokes.add_constraint_bc(0.0, "Bottom") + + with pytest.raises(RuntimeError, match=r"add_rotated_freeslip_bc.*#464"): + stokes.add_rotated_freeslip_bc(0.0, "Left") + + +def test_several_rotated_boundaries_are_still_allowed(stokes): + """The control: the guard is about mixing mechanisms, not about counting.""" + + stokes.add_rotated_freeslip_bc(0.0, "Left") + stokes.add_rotated_freeslip_bc(0.0, "Right") + + assert len(stokes._rotated_freeslip_bcs) == 2 + + +def test_several_block_constraints_are_still_allowed(stokes): + """The other control.""" + + stokes.add_constraint_bc(0.0, "Bottom") + stokes.add_constraint_bc(0.0, "Top") + + assert len(stokes._multipliers) == 2 From aa2b8d875cf7f196cd21f171c963e0fff8032886 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Sun, 16 Aug 2026 22:25:46 +1000 Subject: [PATCH 2/2] Check the constraint pair at the solve as well as at registration The registration checks only cover what goes through the solver's own methods, and fault_contact writes _fault_contact_faults directly rather than through add_fault_bc. The dispatch is where both lists are read together, so it carries the check as well: registration gives the better message, the dispatch gives the guarantee. It runs as the first statement of solve(), before any setup touches either list, so an unsupported pair costs nothing before it is refused. Underworld development team with AI support from Claude Code --- .../cython/petsc_generic_snes_solvers.pyx | 23 +++++++++++++++++++ ...est_0062_rotated_constraint_exclusivity.py | 17 ++++++++++++++ 2 files changed, 40 insertions(+) diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index 4d3b3386d..dcadfe232 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -6126,6 +6126,21 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): ) multipliers = list(getattr(self, "_multipliers", None) or []) + if adding == "solve": + # The dispatch reads both lists, so it can only report the pair. + if not (rotated and multipliers): + return + raise RuntimeError( + f"solve(): this solver carries {len(rotated)} rotated " + f"(free-slip or fault contact) and {len(multipliers)} " + f"block-constraint boundary condition(s). The rotated solve " + f"splits velocity and pressure by field number and the " + f"multiplier fields lie outside that split, so the " + f"preconditioner would cover only part of the operator. Both " + f"impose the same wall-normal condition — use one of them " + f"(issue #464)." + ) + if adding == "add_constraint_bc": if not rotated: return @@ -9387,6 +9402,14 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): """ + # Checked here as well as at registration (#464). The registration check + # gives the better message — it knows which call was refused — but it + # only covers what goes through the solver's own methods, and + # `fault_contact` writes `_fault_contact_faults` directly. This runs + # before any setup reads either list, so an unsupported pair costs + # nothing before it is refused. + self._reject_mixed_constraint_mechanisms("solve") + if homotopy: # The march runs a SEQUENCE of ordinary solves at successively sharper # yield surfaces; each one re-enters this method with homotopy=False. diff --git a/tests/test_0062_rotated_constraint_exclusivity.py b/tests/test_0062_rotated_constraint_exclusivity.py index bd713313c..0d75fabfb 100644 --- a/tests/test_0062_rotated_constraint_exclusivity.py +++ b/tests/test_0062_rotated_constraint_exclusivity.py @@ -63,3 +63,20 @@ def test_several_block_constraints_are_still_allowed(stokes): stokes.add_constraint_bc(0.0, "Top") assert len(stokes._multipliers) == 2 + + +def test_the_dispatch_refuses_the_pair_even_if_registration_was_bypassed(stokes): + """The guarantee, as opposed to the message. + + The registration checks only cover what goes through the solver's own + methods, and `fault_contact` writes `_fault_contact_faults` directly. The + solve dispatch is where both lists are read together, so it carries the + check as well. Reaching it here means bypassing registration, which is + exactly the case the dispatch exists to catch. + """ + + stokes.add_rotated_freeslip_bc(0.0, "Left") + stokes._multipliers.append(object()) + + with pytest.raises(RuntimeError, match=r"solve\(\).*#464"): + stokes.solve()