Make mesh.cell_size() local and partition independent

docs/developer/CHANGELOG.md

@@ -506,6 +506,11 @@ in `Stokes_Constrained` (#224), then made parallel-correct.
 minimum radius, restoring correct stiffness on graded and adapted meshes
 (#275).
 
+- The local size now comes from each cell's own geometry instead of a kd-tree
+  over the centroids held by the current MPI rank. The old field changed at
+  partition boundaries and moved the default ``local_h=True`` Nitsche velocity
+  answer by 6.6e-3 between rank counts; the replacement is cell-by-cell
+  identical from one to eight ranks (#569, #687).
 - `mesh.boundary_slip` API with `BoundingSurface` objects for boundary
   tangent-slip (#225); `Surface.influence_function` respects finite edges
   (#241).

src/underworld3/discretisation/discretisation_mesh.py

@@ -3224,7 +3224,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
+        RMS distance of its vertices from their own centroid). This is a
+        purely cell-local quantity, 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,34 +3291,19 @@ 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
+        Uses the cell-geometry characteristic lengths ``self._cell_radii`` computed by
         :meth:`_get_mesh_sizes` on the *current* geometry. A degree-0
-        discontinuous variable's local DOFs and ``self._radii`` are BOTH
+        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.
 
         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``."""
-        # TODO(BUG): this field is PARTITION-DEPENDENT, and so therefore is the
-        # Nitsche penalty gamma*mu/h that consumes it (local_h=True, the default).
-        # Not the indexing here — the values. `_get_mesh_sizes` measures a cell by
-        # the distance from its vertices to the NEAREST CENTROID in a kd-tree built
-        # from THIS RANK's centroids, so near a partition seam the nearest centroid
-        # may simply be absent. Measured on Annulus(cellSize=0.12): the field's sum
-        # is 26.0822 at np=1, 26.1211 at np=2 and 26.1386 at np=4, and its max moves
-        # at np=4. End to end that is 6.6e-03 in the velocity of a Nitsche free-slip
-        # annulus and it does NOT shrink with solver tolerance.
-        # This is a DIFFERENT defect from the boundary normal fixed for #564 (which
-        # is now clean: the same solve with local_h=False agrees to 3.6e-10 at
-        # np=1..4). It is the local h that is left, and it also reaches every other
-        # consumer of `cell_size()`. Not fixed here because `_get_mesh_sizes` also
-        # feeds `get_min_radius`, the adaptivity metrics and the free-surface
-        # 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)
+        # Own-cell radii fix #687 without changing the legacy kd-tree radii
+        # used by global timestep estimates, adaptivity, and mesh relaxation.
+        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
@@ -6895,8 +6881,11 @@ def _eval_use_robust_location(self) -> bool:
 
     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
+        Cache own-cell radii for cell_size and return legacy kd-tree radii.
+
+        Own-cell sizes use current DM vertices, so neither partition-local
+        neighbours nor stale coordinate views affect stabilization (#687).
+        Legacy radii remain unchanged for their other consumers.
         """
 
         centroids = self._get_coords_for_basis(0, False)
@@ -6909,6 +6898,9 @@ def _get_mesh_sizes(self, verbose=False):
         cell_length = np.empty(centroids.shape[0])
         cell_min_r = np.empty(centroids.shape[0])
         cell_r = np.empty(centroids.shape[0])
+        cell_radii = np.empty(centroids.shape[0])
+        coordinate_section = self.dm.getCoordinateDM().getLocalSection()
+        vertex_coordinates = self.dm.getCoordinatesLocal().array
 
         for cell in range(cEnd - cStart):
             cell_num_points = self.dm.getConeSize(cell)
@@ -6922,6 +6914,16 @@ def _get_mesh_sizes(self, verbose=False):
             cell_r[cell] = np.sqrt(distsq.mean())
             cell_min_r[cell] = np.sqrt(distsq.min())
 
+            # A hex has six faces but eight vertices: select the vertex
+            # stratum, not a cone-sized suffix of its transitive closure.
+            closure = self.dm.getTransitiveClosure(cStart + cell)[0]
+            vertices = closure[(closure >= pStart) & (closure < pEnd)]
+            offsets = np.array([coordinate_section.getOffset(int(v)) for v in vertices])
+            own_coords = vertex_coordinates[offsets[:, None] + np.arange(self.cdim)]
+            delta = own_coords - own_coords.mean(axis=0)
+            cell_radii[cell] = np.sqrt(np.mean(np.sum(delta ** 2, axis=1)))
+
+        self._cell_radii = cell_radii
         return cell_min_r, cell_r, centroids, cell_length
 
     # ==========

tests/parallel/test_1069_boundary_normal_parallel.py

@@ -333,14 +333,12 @@ 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
-    TODO(BUG) on ``Mesh._assemble_cell_size``.
+    This test now leaves ``local_h`` at its default ``True``. Before #569/#687,
+    doing so mixed the boundary-normal regression with a second partition-dependent
+    input from ``mesh.cell_size()``; this test therefore had to disable the public
+    default. The cell-local geometric size is now partition independent, so retaining
+    the default jointly guards the normal assembly and the Nitsche penalty path users
+    actually run.
     """
     RI, RO = 0.5, 1.0
     mesh = uw.meshing.Annulus(radiusInner=RI, radiusOuter=RO, cellSize=0.12, qdegree=3)
@@ -357,7 +355,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()

tests/parallel/test_1077_cell_size_partition_independence.py

@@ -0,0 +1,47 @@
+"""Rank-count regression for the cell-local stabilization length.
+
+The parallel result is compared with a fresh single-rank run on the same Gmsh
+mesh.  This checks the complete cell geometry table, not merely a reduction or
+a within-rank geometric identity.
+"""
+
+import numpy as np
+import pytest
+
+import underworld3 as uw
+
+from serial_reference import emit, mesh_fingerprint, serial_reference
+
+
+pytestmark = [pytest.mark.mpi(min_size=2), pytest.mark.timeout(300)]
+
+
+def _cell_geometry_table():
+    mesh = uw.meshing.UnstructuredSimplexBox(cellSize=0.12, qdegree=2)
+    mesh.cell_size()
+
+    local = np.column_stack((mesh._centroids, mesh._cell_radii))
+    local = local[mesh._get_owned_cells_mask()]
+    gathered = uw.mpi.comm.allgather(local)
+    table = np.vstack(gathered)
+    order = np.lexsort(tuple(table[:, axis] for axis in reversed(range(mesh.dim))))
+    return table[order].reshape(-1), mesh_fingerprint(mesh)
+
+
+def test_cell_size_matches_single_rank_cell_by_cell():
+    values, fingerprint = _cell_geometry_table()
+    reference = serial_reference(__file__, "cell_size")
+
+    assert int(fingerprint[0]) == int(reference["fingerprint"][0])
+    assert np.isclose(
+        fingerprint[1], reference["fingerprint"][1], rtol=1.0e-12, atol=0.0
+    )
+
+    expected = np.asarray(reference["values"])
+    assert values.shape == expected.shape
+    np.testing.assert_allclose(values, expected, rtol=0.0, atol=1.0e-14)
+
+
+if __name__ == "__main__":
+    _values, _fingerprint = _cell_geometry_table()
+    emit(_values, _fingerprint)

tests/test_0010_cell_size_geometry.py

@@ -0,0 +1,66 @@
+"""Issue #687: cell_size is an own-cell geometric quantity, including after deform.
+
+The independent oracle reads vertex coordinates through the coordinate section;
+it does not use the mesh's cached radii or centroid kd-tree. Run serial and MPI.
+"""
+
+import numpy as np
+import pytest
+
+import underworld3 as uw
+
+pytestmark = [pytest.mark.level_1, pytest.mark.tier_b]
+
+
+def _vertex_rms(mesh):
+    dm = mesh.dm
+    section = dm.getCoordinateDM().getLocalSection()
+    coordinates = dm.getCoordinatesLocal().array
+    start, end = dm.getHeightStratum(0)
+    first_vertex, last_vertex = dm.getDepthStratum(0)
+    radii = []
+    for cell in range(start, end):
+        vertices = [int(point) for point in dm.getTransitiveClosure(cell)[0]
+                    if first_vertex <= point < last_vertex]
+        points = np.array([coordinates[section.getOffset(v):section.getOffset(v) + mesh.cdim]
+                           for v in vertices])
+        radii.append(np.sqrt(np.mean(np.sum((points - points.mean(axis=0)) ** 2, axis=1))))
+    return np.asarray(radii)
+
+
+@pytest.mark.parametrize("dim", [2, 3])
+@pytest.mark.parametrize("simplex", [True, False], ids=["simplex", "tensor"])
+def test_cell_size_matches_own_vertices_and_tracks_deform(dim, simplex):
+    geometry = dict(minCoords=(0.0,) * dim, maxCoords=(1.0,) * dim, qdegree=3)
+    mesh = (uw.meshing.UnstructuredSimplexBox(**geometry, cellSize=0.25, regular=False)
+            if simplex else uw.meshing.StructuredQuadBox(**geometry, elementRes=(4,) * dim))
+    mesh.cell_size()
+    field = mesh._cell_size_variable
+    errors = []
+    for phase in ("initial", "deformed"):
+        if phase == "deformed":
+            coordinates = np.array(mesh.X.coords)
+            coordinates[:, 0] = 1.7 * coordinates[:, 0] + 0.2 * coordinates[:, 1]
+            mesh.deform(coordinates)
+        expected = _vertex_rms(mesh)
+        actual = np.asarray(field.array[:, 0, 0])
+        shapes_match = actual.shape == expected.shape
+        assert all(uw.mpi.comm.allgather(shapes_match)), (actual.shape, expected.shape)
+        local_error = float(np.abs(actual - expected).max(initial=0.0))
+        error = max(uw.mpi.comm.allgather(local_error))
+        errors.append(error)
+        uw.pprint(f"CELL_SIZE_GEOMETRY dim={dim} simplex={simplex} phase={phase} "
+                  f"ranks={uw.mpi.size} max_error={error:.12g}")
+    assert max(errors) < 1e-12, errors
+
+
+def test_regular_square_cell_size_keeps_global_radius():
+    mesh = uw.meshing.StructuredQuadBox(elementRes=(4, 4), qdegree=2)
+    legacy = np.array(mesh._radii)
+    global_radius = mesh.get_min_radius()
+    mesh.cell_size()
+    expected = np.sqrt(2.0) / 8.0
+    error = float(np.abs(np.asarray(mesh._cell_size_variable.array) - expected).max(initial=0.0))
+    assert max(uw.mpi.comm.allgather(error)) < 1e-12
+    assert global_radius == pytest.approx(expected, rel=1e-12)
+    assert all(uw.mpi.comm.allgather(np.array_equal(mesh._radii, legacy)))

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."""
+    own-vertex RMS 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)
+    field = np.asarray(mesh._cell_size_variable.array[:, 0, 0]).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)
@@ -168,22 +168,23 @@ def test_cell_size_is_local_per_cell():
     gfmin, gfmax = _gmin(field), _gmax(field)
     assert gfmax / gfmin > 3.0
 
-    # the global scalar that global-h would use is just the minimum cell size
-    assert np.isclose(mesh.get_min_radius(), gfmin, rtol=1e-6)
+    # The unchanged nearest-centroid minimum cannot exceed the own-cell
+    # minimum; these are no longer the same definition on an irregular mesh.
+    assert 0.0 < mesh.get_min_radius() <= gfmin * (1.0 + 1e-12)
 
 
 def test_local_h_at_coarse_freeslip_boundary_exceeds_global_min():
     """At the COARSE Top free-slip boundary the LOCAL penalty size is many
     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
+    # from the field / _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_size_variable.array[:, 0, 0]).reshape(-1)
 
     near_top = cen[:, 1] > 0.85          # cells adjacent to the Top free-slip edge
     h_top = radii[near_top]
@@ -204,14 +205,14 @@ def test_cell_size_tracks_deformation():
     the Nitsche mis-scaling on the free surface."""
     mesh = _graded_box()
     _ = mesh.cell_size()
-    h_before = mesh._cell_size_variable.data[:, 0].copy()
+    h_before = np.array(mesh._cell_size_variable.array[:, 0, 0])
 
     X = np.asarray(mesh.X.coords).copy()
     moved = mesh.deform(_box_wobble(X, amp=0.04))
     assert moved  # geometry actually changed
 
-    h_after = mesh._cell_size_variable.data[:, 0].copy()
-    radii_after = np.asarray(mesh._radii).reshape(-1)
+    h_after = np.array(mesh._cell_size_variable.array[:, 0, 0])
+    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])