diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index 3ff11af0c..dcadfe232 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,69 @@ 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 == "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 + 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 +6217,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 +6297,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; " @@ -9344,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/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..0d75fabfb --- /dev/null +++ b/tests/test_0062_rotated_constraint_exclusivity.py @@ -0,0 +1,82 @@ +"""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 + + +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()