Remove the Null_Boundary fake natural BC and the label it was built on

src/underworld3/cython/petsc_generic_snes_solvers.pyx

@@ -1961,17 +1961,6 @@ class SolverBaseClass(uw_object):
             if hasattr(self, "_constant_nullspace_obj"):
                 self._constant_nullspace_obj = None
 
-        # This is a workaround for some problem in the PETSc machinery
-        # where we need a surface integral term somewhere on every process
-        # if we have a contribution from anywhere. We add a fake one here
-        # which just integrates nothing over a bunch of points. It's enough
-        # to let the rest of the machinery work.
-
-        if len(self.natural_bcs) > 0:
-            if not any(bc.boundary == "Null_Boundary" for bc in self.natural_bcs):
-                bc = (0,)*self.Unknowns.u.shape[1]
-                self.add_natural_bc(bc, "Null_Boundary")
-
         if verbose:
             uw.pprint("Build pointwise functions")
         self._setup_pointwise_functions(verbose, debug=debug, debug_name=debug_name)

src/underworld3/discretisation/discretisation_mesh.py

@@ -884,13 +884,24 @@ def _patch_boundary_enum(
     ):
         """Extend the boundary enum and rebuild the DM's boundary labels.
 
-        Patches the user-supplied boundary enum with the members every UW
-        mesh needs (``Null_Boundary`` — every vertex, value 666 — and
-        ``All_Boundaries``, value 1001), records the boundary / region
-        metadata attributes on the mesh, rebuilds named boundary labels for
-        wrapped DMPlex imports that only expose stacked Gmsh label sets,
-        and builds the ``Null_Boundary`` and stacked ``UW_Boundaries``
-        labels used by the boundary-condition machinery.
+        Patches the user-supplied boundary enum with ``All_Boundaries``
+        (value 1001, every exterior facet, populated by ``markBoundaryFaces``),
+        records the boundary / region metadata attributes on the mesh, rebuilds
+        named boundary labels for wrapped DMPlex imports that only expose
+        stacked Gmsh label sets, and builds the stacked ``UW_Boundaries`` label
+        used by the boundary-condition machinery.
+
+        ``Null_Boundary`` (value 666, every VERTEX of the mesh) used to be
+        added here as well. Nothing needs it: it marks no facet, so it
+        integrates nothing, and its one functional consumer — a fake natural BC
+        the solver manufactured to guarantee "a surface integral term on every
+        process" — has been removed and measured to change no answer at any
+        rank count. What it did do was mislead: any pass that reads *labelled
+        points* to find material interfaces saw every vertex of the mesh
+        labelled, which silently refused 1114 of 1114 mesh-repair candidates.
+        A boundaries enum supplied by a caller may still declare it, and old
+        checkpoints still carry the value, so the sentinel skips downstream
+        remain.
         """
         ## Patch up the boundaries to include the additional
         ## definitions that we do / might need. Note: the
@@ -900,15 +911,13 @@ def _patch_boundary_enum(
         if boundaries is None:
 
             class replacement_boundaries(Enum):
-                Null_Boundary = 666
                 All_Boundaries = 1001
 
             boundaries = replacement_boundaries
         else:
 
             @extend_enum([boundaries])
             class replacement_boundaries(Enum):
-                Null_Boundary = 666
                 All_Boundaries = 1001
 
             boundaries = replacement_boundaries
@@ -948,20 +957,6 @@ class replacement_boundaries(Enum):
                         uw.adaptivity._dm_unstack_bcs(self.dm, self.boundaries, stacked_label_name)
                         break
 
-        # Null_Boundary marks EVERY vertex with the reserved value 666 —
-        # the catch-all stratum used when a condition applies mesh-wide.
-        # Note: getStratumIS(0) on the "depth" label fetches the depth-0
-        # points, i.e. VERTICES (the old name "all_edges" was wrong).
-        depth_label = self.dm.getLabel("depth")
-        self.dm.createLabel("Null_Boundary")
-        null_boundary_label = self.dm.getLabel("Null_Boundary")
-        if depth_label and null_boundary_label:
-            vertex_stratum_is = depth_label.getStratumIS(0)
-            if vertex_stratum_is:
-                null_boundary_label.setStratumIS(
-                    boundaries.Null_Boundary.value, vertex_stratum_is
-                )
-
         ## --- UW_Boundaries label
         if self.boundaries is not None:
 

src/underworld3/meshing/segmented.py

@@ -775,7 +775,6 @@ class boundaries(Enum):
         UpperPlus = 31
         Centre = 1
         Slices = 40
-        Null_Boundary = 666
 
     meshRes = cellSize
     num_segments = numSegments

tests/parallel/test_0768_unowned_boundary_natural_bc_mpi.py

@@ -0,0 +1,190 @@
+"""Natural BCs and boundary integrals on a boundary a rank owns NO part of.
+
+This is the case the removed ``Null_Boundary`` workaround claimed to need. It
+added a fake natural BC on a label marking every vertex, because "we need a
+surface integral term somewhere on every process if we have a contribution from
+anywhere". The workaround was removed after measuring that it changed no answer
+at any rank count — value 666 marks no facet, so it integrated nothing, and
+natural BCs are registered against ``UW_Boundaries`` unconditionally on every
+rank anyway.
+
+**What these tests are.** A tripwire on current behaviour, not a proof that the
+PETSc problem the workaround was aimed at is solved. Removing a failed fix does
+not establish that the thing it failed to fix is gone. If it resurfaces, these
+tests are what will show it — and the answer then is a real fix, not a sentinel
+label that integrates nothing.
+
+**Why the suite needed them.** Before this file, the only parallel test creating
+a natural BC was ``test_0770_submesh_extract_mpi``, whose annulus boundaries are
+owned by every rank at np=2 and np=4 (``Internal`` [26,26] and [12,14,13,13]).
+No test exercised a rank owning none, so nothing would catch a rank-local skip
+being added to the DS registration — which the dead
+``bc_is = bc_label.getStratumIS(value)`` in ``_setup_discretisation`` invites.
+
+The geometry is a LONG FLAT box, so the partitioner cuts along x and the short
+end faces land on one rank each. ``test_the_unowned_case_is_exercised`` asserts
+that actually happened: without it a partitioner change could make every other
+test here pass vacuously.
+"""
+
+import numpy as np
+import pytest
+import sympy
+
+import underworld3 as uw
+
+pytestmark = [
+    pytest.mark.level_2,
+    pytest.mark.tier_a,
+    pytest.mark.mpi(min_size=2),
+    pytest.mark.timeout(300),
+]
+
+LX, LY = 4.0, 0.25
+CELL = 0.05
+
+
+def _mesh():
+    return uw.meshing.UnstructuredSimplexBox(
+        minCoords=(0.0, 0.0), maxCoords=(LX, LY), cellSize=CELL,
+        regular=False, qdegree=3)
+
+
+def _owned_facets(mesh, name):
+    """Facets of ``name`` this rank OWNS, gathered to a per-rank list on rank 0.
+
+    Owned, not merely present: a seam facet sits on every rank sharing it, and a
+    local count reports it once per sharer. A point is a ghost exactly when it is
+    a leaf of the point star-forest.
+    """
+    dm = mesh.dm
+    ghosts = set()
+    if uw.mpi.size > 1:
+        _nroots, ilocal, iremote = dm.getPointSF().getGraph()
+        ghosts = (set(range(0 if iremote is None else len(iremote)))
+                  if ilocal is None else set(int(p) for p in ilocal))
+
+    fS, fE = dm.getHeightStratum(1)
+    n = 0
+    label = dm.getLabel(name)
+    if label is not None:
+        ist = label.getStratumIS(int(mesh.boundaries[name].value))
+        # A rank owning NO point with this value gets back a valid IS of SIZE
+        # ZERO, and `.getIndices()` on it segfaults (issue #291). `is not None`
+        # does not catch it, and the rank that trips it is precisely the one
+        # this module is about.
+        if ist and ist.getSize() > 0:
+            pts = np.asarray(ist.getIndices(), dtype=np.int64)
+            pts = pts[(pts >= fS) & (pts < fE)]
+            n = sum(1 for p in pts if int(p) not in ghosts)
+    return uw.mpi.comm.allgather(n)
+
+
+def test_the_unowned_case_is_exercised():
+    """Some rank must own no facet of `Right`, or the rest of this file is vacuous."""
+    counts = _owned_facets(_mesh(), "Right")
+    assert sum(counts) > 0, f"nobody owns the boundary at all: {counts}"
+    assert min(counts) == 0, (
+        f"every rank owns part of Right ({counts}); this file is meant to test "
+        f"the case where at least one rank owns none, and no longer does")
+
+
+def test_boundary_integral_over_an_unowned_boundary():
+    """`BdIntegral` is exact on a boundary most ranks hold no part of."""
+    mesh = _mesh()
+    x, y = mesh.X
+    _T = uw.discretisation.MeshVariable("Tbi", mesh, 1, degree=2)
+
+    assert min(_owned_facets(mesh, "Left")) == 0
+
+    assert np.isclose(uw.maths.BdIntegral(mesh, 1.0, "Left").evaluate(),
+                      LY, rtol=1e-12)
+    assert np.isclose(uw.maths.BdIntegral(mesh, y, "Left").evaluate(),
+                      0.5 * LY ** 2, rtol=1e-12)
+    # A boundary every rank DOES own, as the control on the same mesh.
+    assert np.isclose(uw.maths.BdIntegral(mesh, 1.0, "Top").evaluate(),
+                      LX, rtol=1e-12)
+
+
+def test_natural_bc_on_an_unowned_boundary_is_assembled():
+    """A Neumann term on a boundary most ranks own nothing of reaches the residual.
+
+    Checked against the manufactured solution ``T = x^2``, which P2 reproduces
+    EXACTLY, with the tolerance tightened well below the bundle default of
+    ``snes_rtol`` 1e-4. The error then sits at the solver tolerance and a dropped
+    surface term is orders of magnitude, not a wobble: zeroing the Right flux
+    moves it from 3.6e-11 to 8.0e-03.
+
+    Convergence alone would prove nothing — a silently dropped surface term
+    still converges, just to the wrong field.
+    """
+    mesh = _mesh()
+    x, y = mesh.X
+    T = uw.discretisation.MeshVariable("Tnbc", mesh, 1, degree=2)
+
+    assert min(_owned_facets(mesh, "Right")) == 0
+
+    poisson = uw.systems.Poisson(mesh=mesh, u_Field=T, degree=2)
+    poisson.constitutive_model = uw.constitutive_models.DiffusionModel
+    poisson.constitutive_model.Parameters.diffusivity = 1.0
+    poisson.f = -2.0
+    poisson.add_dirichlet_bc([x ** 2], "Bottom")
+    poisson.add_dirichlet_bc([x ** 2], "Top")
+    poisson.add_natural_bc([-2 * x], "Right")
+    poisson.add_natural_bc(0.0, "Left")
+    # `petsc_options["ksp_rtol"]` is silently overridden; use `.tolerance`.
+    poisson.tolerance = 1.0e-10
+    poisson.solve()
+
+    assert poisson.snes.getConvergedReason() > 0
+
+    exact = x ** 2
+    err = uw.maths.Integral(mesh, (T.sym[0] - exact) ** 2).evaluate() ** 0.5
+    norm = uw.maths.Integral(mesh, exact ** 2).evaluate() ** 0.5
+    assert err / norm < 1.0e-8, f"relative L2 error {err / norm:.3e}"
+
+    # No condition is manufactured behind the user's back.
+    assert [bc.boundary for bc in poisson.natural_bcs] == ["Right", "Left"]
+
+
+def test_vector_natural_bc_on_an_unowned_boundary_changes_the_answer():
+    """The VECTOR path: a Stokes traction on an unowned boundary is live.
+
+    The removed workaround built its fake condition as ``(0,)*u.shape[1]``, so
+    the vector case is the one it was shaped for. Asserted by DIFFERENCE against
+    the same problem with no condition on Right, rather than against a stored
+    number: if the traction on the boundary most ranks own nothing of never
+    reached the residual, the two solves would agree.
+
+    Lid-driven, and the lid profile vanishes at both ends. A uniform body force
+    in a closed incompressible box is hydrostatic — v is identically zero, and a
+    probe whose answer is zero by construction cannot detect a dropped term
+    however many ranks own nothing of the boundary.
+    """
+    def _vrms(with_traction):
+        mesh = _mesh()
+        x, y = mesh.X
+        tag = "t" if with_traction else "n"
+        v = uw.discretisation.MeshVariable(f"v{tag}", mesh, mesh.dim, degree=2)
+        p = uw.discretisation.MeshVariable(f"p{tag}", mesh, 1, degree=1,
+                                           continuous=True)
+        stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p)
+        stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel
+        stokes.constitutive_model.Parameters.shear_viscosity_0 = 1.0
+        stokes.bodyforce = sympy.Matrix([0.0, 0.0])
+        stokes.add_dirichlet_bc((sympy.sin(sympy.pi * x / LX), 0.0), "Top")
+        stokes.add_dirichlet_bc((0.0, 0.0), "Bottom")
+        stokes.add_dirichlet_bc((0.0, 0.0), "Left")
+        if with_traction:
+            stokes.add_natural_bc((0.0, -0.5), "Right")
+        stokes.tolerance = 1.0e-10
+        stokes.solve()
+        assert stokes.snes.getConvergedReason() > 0
+        return (uw.maths.Integral(mesh, v.sym.dot(v.sym)).evaluate()
+                / (LX * LY)) ** 0.5
+
+    free = _vrms(False)
+    loaded = _vrms(True)
+    assert abs(loaded - free) / free > 1.0e-5, (
+        f"the traction on Right did not reach the residual: "
+        f"{free:.12e} vs {loaded:.12e}")

tests/test_0001_meshes.py

@@ -177,21 +177,37 @@ def test_create_solid_usbIB_3d_mesh():
 
 
 @pytest.mark.tier_a
-def test_null_boundary_marks_every_vertex():
-    """Null_Boundary (value 666) covers the depth-0 (vertex) stratum.
+def test_no_null_boundary_label_is_manufactured():
+    """A mesh does not label every vertex with the ``Null_Boundary`` sentinel.
 
-    Regression for D-37/READ-41: the label was built from a variable that
-    was only assigned inside an `if` guard (NameError-risk when the DM has
-    no "depth" label) and misnamed "all_edges" when getStratumIS(0)
-    actually fetches vertices. The rewrite must still mark every vertex.
+    ``Null_Boundary`` (666) used to be created on every mesh, marking the whole
+    depth-0 stratum. It marked no facet, so it integrated nothing, and its only
+    consumer was a fake natural BC the solver manufactured — since removed,
+    with no change to any answer at any rank count.
+
+    Asserting its ABSENCE is what protects the fix that motivated the removal:
+    a pass that looks for material interfaces by reading labelled *points* saw
+    every vertex of every UW3 mesh labelled, and silently declined all of them.
     """
     from underworld3.meshing import StructuredQuadBox
 
     mesh = StructuredQuadBox(elementRes=(4, 4))
-
-    label = mesh.dm.getLabel("Null_Boundary")
-    assert label is not None
-
-    p_start, p_end = mesh.dm.getDepthStratum(0)
-    n_vertices = p_end - p_start
-    assert label.getStratumSize(mesh.boundaries.Null_Boundary.value) == n_vertices
+    dm = mesh.dm
+
+    # Asserted against the DM's own list of label NAMES. `getLabel` on a name
+    # the DM does not have returns a non-None DMLabel wrapper with a null
+    # handle, so `is not None` reports every absent label as present; the
+    # codebase's `if label:` idiom exists for exactly this reason.
+    names = [dm.getLabelName(i) for i in range(dm.getNumLabels())]
+    assert "Null_Boundary" not in names, names
+    assert not dm.getLabel("Null_Boundary")
+    assert "Null_Boundary" not in [b.name for b in mesh.boundaries]
+
+    # The boundary that DOES mean "the whole outside" is still there, and it
+    # is a different animal: exterior FACETS, so it integrates.
+    all_bd = mesh.dm.getLabel("All_Boundaries")
+    assert all_bd is not None
+    fS, fE = mesh.dm.getHeightStratum(1)
+    pts = all_bd.getStratumIS(1001).getIndices()
+    assert len(pts) > 0
+    assert ((pts >= fS) & (pts < fE)).all(), "All_Boundaries must mark facets"

tests/test_1001_poisson_constants.py

@@ -260,9 +260,13 @@ def test_poisson_constant_in_essential_bc():
 def test_stokes_natural_bc_constant_no_recompile():
     """Stokes natural BC with constant traction — change without recompile.
 
-    Also verifies the Null_Boundary bug fix: _build() must not re-add
-    Null_Boundary on every call, which would reset is_setup and force
-    recompilation.
+    This test used to guard a bug in which ``_build()`` re-added a fake natural
+    BC on ``Null_Boundary`` on every call, resetting ``is_setup`` and forcing a
+    recompile. That fake BC has since been removed altogether — it was a
+    workaround for a PETSc problem that no longer exists, and it integrated
+    nothing (value 666 marks vertices, never a facet). So the assertion is now
+    that no such condition is manufactured at all: a solver's ``natural_bcs``
+    holds exactly what the user asked for.
     """
 
     mesh = uw.meshing.UnstructuredSimplexBox(cellSize=0.25)
@@ -291,9 +295,10 @@ def test_stokes_natural_bc_constant_no_recompile():
         f"tau_nbc not found in constants manifest: {const_names}"
     )
 
-    # Null_Boundary should appear exactly once
+    # No fake condition is manufactured: the user asked for one natural BC.
     null_count = sum(1 for bc in stokes.natural_bcs if bc.boundary == "Null_Boundary")
-    assert null_count == 1, f"Null_Boundary appears {null_count} times (expected 1)"
+    assert null_count == 0, f"Null_Boundary appears {null_count} times (expected 0)"
+    assert [bc.boundary for bc in stokes.natural_bcs] == ["Top"]
 
     # Change traction and re-solve — no recompilation
     n_before = len(_ext_dict)
@@ -307,8 +312,8 @@ def test_stokes_natural_bc_constant_no_recompile():
         f"Null_Boundary count: {sum(1 for bc in stokes.natural_bcs if bc.boundary == 'Null_Boundary')}"
     )
 
-    # Null_Boundary should still be exactly one
+    # ... and a second solve still does not manufacture one.
     null_count_2 = sum(1 for bc in stokes.natural_bcs if bc.boundary == "Null_Boundary")
-    assert null_count_2 == 1, f"Null_Boundary duplicated to {null_count_2} after second solve"
+    assert null_count_2 == 0, f"Null_Boundary appeared ({null_count_2}) on re-solve"
 
     del stokes