From 866ae55d4d2b792bc64bcdd2acdeeb120d0af564 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Tue, 4 Aug 2026 20:37:40 +1000 Subject: [PATCH 1/3] Remove the fake Null_Boundary natural BC: it fixed nothing and integrated nothing MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Every solver with a natural BC quietly gained one it never asked for, on a label marking every vertex of the mesh, justified by a comment that named no PETSc version, no error and no test: 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 Measured on the case that describes — a Neumann condition on one face of a long flat box, at rank counts where up to 7 of 8 ranks own no facet of it, checked against a manufactured solution rather than against convergence, because a silently dropped surface term still converges, just to the wrong field. The solution is x^2, which P2 reproduces exactly, and the solver tolerance is tightened to 1e-10, so the error sits at the tolerance and a missing term is unmissable rather than buried under discretisation error. The answers are bit-identical to every digit, with and without the workaround, at np = 1, 2, 3, 4, 6 and 8. A negative control that zeroes the flux term — what a dropped surface integral amounts to — moves the error from 3.6e-11 to 8.0e-03, so the probe demonstrably can fail. Boundary integrals over the same geometry are byte-for-byte identical too. Why there was nothing to fix: natural BCs are registered by PetscDSAddBoundary_UW against the consolidated UW_Boundaries label, with the boundary's own value, unconditionally on every rank. There is no rank-local skip, so the DS boundary list was already identical everywhere. (The bc_is = bc_label.getStratumIS(value) alongside it is dead — computed, never read.) And the fake BC could not have contributed anyway: value 666 marks 198 points on a test mesh, all of them VERTICES and none a facet, measured after a solve, so PETSc never completes it into something integrable. What it did do was mislead. Reading labelled points to find material interfaces saw every vertex of every UW3 mesh labelled and silently refused 1114 of 1114 mesh-repair candidates. test_stokes_natural_bc_constant_no_recompile asserted the workaround's presence — it was guarding an older bug where _build() re-added the BC on every call, resetting is_setup and forcing a recompile. That bug class cannot exist once nothing is manufactured, so the assertion becomes its opposite: a solver's natural_bcs holds exactly what the user asked for. The test's real subject, that changing a natural-BC constant does not recompile, is untouched and still passes. Probes and logs: ~/+Simulations/null_boundary_bc_hack/ Underworld development team with AI support from Claude Code --- .../cython/petsc_generic_snes_solvers.pyx | 11 ----------- tests/test_1001_poisson_constants.py | 19 ++++++++++++------- 2 files changed, 12 insertions(+), 18 deletions(-) diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index 32f46d5a6..d2b35638b 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/tests/test_1001_poisson_constants.py b/tests/test_1001_poisson_constants.py index c6beb0b3c..27dd637b0 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 From feb8640141b7c33d4242f54695ac2372ccfb76ef Mon Sep 17 00:00:00 2001 From: lmoresi Date: Wed, 5 Aug 2026 07:17:46 +1000 Subject: [PATCH 2/3] Stop manufacturing the Null_Boundary label: nothing reads it, and it misleads MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Every UW3 mesh carried a label marking EVERY VERTEX in its chart with the reserved value 666. It marks no facet, so it integrates nothing, and its only functional consumer — a fake natural BC the solver manufactured — was removed in the previous commit, measured to change no answer at any rank count. What it cost is legibility. A label's presence on a point says nothing about whether it means a material interface, and this label put a blanket over the whole vertex stratum: a mesh-repair pass that looked for interfaces by reading labelled points saw every vertex labelled and silently declined 1114 of 1114 candidates. Nothing errored; the feature was simply a no-op on every real mesh. Removed from the two injected boundaries enums, from the label-creation block, and from segmented.py's own enum. KEPT DELIBERATELY: the sentinel skips in graph.py, discretisation_mesh.py and adaptivity.py. A caller may still supply a boundaries enum declaring Null_Boundary = 666 (tests/test_1012_stokesImportedDMPlex and several docs/examples do), and a mesh reloaded from an older checkpoint still carries both the label and the enum member. Downstream readers therefore cannot assume it is absent, and _dm_unstack_bcs restores it only when the enum names it, so it becomes a no-op for our meshes and stays correct for a caller's. Backward compatibility measured rather than assumed: a checkpoint written before the removal reads at np=1 and np=3 with area exactly 1.0000000000, keeping its own label and enum member from the file. Old data keeps its labels; new meshes stop manufacturing one. Verification, against a baseline recorded before either commit: all three probe sweeps byte-for-byte identical at np = 1, 2, 3, 4, 6, 8; 93 serial tests pass, including the caller-declared-enum case; 23 parallel tests pass at np=2 and np=4, including test_0766_box_internal_boundary_mpi — the issue-#162 case, where losing these labels once made BoxInternalBoundary's natural BCs contribute nothing. The heavy solver suite matches baseline exactly, down to the same single pre-existing unrelated failure. test_null_boundary_marks_every_vertex becomes test_no_null_boundary_label_is_manufactured, which also pins the distinction from All_Boundaries: the two are not synonyms and only one of them is removable. All_Boundaries (1001) is a real geometric boundary — every exterior FACET, from markBoundaryFaces — and is itself one of the values inside UW_Boundaries, the single consolidated label the DS is pointed at. Measured on a unit box: All_Boundaries 32 edges / 0 verts, Null_Boundary 0 edges / 98 verts. The absence is asserted against the DM's list of label NAMES, not against `getLabel(...) is None`: getLabel returns a non-None DMLabel wrapper with a null handle for a name the DM does not have, so the `is None` form reports every absent label as present. Same family as issue #291. Probes and logs: ~/+Simulations/null_boundary_bc_hack/ Underworld development team with AI support from Claude Code --- .../discretisation/discretisation_mesh.py | 41 ++++++++---------- src/underworld3/meshing/segmented.py | 1 - tests/test_0001_meshes.py | 42 +++++++++++++------ 3 files changed, 47 insertions(+), 37 deletions(-) diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py index a032ca6ba..3c9f4ed84 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 662c5a15a..fdf9433bf 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/test_0001_meshes.py b/tests/test_0001_meshes.py index c88403746..3b3cd6675 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" From 2aec792325ef3890da4f1465bed846d2127a4b65 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Thu, 6 Aug 2026 18:09:04 +1000 Subject: [PATCH 3/3] Pin the case the Null_Boundary workaround claimed to need The removal of the fake Null_Boundary natural BC is pinned by tests; the claim that it was SAFE to remove was pinned by nothing. No test in the repo, before or after, exercises a natural BC on a boundary some rank owns no part of. The only parallel test that creates a natural BC at all is test_0770_submesh_extract_mpi, and there every rank owns part of both natural boundaries (Internal [26,26] at np=2, [12,14,13,13] at np=4). A long flat box (4.0 x 0.25) gets the discriminating case reliably: the partitioner cuts along x, so the short end faces land on one rank each. Measured owned facets of Right: [5,0] at np=2, [0,0,5,0] at np=4, [5,0,0,0,0,0,0,0] at np=8 -- 7 of 8 ranks owning none. Four tests: a vacuity guard that asserts some rank really does own none (without it a partitioner change makes the rest pass for no reason), a BdIntegral case, a scalar Poisson natural BC against the P2-exact manufactured solution T = x^2, and a vector Stokes traction asserted by difference rather than against a stored number. Both negative controls fire. Dropping the scalar flux term moves the relative L2 error from 3.6e-11 to 8.010e-03. Dropping the vector traction makes the two Stokes solves agree to every digit (2.582851221410e-01), which is what a term that never reached the residual looks like. This is a tripwire on current behaviour, not a proof that the PETSc problem the workaround was aimed at is gone -- removing a failed fix does not establish that. The live hazard it guards is the dead `bc_is = bc_label.getStratumIS(value)` sitting beside the DS registration loop, which invites an `if bc_is is None: continue` optimisation that would reintroduce exactly the failure the workaround claimed to prevent. Underworld development team with AI support from Claude Code --- ...st_0768_unowned_boundary_natural_bc_mpi.py | 190 ++++++++++++++++++ 1 file changed, 190 insertions(+) create mode 100644 tests/parallel/test_0768_unowned_boundary_natural_bc_mpi.py 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 000000000..3da88e679 --- /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}")