fix: make stabilization cell sizes local to each cell (#687)

Adapt only the mesh-size correction from lmoresi's 68e545f on feature/navier-stokes-supg; do not import Navier-Stokes or other branch changes. Cache _radii_own from current DM vertex coordinates and use it for mesh.cell_size(). Preserve the legacy kd-tree radius arrays and global timestep/mesh-motion consumers.

Use coordinate-section offsets and the full vertex stratum so the own-cell RMS definition also handles hexahedra, which have eight vertices but six faces. Correct the field documentation and Nitsche mechanism tests for the new definition; retain physical solve tolerances and use the exact nearest-centroid <= own-centroid ordering instead of an arbitrary approximate-equality tolerance.

Add a first-failing independent geometry/deformation regression for triangles, tetrahedra, quadrilaterals and hexahedra plus a regular-square analytical control. Before: four failures in serial and on eight ranks. After rebuild: 21 passed/one expected skip serial (22.90 s), 22 passed on eight ranks (40.45 s), covering Nitsche solves, radius accessors, frozen PC2 migration and memory/disk snapshots. Own-cell geometry error is zero in these tests; style and whitespace gates pass.

Author: gthyagi

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 own-cell characteristic lengths ``self._radii_own`` 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._radii_own`` 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._radii_own).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_r_own = 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_r_own[cell] = np.sqrt(np.mean(np.sum(delta ** 2, axis=1)))
+
+        self._radii_own = cell_r_own
         return cell_min_r, cell_r, centroids, cell_length
 
     # ==========

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._radii_own``), 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._radii_own).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._radii_own (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._radii_own).reshape(-1)
 
     # not stale: the field changed with the geometry SOMEWHERE (global OR) ...
     nb = min(h_after.shape[0], h_before.shape[0])