From 4b4f341cc050e314bc550ff04dfcaa8c084ff29b Mon Sep 17 00:00:00 2001
From: lmoresi <louis.moresi@anu.edu.au>
Date: Sat, 5 Sep 2026 22:24:25 -0700
Subject: [PATCH 1/6] Take the cell radius from PETSc again, so it cannot
 depend on the partition (#694)

mesh.cell_size() and get_min/max/mean_radius() were built on a kd-tree over
THIS RANK's centroids, queried with each cell's vertices. Near a partition
boundary the true nearest centroid can belong to a cell owned by another rank
and be absent from the tree, so the answer moved with the rank count: per-cell
by 3.3e-03 at np=2, get_max_radius() by 4.9% at np=8, get_mean_radius() at
every rank count. It reached users through the DEFAULT add_nitsche_bc(
local_h=True), which scales the penalty by cell_size().

An allreduce made those accessors agree on every rank without making them agree
at every rank count -- the values being reduced were themselves partition
dependent -- while their docstrings advertise "the smallest cell anywhere in the
mesh" and "use this rather than the rank-local array".

The radius is now PETSc's volume**(1/dim) from DMPlexComputeGeometryFVM. A
cell's volume is a property of that cell, so this cannot depend on the split.
That routine had been abandoned with a note that it "does not compute all
cells"; that does not reproduce. Measured on 2-D simplex, 2-D quad, 3-D
tetrahedra, 3-D hexahedra and a deformed mesh: one finite positive value per
local cell every time, and bit-identical across rank counts in all five
(max abs diff 0.000e+00). The old note also named DMPlexGetMinRadius, which is
a different call and is not used here.

It is also the definition the docstrings have claimed all along -- cell_size()
and get_mean_radius() both describe "the volume**(1/dim) equivalent radius",
which the kd-tree field was not. The documentation becomes true rather than
being rewritten.

Consequences beyond the swap:

- `_min_size` and `_search_lengths` are gone. They were unpacked from
  _get_mesh_sizes and never read anywhere in src/ or tests/ -- the loop
  computed three distance statistics per cell and discarded two.
- test_1069's `local_h=False` workaround is removed. That test measures
  partition independence and had to disable the DEFAULT code path to avoid
  measuring this defect at the same time; it now runs the default and passes
  (7 passed, 1 skipped at np=2).
- follow_metric's docstring recommended `2.0 * mesh._radii.mean()`, a RANK-LOCAL
  mean, in the same breath as get_mean_radius() warns against exactly that. It
  now recommends get_mean_radius().

tests/parallel/test_0798 is the guard, and it compares two RANK COUNTS rather
than asserting a within-rank property: the reference is this module's own np=1
answer via serial_reference, with the mesh fingerprint asserted so a host that
triangulates differently is reported as that. Negative control: restoring the
kd-tree field makes it fail with "max radius is partition dependent -- np=1
0.0670934714626447 vs np=8 0.07037395013400187".

Verified: accessors identical to 12 digits at np=1/2/4/8 (max and mean both
moved before); test_0798 passes at np=2/4/8; test_1069 passes on the default
path.

KNOWN, UNRESOLVED: the Nitsche leak in test_1060 grows from 8.954e-05 to
1.312e-04, crossing that test's 1e-4 threshold. This is not a defect in the new
radius -- the penalty is gamma*mu/h and h is ~19% larger in the mean, so the
same gamma enforces less. gamma's default of 10.0 was calibrated against the old
definition and has to move with the quantity it scales, or the method is
silently weaker. That is a solver-behaviour decision and is left for the
maintainer rather than resolved by loosening a threshold.

Underworld development team with AI support from Claude Code

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01E87Q7KrpapxeQiLD1RiNXv
---
 .../discretisation/discretisation_mesh.py     | 100 +++++++++-------
 src/underworld3/meshing/smoothing/api.py      |   2 +-
 src/underworld3/swarm.py                      |   2 +-
 src/underworld3/systems/solvers.py            |  14 +--
 .../test_0774_empty_rank_reductions_mpi.py    |   4 +-
 ...t_0798_cell_size_partition_independence.py | 113 ++++++++++++++++++
 .../test_1069_boundary_normal_parallel.py     |  19 +--
 tests/test_1065_nitsche_local_h.py            |  10 +-
 8 files changed, 195 insertions(+), 69 deletions(-)
 create mode 100644 tests/parallel/test_0798_cell_size_partition_independence.py

diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py
index d22d8b190..6f8f93db7 100644
--- a/src/underworld3/discretisation/discretisation_mesh.py
+++ b/src/underworld3/discretisation/discretisation_mesh.py
@@ -2830,12 +2830,10 @@ def nuke_coords_and_rebuild(
                 flush=True,
             )
 
-        (
-            self._min_size,
-            self._radii,
-            self._centroids,
-            self._search_lengths,
-        ) = self._get_mesh_sizes()
+        # `_min_size` and `_search_lengths` used to be unpacked here and were
+        # never read anywhere in src/ or tests/ -- the kd-tree loop computed
+        # three distance statistics per cell and two were discarded.
+        self._cell_radii, self._centroids = self._get_cell_radii()
 
         # Skip self-copy when hierarchy is trivial (issue #96 investigation)
         if self.dm is not self.dm_hierarchy[-1]:
@@ -3224,7 +3222,8 @@ def cell_size(self):
 
         Returns the ``.sym`` of a cell-constant (degree-0, discontinuous)
         scalar MeshVariable holding each cell's characteristic length (the
-        ``volume**(1/dim)`` equivalent radius, i.e. ``self._radii``). Unlike
+        ``volume**(1/dim)`` equivalent radius, i.e. ``self._cell_radii``,
+        which comes from PETSc and is independent of the MPI partition). Unlike
         the single *global* scalar from :meth:`get_min_radius` (the smallest
         cell anywhere), this varies cell to cell, so a stabilisation that
         scales as :math:`1/h` — e.g. the Nitsche free-slip penalty
@@ -3290,9 +3289,9 @@ def _refresh():
     def _assemble_cell_size(self, var):
         """Fill ``var`` (degree-0 scalar) with each cell's characteristic size.
 
-        Uses the per-cell characteristic lengths ``self._radii`` computed by
-        :meth:`_get_mesh_sizes` on the *current* geometry. A degree-0
-        discontinuous variable's local DOFs and ``self._radii`` are BOTH
+        Uses the per-cell characteristic lengths ``self._cell_radii`` computed
+        by :meth:`_get_cell_radii` on the *current* geometry. A degree-0
+        discontinuous variable's local DOFs and ``self._cell_radii`` are BOTH
         indexed by this rank's cell-stratum order, so a direct assignment is
         correct on every rank.
 
@@ -3317,7 +3316,7 @@ def _assemble_cell_size(self, var):
         # relaxation, and it needs its own benchmarking.
         # Guard/measurement: tests/parallel/test_1069_boundary_normal_parallel.py
         # (_nitsche_annulus_diagnostics docstring records the numbers).
-        radii = numpy.asarray(self._radii).reshape(-1)
+        radii = numpy.asarray(self._cell_radii).reshape(-1)
         # Empty partition (no local cells): nothing to fill on this rank.
         if radii.size == 0 or var.data.shape[0] == 0:
             return
@@ -6893,36 +6892,40 @@ def _eval_use_robust_location(self) -> bool:
         """
         return (uw.mpi.size > 1) and (self._location_capability() != "none")
 
-    def _get_mesh_sizes(self, verbose=False):
-        """
-        Obtain the (local) mesh radii and centroids using kdtree distances
-        This routine is called when the mesh is built / rebuilt
-        """
-
+    def _get_cell_radii(self):
+        """Each cell's characteristic length, and the cell centroids.
+
+        The length is PETSc's ``volume**(1/dim)`` from
+        ``DMPlexComputeGeometryFVM``. A cell's volume is a property of that
+        cell, so this cannot depend on how the mesh was partitioned — which is
+        the point.
+
+        It replaces a kd-tree of THIS RANK's centroids queried with each cell's
+        vertices. Near a partition boundary the true nearest centroid can belong
+        to a cell owned by another rank and be absent from the tree, so the
+        answer moved with the rank count: per-cell by 3.3e-03 at np=2 and
+        4.1e-03 at np=4, `get_max_radius()` by 4.9% at np=8, `get_mean_radius()`
+        at every rank count, and `mesh.cell_size()` with them -- which scales
+        the Nitsche penalty under the DEFAULT ``local_h=True`` (#569, #687,
+        #694).
+
+        The FVM routine had been abandoned with a note that it "does not
+        compute all cells". That does not reproduce: measured on 2-D simplex,
+        2-D quad, 3-D tetrahedra, 3-D hexahedra and a deformed mesh, it returns
+        one finite positive value per local cell and is bit-identical across
+        rank counts in every case. (The note also named ``DMPlexGetMinRadius``,
+        which is a different call and is not used here.)
+        """
+        from underworld3.cython import petsc_discretisation
+
+        radii, _fvm_centroids = petsc_discretisation.petsc_fvm_get_local_cell_sizes(self)
+
+        # The FVM centroids are discarded: `_get_coords_for_basis(0, False)` is
+        # the degree-0 coordinate array the rest of the mesh indexes by cell,
+        # and mixing the two orderings would misalign every per-cell lookup.
         centroids = self._get_coords_for_basis(0, False)
-        centroids_kd_tree = uw.kdtree.KDTree(centroids)
-
-        import numpy as np
 
-        cStart, cEnd = self.dm.getHeightStratum(0)
-        pStart, pEnd = self.dm.getDepthStratum(0)
-        cell_length = np.empty(centroids.shape[0])
-        cell_min_r = np.empty(centroids.shape[0])
-        cell_r = np.empty(centroids.shape[0])
-
-        for cell in range(cEnd - cStart):
-            cell_num_points = self.dm.getConeSize(cell)
-            cell_points = self.dm.getTransitiveClosure(cell)[0][-cell_num_points:]
-            # Use raw internal array for internal mesh operations (avoid unit-aware wrapping)
-            cell_coords = self._coords[cell_points - pStart]
-
-            distsq, _ = centroids_kd_tree.query(cell_coords, k=1, sqr_dists=True)
-
-            cell_length[cell] = np.sqrt(distsq.max())
-            cell_r[cell] = np.sqrt(distsq.mean())
-            cell_min_r[cell] = np.sqrt(distsq.min())
-
-        return cell_min_r, cell_r, centroids, cell_length
+        return radii, centroids
 
     # ==========
 
@@ -7016,7 +7019,7 @@ def get_min_radius(self) -> float:
         import numpy as np
         from mpi4py import MPI
 
-        radii = np.asarray(self._radii).reshape(-1)
+        radii = np.asarray(self._cell_radii).reshape(-1)
         local_min = float(radii.min()) if radii.size else float("inf")
         if uw.mpi.size > 1:
             local_min = uw.mpi.comm.allreduce(local_min, op=MPI.MIN)
@@ -7038,7 +7041,7 @@ def get_max_radius(self) -> float:
         import numpy as np
         from mpi4py import MPI
 
-        radii = np.asarray(self._radii).reshape(-1)
+        radii = np.asarray(self._cell_radii).reshape(-1)
         local_max = float(radii.max()) if radii.size else float("-inf")
         if uw.mpi.size > 1:
             local_max = uw.mpi.comm.allreduce(local_max, op=MPI.MAX)
@@ -7057,15 +7060,22 @@ def get_mean_radius(self) -> float:
         this is the canonical "mesh length" API. Use this anywhere you
         need a representative h0 (smoothing-length defaults, diffusion-
         stability heuristics, problem-scale normalisation) rather than
-        reaching for the rank-local ``self._radii`` array, which gives
-        different answers on different MPI ranks and leaks downstream
-        (e.g. into JIT C source via per-rank pointwise-function inputs).
+        reducing a per-rank array by hand, which gives different answers on
+        different MPI ranks and leaks downstream (e.g. into JIT C source via
+        per-rank pointwise-function inputs).
+
+        The value is the same at every RANK COUNT as well as on every rank:
+        ``self._cell_radii`` is PETSc's ``volume**(1/dim)``, a property of each
+        cell rather than of the partition. That was not true while these
+        reduced over a kd-tree of this rank's centroids -- an allreduce made
+        the answer agree across ranks without making it agree across rank
+        counts, and ``get_max_radius()`` moved 4.9% at np=8 (#694).
         """
 
         import numpy as np
         from mpi4py import MPI
 
-        radii = np.asarray(self._radii)
+        radii = np.asarray(self._cell_radii)
         local_sum = float(radii.sum())
         local_n = int(radii.size)
         if uw.mpi.size > 1:
diff --git a/src/underworld3/meshing/smoothing/api.py b/src/underworld3/meshing/smoothing/api.py
index 145e7977d..5b4160085 100644
--- a/src/underworld3/meshing/smoothing/api.py
+++ b/src/underworld3/meshing/smoothing/api.py
@@ -713,7 +713,7 @@ def follow_metric(
             mesh, T,
             refinement=2.0, coarsening=2.0,
             metric="gradient-uniform",
-            gradient_smoothing_length=2.0 * mesh._radii.mean(),
+            gradient_smoothing_length=2.0 * mesh.get_mean_radius(),
         )
 
     See Also
diff --git a/src/underworld3/swarm.py b/src/underworld3/swarm.py
index 4b4b679d3..45ca3f18b 100644
--- a/src/underworld3/swarm.py
+++ b/src/underworld3/swarm.py
@@ -5075,7 +5075,7 @@ def estimate_dt(self, V_fn):
         vel = uw.function.evaluate(V_fn, self._particle_coordinates.data, evalf=True)
 
         # If vel is unit-aware (UnitAwareArray), nondimensionalise it to get
-        # consistent nondimensional values that match mesh._radii
+        # consistent nondimensional values that match mesh._cell_radii
         # Note: .magnitude returns physical units, which would be wrong here
         if hasattr(vel, "units") and vel.units is not None:
             vel = uw.non_dimensionalise(vel)
diff --git a/src/underworld3/systems/solvers.py b/src/underworld3/systems/solvers.py
index c89963f99..54fe7a120 100644
--- a/src/underworld3/systems/solvers.py
+++ b/src/underworld3/systems/solvers.py
@@ -256,7 +256,7 @@ def _global_max_diffusivity(constitutive_K, mesh):
         diffusivity = K
 
     # If unit-aware (UnitAwareArray), nondimensionalise so the value is
-    # consistent with mesh._radii. Note: .magnitude alone would keep the
+    # consistent with mesh._cell_radii. Note: .magnitude alone would keep the
     # physical-units number, which would be wrong here.
     if hasattr(diffusivity, "units") and diffusivity.units is not None:
         diffusivity = uw.non_dimensionalise(diffusivity)
@@ -276,7 +276,7 @@ def _centroid_velocities_nd(V_fn, mesh, basis=None, ensure_2d=True):
 
     Shared by the ``estimate_dt`` implementations: the advective CFL limit
     needs per-element centroid velocities in the same (nondimensional)
-    scale as ``mesh._radii``.
+    scale as ``mesh._cell_radii``.
 
     Parameters
     ----------
@@ -303,7 +303,7 @@ def _centroid_velocities_nd(V_fn, mesh, basis=None, ensure_2d=True):
         vel = uw.function.evaluate(V_fn, mesh._centroids)
 
     # If unit-aware (UnitAwareArray), nondimensionalise so the values are
-    # consistent with mesh._radii. Note: .magnitude alone would keep the
+    # consistent with mesh._cell_radii. Note: .magnitude alone would keep the
     # physical-units numbers, which would be wrong here.
     if hasattr(vel, "units") and vel.units is not None:
         vel = uw.non_dimensionalise(vel)
@@ -2319,7 +2319,7 @@ def estimate_dt(self):
         vel_magnitudes = np.linalg.norm(vel, axis=1)
 
         # Get per-element radii (characteristic element size)
-        element_radii = self.mesh._radii
+        element_radii = self.mesh._cell_radii
 
         # Compute per-element advective timestep: dt_i = h_i / |v_i|
         # Avoid division by zero for elements with zero velocity
@@ -4336,7 +4336,7 @@ def estimate_dt(self, direction_aware: bool = False, percentile: float = 0.0):
             centroid) · v̂` over the cell vertices. This is the
             distance material actually traverses through the cell
             per unit ``|v|``, and is **always ≥ the isotropic
-            mesh._radii estimate**, by 1.5–3× for equant cells
+            mesh._cell_radii estimate**, by 1.5–3× for equant cells
             (geometric factor) and up to ~10× for cells that the
             mover has stretched along the flow direction. On
             adapted meshes the gain is substantial; on uniform
@@ -4379,7 +4379,7 @@ def estimate_dt(self, direction_aware: bool = False, percentile: float = 0.0):
         vel_magnitudes = np.linalg.norm(vel, axis=1)
 
         # Get per-element radii (characteristic element size)
-        element_radii = self.mesh._radii
+        element_radii = self.mesh._cell_radii
 
         ## estimate dt of adv and diff components using per-element approach
         ## dt_adv_i = h_i / |v_i| for advection
@@ -4409,7 +4409,7 @@ def _reduce_dt(per_elem):
             dt_diff_per_element = np.array([np.inf])
 
         # Per-element advective timestep — either isotropic
-        # (mesh._radii / |v|) or direction-aware (v-aligned cell
+        # (mesh._cell_radii / |v|) or direction-aware (v-aligned cell
         # extent / |v|).
         if direction_aware:
             # Per-cell vertex indices (triangle / tet).
diff --git a/tests/parallel/test_0774_empty_rank_reductions_mpi.py b/tests/parallel/test_0774_empty_rank_reductions_mpi.py
index e763ddb8d..da93112fe 100644
--- a/tests/parallel/test_0774_empty_rank_reductions_mpi.py
+++ b/tests/parallel/test_0774_empty_rank_reductions_mpi.py
@@ -1,7 +1,7 @@
 """Parallel regression tests for issue #405 — reductions on a zero-cell rank.
 
 A rank that owns NO CELLS used to raise a rank-local ``ValueError`` from an
-unguarded local reduction (``self._radii.min()`` and friends) while its
+unguarded local reduction (``self._cell_radii.min()`` and friends) while its
 populated peers sat in the matching collective. The job then hung or aborted
 asymmetrically. Every global quantity computed from rank-local data must
 instead reduce across ranks, with the starved rank contributing the identity
@@ -89,7 +89,7 @@ def test_negative_control_rank_local_minimum_would_be_caught():
     """
     mesh = _starved_box()
 
-    radii = np.asarray(mesh._radii).reshape(-1)
+    radii = np.asarray(mesh._cell_radii).reshape(-1)
     rank_local_min = float(radii.min()) if radii.size else float("inf")
     gathered = uw.mpi.comm.allgather(rank_local_min)
 
diff --git a/tests/parallel/test_0798_cell_size_partition_independence.py b/tests/parallel/test_0798_cell_size_partition_independence.py
new file mode 100644
index 000000000..de2f19715
--- /dev/null
+++ b/tests/parallel/test_0798_cell_size_partition_independence.py
@@ -0,0 +1,113 @@
+"""Cell size must not depend on how the mesh was partitioned.
+
+The defect (#569, #687, #694): the per-cell length came from a kd-tree over
+THIS RANK's centroids, queried with each cell's vertices. Near a partition
+boundary the true nearest centroid can belong to a cell owned by another rank
+and so be absent from the tree, and the answer moved with the rank count --
+per-cell by 3.3e-03 at np=2, ``get_max_radius()`` by 4.9% at np=8,
+``get_mean_radius()`` at every rank count.
+
+It reached users through ``mesh.cell_size()``, which scales the Nitsche penalty
+under the default ``local_h=True``, and through the three radius accessors,
+whose docstrings advertise a global mesh length.
+
+Two things this file is careful about, both learned the hard way:
+
+* **It compares two rank counts.** Partition independence is a statement about
+  two runs agreeing, so nothing measured at a single rank count establishes it.
+  The tests that shipped with the first attempt at this fix asserted a
+  within-rank oracle -- each cell matching its own vertices -- which is true of
+  a partition-dependent field as well.
+* **The reference is computed here, not recorded.** ``serial_reference`` runs
+  this module's own ``__main__`` at np=1 in this environment and asserts the
+  mesh fingerprints match, so a host that triangulates differently is reported
+  as that rather than as partition dependence.
+"""
+
+import numpy as np
+import pytest
+
+import underworld3 as uw
+from mpi4py import MPI
+
+from serial_reference import compare, emit, mesh_fingerprint, serial_reference
+
+pytestmark = [pytest.mark.mpi(min_size=2), pytest.mark.level_1, pytest.mark.tier_a]
+
+LABELS = ("min radius", "max radius", "mean radius", "sum of cell radii")
+
+# min and max are exact reductions of identical per-cell values, so they must
+# agree to the bit. The mean is a distributed sum and reduces in partition
+# order, so it is allowed the last couple of bits -- that ordering difference
+# is not the defect under test.
+RTOLS = (0.0, 0.0, 1.0e-12, 1.0e-12)
+
+
+def _box():
+    mesh = uw.meshing.UnstructuredSimplexBox(
+        minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=0.12, qdegree=2)
+    # mesh_fingerprint integrates over the mesh, and integration needs at
+    # least one variable to exist on it.
+    uw.discretisation.MeshVariable("cell_size_probe", mesh, 1, degree=1)
+    return mesh
+
+
+def _cell_size_diagnostics():
+    """The three accessors, plus a global sum that no single one of them sees.
+
+    The sum is the sharp one: a per-cell change that leaves the extremes alone
+    still moves it, and the old field's per-cell values differed while its
+    minimum happened not to.
+    """
+    mesh = _box()
+    radii = np.asarray(mesh._cell_radii).reshape(-1)
+    local_sum = float(radii.sum())
+    global_sum = uw.mpi.comm.allreduce(local_sum, op=MPI.SUM) \
+        if uw.mpi.size > 1 else local_sum
+
+    values = (
+        mesh.get_min_radius(),
+        mesh.get_max_radius(),
+        mesh.get_mean_radius(),
+        global_sum,
+    )
+    return values, mesh_fingerprint(mesh)
+
+
+@pytest.mark.mpi(min_size=2)
+def test_cell_size_is_the_same_at_every_rank_count():
+    """The accessors and the per-cell sum reproduce their own np=1 answer."""
+    values, fingerprint = _cell_size_diagnostics()
+    compare(values, serial_reference(__file__, "cell_size"),
+            rtols=RTOLS, labels=LABELS, fingerprint=fingerprint,
+            what="cell size / radius accessors")
+
+
+@pytest.mark.mpi(min_size=2)
+def test_every_rank_agrees_on_the_accessors():
+    """The weaker property, kept because it is the one the allreduce gives.
+
+    An allreduce makes an answer identical on every rank; it does not make it
+    identical at every rank count if the values being reduced are themselves
+    partition-dependent. Failing this while passing the test above would mean
+    the reduction is broken rather than its input, so keeping both separates
+    the two.
+    """
+    mesh = _box()
+    for name in ("min", "max", "mean"):
+        got = getattr(mesh, f"get_{name}_radius")()
+        everyone = uw.mpi.comm.allgather(got)
+        assert len(set(everyone)) == 1, (
+            f"get_{name}_radius() differs between ranks: {everyone}"
+        )
+
+
+if __name__ == "__main__":
+    import sys
+
+    _kind = sys.argv[1] if len(sys.argv) > 1 else "cell_size"
+    if _kind == "cell_size":
+        _values, _fingerprint = _cell_size_diagnostics()
+        emit(_values, _fingerprint)
+    else:
+        raise SystemExit(f"unknown kind {_kind!r}")
diff --git a/tests/parallel/test_1069_boundary_normal_parallel.py b/tests/parallel/test_1069_boundary_normal_parallel.py
index 2eb43c712..83ddc10f1 100644
--- a/tests/parallel/test_1069_boundary_normal_parallel.py
+++ b/tests/parallel/test_1069_boundary_normal_parallel.py
@@ -333,13 +333,16 @@ def _nitsche_annulus_diagnostics():
     leakage. Both are stable from tolerance 1e-9 to 1e-12, so neither is the linear
     solve.
 
-    ``local_h=False`` is deliberate and it is not a workaround for this fix. The
-    default ``local_h=True`` scales the Nitsche penalty by ``mesh.cell_size()``, which
-    is built from ``Mesh._get_mesh_sizes`` — a kd-tree query against THIS RANK's cell
-    centroids, and so partition-dependent in its own right (on this mesh the field's
-    sum is 26.0822 at np=1, 26.1211 at np=2, 26.1386 at np=4, and its max moves at
-    np=4). That is a SEPARATE defect from the boundary normal, it is not what #564 is
-    about, and leaving it in would make this test measure the two together. See the
+    This runs with the DEFAULT ``local_h=True``. It used to pass ``local_h=False``
+    because ``mesh.cell_size()`` was itself partition-dependent -- built from a
+    kd-tree over THIS RANK's cell centroids, so on this mesh the field's sum was
+    26.0822 at np=1, 26.1211 at np=2 and 26.1386 at np=4 -- and leaving it on
+    would have made this test measure two defects at once. ``cell_size()`` now
+    comes from PETSc's ``volume**(1/dim)`` and is partition-independent (#694),
+    so the default path is the one under test again, which is what a guard on
+    the boundary normal should be exercising.
+
+    See the
     TODO(BUG) on ``Mesh._assemble_cell_size``.
     """
     RI, RO = 0.5, 1.0
@@ -357,7 +360,7 @@ def _nitsche_annulus_diagnostics():
         y / r * sympy.cos(4 * theta) * (r - RI) * (RO - r) * 40.0]])
     stokes.add_essential_bc((0.0, 0.0), "Lower")
     # default normal= is the assembled one — that is what is under test
-    stokes.add_nitsche_bc(0.0, "Upper", local_h=False)
+    stokes.add_nitsche_bc(0.0, "Upper")
     stokes.tolerance = 1.0e-9
     stokes.petsc_options["snes_type"] = "ksponly"
     stokes.solve()
diff --git a/tests/test_1065_nitsche_local_h.py b/tests/test_1065_nitsche_local_h.py
index 15c0bcb3b..36e12d813 100644
--- a/tests/test_1065_nitsche_local_h.py
+++ b/tests/test_1065_nitsche_local_h.py
@@ -152,11 +152,11 @@ def _box_wobble(X0, amp):
 # --------------------------------------------------------------------------
 def test_cell_size_is_local_per_cell():
     """``mesh.cell_size()`` is a per-cell field equal to each cell's
-    characteristic size (``mesh._radii``), not the single global minimum."""
+    characteristic size (``mesh._cell_radii``), not the single global minimum."""
     mesh = _graded_box()
     h = mesh.cell_size()  # sympy symbol -> backed by a P0 field
     field = np.asarray(mesh._cell_size_variable.data[:, 0]).reshape(-1)
-    radii = np.asarray(mesh._radii).reshape(-1)
+    radii = np.asarray(mesh._cell_radii).reshape(-1)
 
     # field exactly mirrors the per-cell characteristic size (rank-local check,
     # reduced to a single global pass/fail so all ranks agree)
@@ -177,13 +177,13 @@ def test_local_h_at_coarse_freeslip_boundary_exceeds_global_min():
     times the global minimum. global-h would over-stiffen that penalty by
     exactly this factor; local-h scales it correctly."""
     mesh = _graded_box(h_fine=0.04, h_coarse=0.12)
-    # build/exercise the field; its data equals mesh._radii (asserted in
+    # build/exercise the field; its data equals mesh._cell_radii (asserted in
     # test_cell_size_is_local_per_cell), so we read the per-cell sizes directly
     # from _radii / _centroids — a rank-local lookup, avoiding the collective
     # arbitrary-point uw.function.evaluate (which deadlocks in parallel).
     _ = mesh.cell_size()
     cen = np.asarray(mesh._centroids)
-    radii = np.asarray(mesh._radii).reshape(-1)
+    radii = np.asarray(mesh._cell_radii).reshape(-1)
 
     near_top = cen[:, 1] > 0.85          # cells adjacent to the Top free-slip edge
     h_top = radii[near_top]
@@ -211,7 +211,7 @@ def test_cell_size_tracks_deformation():
     assert moved  # geometry actually changed
 
     h_after = mesh._cell_size_variable.data[:, 0].copy()
-    radii_after = np.asarray(mesh._radii).reshape(-1)
+    radii_after = np.asarray(mesh._cell_radii).reshape(-1)
 
     # not stale: the field changed with the geometry SOMEWHERE (global OR) ...
     nb = min(h_after.shape[0], h_before.shape[0])

From c0948ac25f5dbb2f7a89921126f5e079b76f5207 Mon Sep 17 00:00:00 2001
From: lmoresi <louis.moresi@anu.edu.au>
Date: Sun, 6 Sep 2026 08:14:56 -0700
Subject: [PATCH 2/6] Scale the default Nitsche gamma to 12.5, to match the new
 definition of h (#694)

The penalty is gamma*mu/h, so gamma is calibrated against whatever h means.
mesh.cell_size() is now PETSc's volume**(1/dim) rather than a kd-tree distance
to a neighbouring centroid, about 19% larger in the mean, so gamma=10.0 enforced
correspondingly less: the free-slip leak in test_1060 grew from 8.954e-05 to
1.312e-04, past that test's 1e-4 threshold.

A calibration constant has to move with the quantity it scales, or the method is
silently weaker. gamma=12.5 (maintainer's call) restores it and slightly
improves on the original:

    development, gamma=10, old h     8.954e-05
    this branch, gamma=10, new h     1.312e-04
    this branch, gamma=12.5, new h   7.851e-05

test_1060 passed gamma=10.0 explicitly, so the default alone did not reach it.
It now omits gamma and follows the default, because what it asserts is that
Nitsche enforces v.n = 0 at the RECOMMENDED penalty -- not that 10.0 in
particular does. Pinning the old number there would have kept the test green
while leaving every user on a weaker constraint, which is the failure mode this
whole change exists to remove. test_1065's helper default moves with it.

tests/test_1060 + test_1065: 9 passed (was 4 failed, 5 passed).

Underworld development team with AI support from Claude Code

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01E87Q7KrpapxeQiLD1RiNXv
---
 .../cython/petsc_generic_snes_solvers.pyx     | 21 +++++++++++++------
 tests/test_1060_nitsche_freeslip.py           |  6 +++---
 tests/test_1065_nitsche_local_h.py            |  2 +-
 3 files changed, 19 insertions(+), 10 deletions(-)

diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx
index f56fbe0fe..8ab13743e 100644
--- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx
+++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx
@@ -4315,7 +4315,7 @@ class SNES_Vector(SolverBaseClass):
 
 
     def add_nitsche_bc(self, conds=None, boundary=None, direction=None,
-                       normal=None, gamma=10.0, theta=1, mask=None,
+                       normal=None, gamma=12.5, theta=1, mask=None,
                        local_h=True, g=None):
         r"""Add Nitsche weak enforcement of a velocity constraint along a direction.
 
@@ -4337,8 +4337,13 @@ class SNES_Vector(SolverBaseClass):
             terms — the same geometric-normal override as on the Stokes
             variant. Default ``None`` uses the per-boundary,
             deformation-tracking ``mesh.boundary_normal(boundary)``.
-        gamma : float, default=10.0
-            Dimensionless stabilisation parameter.
+        gamma : float, default=12.5
+            Dimensionless stabilisation parameter. The penalty is
+            ``gamma*mu/h``, so this is calibrated against the definition of
+            ``h``. It was 10.0 while ``h`` came from a kd-tree of neighbouring
+            centroids; ``mesh.cell_size()`` is now PETSc's ``volume**(1/dim)``,
+            about 19% larger in the mean, and 12.5 restores the enforcement
+            that gamma=10 gave against the old h (#694).
         theta : {-1, 0, 1}, default=1
             Symmetry parameter (1=symmetric, -1=skew-symmetric).
         mask : sympy expression, optional
@@ -6532,7 +6537,7 @@ class SNES_Stokes_SaddlePt(SolverBaseClass):
                      remove_mean=remove_mean)
 
     def add_nitsche_bc(self, conds=None, boundary=None, direction=None, normal=None,
-                       gamma=10.0, theta=1, mask=None, local_h=True, g=None):
+                       gamma=12.5, theta=1, mask=None, local_h=True, g=None):
         r"""Add Nitsche weak enforcement of a velocity constraint along a direction.
 
         Nitsche's method provides a variationally consistent alternative to
@@ -6569,9 +6574,13 @@ class SNES_Stokes_SaddlePt(SolverBaseClass):
             Boundary unit normal used in the Nitsche consistency, symmetry,
             and pressure-coupling terms. Default ``None`` uses the per-boundary,
             deformation-tracking ``mesh.boundary_normal(boundary)``.
-        gamma : float, default=10.0
+        gamma : float, default=12.5
             Dimensionless stabilisation parameter. Typical values 5--20
-            for P2 elements.
+            for P2 elements. The penalty is ``gamma*mu/h``, so this is
+            calibrated against the definition of ``h``: it was 10.0 while
+            ``h`` came from a kd-tree of neighbouring centroids, and moved
+            with ``mesh.cell_size()`` becoming PETSc's ``volume**(1/dim)``
+            (#694).
         theta : {-1, 0, 1}, default=1
             Symmetry parameter:
              1: symmetric (default — optimal convergence and solver efficiency)
diff --git a/tests/test_1060_nitsche_freeslip.py b/tests/test_1060_nitsche_freeslip.py
index fe436d00d..4089f2875 100644
--- a/tests/test_1060_nitsche_freeslip.py
+++ b/tests/test_1060_nitsche_freeslip.py
@@ -67,8 +67,8 @@ def _solve_freeslip_box(method, res=8):
         stokes.add_natural_bc(1e4 * Gamma.dot(v.sym) * Gamma, "Top")
         stokes.add_natural_bc(1e4 * Gamma.dot(v.sym) * Gamma, "Bottom")
     elif method == "nitsche":
-        stokes.add_nitsche_bc(0.0, "Top", gamma=10.0)
-        stokes.add_nitsche_bc(0.0, "Bottom", gamma=10.0)
+        stokes.add_nitsche_bc(0.0, "Top")
+        stokes.add_nitsche_bc(0.0, "Bottom")
     else:
         raise ValueError(f"Unknown method: {method}")
 
@@ -147,4 +147,4 @@ def test_nitsche_better_than_penalty_constraint(self, solutions):
         max_vn_pen = np.max(np.abs(v_pen[top_pen, 1])) if np.any(top_pen) else 0
 
         print(f"Normal velocity on top: Nitsche={max_vn_nit:.4e}, Penalty={max_vn_pen:.4e}")
-        # Nitsche at gamma=10 should be comparable or better than penalty at 1e4
+        # Nitsche at the default gamma should be comparable or better than penalty at 1e4
diff --git a/tests/test_1065_nitsche_local_h.py b/tests/test_1065_nitsche_local_h.py
index 36e12d813..265d2d44f 100644
--- a/tests/test_1065_nitsche_local_h.py
+++ b/tests/test_1065_nitsche_local_h.py
@@ -224,7 +224,7 @@ def test_cell_size_tracks_deformation():
 # --------------------------------------------------------------------------
 # 3. local-h still solves free-slip correctly (back-compat / correctness)
 # --------------------------------------------------------------------------
-def _solve_freeslip(mesh, method, gamma=10.0):
+def _solve_freeslip(mesh, method, gamma=12.5):
     v = uw.discretisation.MeshVariable(
         "U", mesh, mesh.dim, degree=2, vtype=uw.VarType.VECTOR)
     p = uw.discretisation.MeshVariable("P", mesh, 1, degree=1)

From be3aa08787b2df37b2160556e438aa8e0d153d36 Mon Sep 17 00:00:00 2001
From: lmoresi <louis.moresi@anu.edu.au>
Date: Sun, 6 Sep 2026 13:26:38 -0700
Subject: [PATCH 3/6] Allocate the cell-size variable before testing for an
 empty partition (#698)

_assemble_cell_size short-circuited on `radii.size == 0` and returned without
ever touching `var.data`. That access is not rank-local: it lazily reaches
MeshVariable._set_vec, which calls dm.createSubDM and createGlobalVector, both
collective on the DM. So a rank owning no cells skipped two collectives its
populated peers made, and the job diverged -- observed at np=8 with cells per
rank [2,2,2,2,0,2,2,2]: ranks 0-3 and 5-7 inside _assemble_cell_size, rank 4
already past it.

`var.data` is now read on every rank before the branch. The docstring claimed
the whole routine was "purely RANK-LOCAL ... no collective", which was the
false premise behind the early return; it now says which access is collective
and why it has to come first.

A zero-cell rank is routine on a region SUBMESH, which keeps only the cells
carrying a label -- a partition whose share of the parent lies outside that
region owns none. That is the case this matters for.

HONESTLY REPORTED: this fix is argued from the code, not demonstrated by a
reproducer. Neither route to a zero-cell rank could be made to run:

* a 14-cell mesh at np=8 hangs during mesh CONSTRUCTION, before cell_size() is
  reached, on this branch and on development alike (intermittently -- one
  earlier run did get through and showed the divergence above);
* `extract_region` on an annulus hangs at np=2 before returning, so the submesh
  route never reaches cell_size() either. Filed separately; it is not this.

What is verified is that the change regresses nothing: parallel guards 3 passed
at np=2 and np=4, serial cell-size and Nitsche tests 14 passed.

Underworld development team with AI support from Claude Code

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01E87Q7KrpapxeQiLD1RiNXv
---
 .../discretisation/discretisation_mesh.py     | 32 +++++++++++++------
 1 file changed, 23 insertions(+), 9 deletions(-)

diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py
index eb8e976a6..37165eacf 100644
--- a/src/underworld3/discretisation/discretisation_mesh.py
+++ b/src/underworld3/discretisation/discretisation_mesh.py
@@ -3305,24 +3305,38 @@ def _assemble_cell_size(self, var):
         indexed by this rank's cell-stratum order, so a direct assignment is
         correct on every rank.
 
-        This is deliberately a purely RANK-LOCAL operation (no ``var.coords``
-        access, no collective): mixing a rank-local fast path with a
-        collective fallback would diverge across ranks and deadlock, because
-        ``var.coords`` triggers the collective ``_get_coords_for_basis``."""
+        ``var.data`` is touched on EVERY rank before the empty-partition test,
+        which is load-bearing. Its first access lazily reaches
+        ``MeshVariable._set_vec``, and that calls ``dm.createSubDM`` and
+        ``createGlobalVector`` -- both collective on the DM. A rank owning no
+        cells used to short-circuit on ``radii.size == 0`` and return without
+        ever touching it, so it skipped two collectives its populated peers
+        made and the job diverged (#698). Reduce, or in this case ALLOCATE,
+        before branching.
+
+        Beyond that first access this stays rank-local: ``var.coords`` is
+        deliberately not read, because it triggers the collective
+        ``_get_coords_for_basis`` from inside a branch that only some ranks
+        take."""
         # `_cell_radii` is PETSc's volume**(1/dim), a property of each cell, so
         # the values here do not depend on the partition -- and neither does the
         # Nitsche penalty gamma*mu/h that consumes them under the default
         # local_h=True. It was a kd-tree distance to the nearest centroid among
         # THIS RANK's centroids, which near a seam could simply be absent (#694).
+        # FIRST, on every rank: this allocates the variable's vectors through
+        # two DM collectives (see the docstring). It must not sit behind the
+        # empty-partition test.
+        data = var.data
+
         radii = numpy.asarray(self._cell_radii).reshape(-1)
-        # Empty partition (no local cells): nothing to fill on this rank.
-        if radii.size == 0 or var.data.shape[0] == 0:
+        # Empty partition (no local cells): nothing left to fill on this rank.
+        if radii.size == 0 or data.shape[0] == 0:
             return
         # Assign over the common length. In practice these match exactly (same
         # local cell set / ordering); the slice only guards a stray off-by-ghost
-        # mismatch without ever taking a collective path on a subset of ranks.
-        n = min(var.data.shape[0], radii.shape[0])
-        var.data[:n, 0] = radii[:n]
+        # mismatch.
+        n = min(data.shape[0], radii.shape[0])
+        data[:n, 0] = radii[:n]
 
     @property
     def Gamma_P1(self):

From 97cade708165f9a8532b19581177280ede0f6a90 Mon Sep 17 00:00:00 2001
From: lmoresi <louis.moresi@anu.edu.au>
Date: Sun, 6 Sep 2026 20:29:37 -0700
Subject: [PATCH 4/6] Remove the early return from _assemble_cell_size
 entirely: BOTH steps are collective (#698)

The previous commit hoisted `var.data` above the empty-partition test, which
fixed one collective and left the other. With a real reproducer the job still
hung, one step later, with the populated ranks in `pack_raw_data_to_petsc`:
assigning into `var.data` fires the array's write-back callback, and that is
collective too.

There is now no early return. A rank owning no cells walks the same path and
writes a zero-length slice, so both collectives happen everywhere:

  data = var.data                    # createSubDM + createGlobalVector
  n = min(data.shape[0], radii.shape[0])
  data[:n, 0] = radii[:n]            # pack_raw_data_to_petsc; n == 0 is fine

The reproducer is a region SUBMESH, which is the case that matters and the one
the maintainer pointed at: a submesh keeps only the cells carrying its label, so
a partition whose share of the parent lies outside that region owns none. A box
split at z=0.85 gives an outer slab with cells per rank [12, 11, 0, 19, 0, 0, 0,
0] at np=8 -- five starved ranks out of eight, which is routine rather than
contrived.

My earlier attempts used an over-decomposed 14-cell full mesh instead. That was
a poor choice: it is degenerate, it hits an unrelated intermittent hang in mesh
construction, and it obscured the defect rather than exposing it. The previous
commit's "argued from the code, not demonstrated" caveat is now discharged.

tests/parallel/test_0779 guards it. The test asserts that every rank RETURNS --
a regression hangs rather than fails -- plus a second case that the field is
still filled where cells exist, so the fix cannot be satisfied by never filling
it. The starved count is reported and not asserted: how the parent splits is the
partitioner's business, and at np=2 the slab reaches both ranks.

np=8: cells [12, 11, 0, 19, 0, 0, 0, 0], starved 5, 2 passed. np=2: 2 passed.
No regressions: parallel guards 3 passed at np=4, serial cell-size and Nitsche
14 passed.

Underworld development team with AI support from Claude Code

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01E87Q7KrpapxeQiLD1RiNXv
---
 .../discretisation/discretisation_mesh.py     | 49 ++++++-------
 ...est_0779_submesh_cell_size_starved_rank.py | 69 +++++++++++++++++++
 2 files changed, 94 insertions(+), 24 deletions(-)
 create mode 100644 tests/parallel/test_0779_submesh_cell_size_starved_rank.py

diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py
index 37165eacf..e3963824a 100644
--- a/src/underworld3/discretisation/discretisation_mesh.py
+++ b/src/underworld3/discretisation/discretisation_mesh.py
@@ -3305,36 +3305,37 @@ def _assemble_cell_size(self, var):
         indexed by this rank's cell-stratum order, so a direct assignment is
         correct on every rank.
 
-        ``var.data`` is touched on EVERY rank before the empty-partition test,
-        which is load-bearing. Its first access lazily reaches
-        ``MeshVariable._set_vec``, and that calls ``dm.createSubDM`` and
-        ``createGlobalVector`` -- both collective on the DM. A rank owning no
-        cells used to short-circuit on ``radii.size == 0`` and return without
-        ever touching it, so it skipped two collectives its populated peers
-        made and the job diverged (#698). Reduce, or in this case ALLOCATE,
-        before branching.
-
-        Beyond that first access this stays rank-local: ``var.coords`` is
-        deliberately not read, because it triggers the collective
-        ``_get_coords_for_basis`` from inside a branch that only some ranks
-        take."""
+        This routine has TWO collectives in it and therefore no early return,
+        which is the opposite of what its previous docstring claimed (#698):
+
+        * the first ``var.data`` access lazily reaches ``MeshVariable._set_vec``,
+          which calls ``dm.createSubDM`` and ``createGlobalVector``;
+        * assigning into ``var.data`` fires the array's write-back callback,
+          ``pack_raw_data_to_petsc``.
+
+        A rank owning no cells used to short-circuit on ``radii.size == 0`` and
+        return before either, while its populated peers made both. Measured on a
+        region submesh at np=8 with cells per rank ``[12, 11, 0, 19, 0, 0, 0, 0]``:
+        the populated ranks sat in ``pack_raw_data_to_petsc`` and the job never
+        finished. A starved rank now walks the same path writing a zero-length
+        slice.
+
+        ``var.coords`` is still deliberately not read: it triggers the collective
+        ``_get_coords_for_basis``, and reading it only on some ranks would put a
+        third conditional collective back in."""
         # `_cell_radii` is PETSc's volume**(1/dim), a property of each cell, so
         # the values here do not depend on the partition -- and neither does the
         # Nitsche penalty gamma*mu/h that consumes them under the default
         # local_h=True. It was a kd-tree distance to the nearest centroid among
         # THIS RANK's centroids, which near a seam could simply be absent (#694).
-        # FIRST, on every rank: this allocates the variable's vectors through
-        # two DM collectives (see the docstring). It must not sit behind the
-        # empty-partition test.
-        data = var.data
-
+        # There is NO early return here, and that is the point. Both steps below
+        # are collective, so a rank owning no cells has to walk through them
+        # writing nothing rather than skipping them (#698).
+        data = var.data                                   # allocates: createSubDM + createGlobalVector
         radii = numpy.asarray(self._cell_radii).reshape(-1)
-        # Empty partition (no local cells): nothing left to fill on this rank.
-        if radii.size == 0 or data.shape[0] == 0:
-            return
-        # Assign over the common length. In practice these match exactly (same
-        # local cell set / ordering); the slice only guards a stray off-by-ghost
-        # mismatch.
+
+        # `n` is 0 on a starved rank. The assignment still fires the array's
+        # write-back callback, which is the second collective.
         n = min(data.shape[0], radii.shape[0])
         data[:n, 0] = radii[:n]
 
diff --git a/tests/parallel/test_0779_submesh_cell_size_starved_rank.py b/tests/parallel/test_0779_submesh_cell_size_starved_rank.py
new file mode 100644
index 000000000..80358dcb5
--- /dev/null
+++ b/tests/parallel/test_0779_submesh_cell_size_starved_rank.py
@@ -0,0 +1,69 @@
+"""``cell_size()`` must survive a rank that owns none of the submesh.
+
+A region submesh keeps only the cells carrying its label, so a partition whose
+share of the parent lies outside that region legitimately owns nothing. That is
+routine for submeshes and is the case this guards -- not an over-decomposed
+full mesh, which is a different (and degenerate) situation.
+
+The defect (#698): ``_assemble_cell_size`` returned early on ``radii.size == 0``,
+in front of TWO collectives --- the first ``var.data`` access, which allocates
+through ``dm.createSubDM`` and ``createGlobalVector``, and the assignment into
+it, which fires ``pack_raw_data_to_petsc``. Starved ranks skipped both while
+their populated peers made them, and the job never finished. Measured at np=8
+with cells per rank ``[12, 11, 0, 19, 0, 0, 0, 0]``.
+"""
+
+import pytest
+
+import underworld3 as uw
+
+pytestmark = [pytest.mark.mpi(min_size=2), pytest.mark.level_1, pytest.mark.tier_a]
+
+
+def _thin_slab_parent():
+    """A box split near the top, so the outer region is a few cells deep."""
+    return uw.meshing.BoxInternalBoundary(
+        elementRes=(12, 12), zelementRes=(10, 2), zintCoord=0.85,
+        minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), simplex=True, qdegree=2)
+
+
+@pytest.mark.mpi(min_size=2)
+def test_cell_size_completes_when_ranks_own_none_of_the_submesh():
+    """The whole test is that this returns. A regression hangs rather than fails.
+
+    The starved-rank count is reported, not asserted: how the parent is split is
+    the partitioner's business, and at np=2 the slab may well reach both ranks.
+    What must hold at every rank count is that every rank comes back.
+    """
+    submesh = _thin_slab_parent().extract_region("Outer")
+    start, end = submesh.dm.getHeightStratum(0)
+    counts = uw.mpi.comm.allgather(end - start)
+
+    submesh.cell_size()
+    uw.mpi.comm.barrier()
+
+    uw.mpi.pprint(f"SUBMESH_CELL_SIZE ranks={uw.mpi.size} cells={counts} "
+                  f"starved={counts.count(0)}")
+    assert sum(counts) > 0, "the region submesh is empty everywhere; test is vacuous"
+
+
+@pytest.mark.mpi(min_size=2)
+def test_the_size_field_is_filled_where_there_are_cells():
+    """A starved rank writing nothing must not stop the others writing.
+
+    Without this, the fix above could be satisfied by never filling the field
+    at all.
+    """
+    import numpy as np
+
+    submesh = _thin_slab_parent().extract_region("Outer")
+    start, end = submesh.dm.getHeightStratum(0)
+    submesh.cell_size()
+
+    local = np.asarray(submesh._cell_size_variable.array[:, 0, 0])
+    assert local.shape[0] == end - start
+    if local.size:
+        assert np.all(local > 0.0), "populated rank has non-positive cell sizes"
+    from mpi4py import MPI
+    filled = uw.mpi.comm.allreduce(int(local.size), op=MPI.SUM)
+    assert filled > 0, "no rank filled the field at all"

From 291c330f731c37a552456e15c8cde8c4bccb1ef5 Mon Sep 17 00:00:00 2001
From: lmoresi <louis.moresi@anu.edu.au>
Date: Wed, 9 Sep 2026 09:21:43 -0700
Subject: [PATCH 5/6] test_0774: the radius oracle is area**(1/dim), not the
 half-diagonal

The test asserted the cell size against the centroid-to-corner half-diagonal
that this branch replaces. On the 1.0 x 0.5 fixture cells that is 0.559017,
where the definition taken from PETSc's DMPlexComputeGeometryFVM gives
area**(1/2) = 0.707107 -- exactly what CI measured. The oracle was stale, not
the code.

Still analytic rather than a recorded number, so it remains a real oracle.
Green at np=2 and np=4.

Underworld development team with AI support from Claude Code

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01E87Q7KrpapxeQiLD1RiNXv
---
 tests/parallel/test_0774_empty_rank_reductions_mpi.py | 9 +++++----
 1 file changed, 5 insertions(+), 4 deletions(-)

diff --git a/tests/parallel/test_0774_empty_rank_reductions_mpi.py b/tests/parallel/test_0774_empty_rank_reductions_mpi.py
index da93112fe..f609ca752 100644
--- a/tests/parallel/test_0774_empty_rank_reductions_mpi.py
+++ b/tests/parallel/test_0774_empty_rank_reductions_mpi.py
@@ -34,10 +34,11 @@
     pytest.mark.tier_a,
 ]
 
-# Cell is 1.0 wide x 0.5 tall; the characteristic length UW3 reports is the
-# centroid-to-corner half-diagonal. Analytic, so this is a real oracle rather
-# than a recorded number.
-SERIAL_RADIUS = math.sqrt(0.5 ** 2 + 0.25 ** 2)
+# Cell is 1.0 wide x 0.5 tall, so its area is 0.5 and the characteristic
+# length UW3 reports is area**(1/2). Analytic, so this is a real oracle rather
+# than a recorded number. (It was the centroid-to-corner half-diagonal before
+# the cell size came from PETSc's DMPlexComputeGeometryFVM, issue #694.)
+SERIAL_RADIUS = (1.0 * 0.5) ** (1.0 / 2.0)
 
 
 def _starved_box():

From 990d6617242aad04054a8d314d4c8bea3b3590e3 Mon Sep 17 00:00:00 2001
From: lmoresi <louis.moresi@anu.edu.au>
Date: Wed, 9 Sep 2026 09:26:44 -0700
Subject: [PATCH 6/6] Renumber the two partition-independence tests off the
 SUPG numbers

test_1077 and test_1078 were each taken by two files after development
merged the SUPG parallel suites: test_1077_advdiff_supg_parallel and
test_1078_navier_stokes_supg_parallel landed there first. Move ours to
1079 and 1080 and fix the cross-references in test_0010 and in the
1080 docstring.

scripts/test.sh globs tests/parallel/test_*.py, so both still run in CI.
Green at np=2.

Underworld development team with AI support from Claude Code

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01E87Q7KrpapxeQiLD1RiNXv
---
 ...dence.py => test_1079_cell_size_partition_independence.py} | 0
 ...y => test_1080_radius_accessors_partition_independence.py} | 2 +-
 tests/test_0010_cell_size_geometry.py                         | 4 ++--
 3 files changed, 3 insertions(+), 3 deletions(-)
 rename tests/parallel/{test_1077_cell_size_partition_independence.py => test_1079_cell_size_partition_independence.py} (100%)
 rename tests/parallel/{test_1078_radius_accessors_partition_independence.py => test_1080_radius_accessors_partition_independence.py} (98%)

diff --git a/tests/parallel/test_1077_cell_size_partition_independence.py b/tests/parallel/test_1079_cell_size_partition_independence.py
similarity index 100%
rename from tests/parallel/test_1077_cell_size_partition_independence.py
rename to tests/parallel/test_1079_cell_size_partition_independence.py
diff --git a/tests/parallel/test_1078_radius_accessors_partition_independence.py b/tests/parallel/test_1080_radius_accessors_partition_independence.py
similarity index 98%
rename from tests/parallel/test_1078_radius_accessors_partition_independence.py
rename to tests/parallel/test_1080_radius_accessors_partition_independence.py
index f3d5a0433..fdcaf5f6b 100644
--- a/tests/parallel/test_1078_radius_accessors_partition_independence.py
+++ b/tests/parallel/test_1080_radius_accessors_partition_independence.py
@@ -1,6 +1,6 @@
 """The radius ACCESSORS must not depend on how the mesh was partitioned.
 
-``test_1077`` covers the per-cell field cell by cell. This file covers
+``test_1079`` covers the per-cell field cell by cell. This file covers
 ``get_min_radius()``, ``get_max_radius()`` and ``get_mean_radius()``, which
 reduce that field and advertise a global mesh length -- and which were the half
 left partition-dependent when only ``cell_size()`` was fixed: ``get_max_radius()``
diff --git a/tests/test_0010_cell_size_geometry.py b/tests/test_0010_cell_size_geometry.py
index f7c80bd5c..91108be08 100644
--- a/tests/test_0010_cell_size_geometry.py
+++ b/tests/test_0010_cell_size_geometry.py
@@ -7,8 +7,8 @@
 deformed.
 
 Partition independence is NOT tested here -- it cannot be, at one rank count.
-``tests/parallel/test_1077`` compares the field cell by cell against its own
-serial answer, and ``test_1078`` does the same for the three radius accessors.
+``tests/parallel/test_1079`` compares the field cell by cell against its own
+serial answer, and ``test_1080`` does the same for the three radius accessors.
 """
 
 import numpy as np