diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index 32f46d5a..d2b35638 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/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) diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py index a032ca6b..3c9f4ed8 100644 --- a/src/underworld3/discretisation/discretisation_mesh.py +++ b/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,7 +911,6 @@ def _patch_boundary_enum( if boundaries is None: class replacement_boundaries(Enum): - Null_Boundary = 666 All_Boundaries = 1001 boundaries = replacement_boundaries @@ -908,7 +918,6 @@ class replacement_boundaries(Enum): @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: diff --git a/src/underworld3/meshing/segmented.py b/src/underworld3/meshing/segmented.py index 662c5a15..fdf9433b 100644 --- a/src/underworld3/meshing/segmented.py +++ b/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 diff --git a/tests/parallel/test_0768_unowned_boundary_natural_bc_mpi.py b/tests/parallel/test_0768_unowned_boundary_natural_bc_mpi.py new file mode 100644 index 00000000..3da88e67 --- /dev/null +++ b/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}") diff --git a/tests/test_0001_meshes.py b/tests/test_0001_meshes.py index c8840374..3b3cd667 100644 --- a/tests/test_0001_meshes.py +++ b/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" diff --git a/tests/test_1001_poisson_constants.py b/tests/test_1001_poisson_constants.py index c6beb0b3..27dd637b 100644 --- a/tests/test_1001_poisson_constants.py +++ b/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