Take the cell radius from PETSc, and fix the radius accessors with it (#694)

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)

src/underworld3/discretisation/discretisation_mesh.py

@@ -2840,12 +2840,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]:
@@ -3234,8 +3232,8 @@ def cell_size(self):
 
         Returns the ``.sym`` of a cell-constant (degree-0, discontinuous)
         scalar MeshVariable holding each cell's characteristic length (the
-        RMS distance of its vertices from their own centroid). This is a
-        purely cell-local quantity, independent of the MPI partition. 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
@@ -3301,27 +3299,45 @@ def _refresh():
     def _assemble_cell_size(self, var):
         """Fill ``var`` (degree-0 scalar) with each cell's characteristic size.
 
-        Uses the cell-geometry characteristic lengths ``self._cell_radii`` computed by
-        :meth:`_get_mesh_sizes` on the *current* geometry. A degree-0
+        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.
 
-        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``."""
-        # Own-cell radii fix #687 without changing the legacy kd-tree radii
-        # used by global timestep estimates, adaptivity, and mesh relaxation.
+        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).
+        # 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 to fill on this rank.
-        if radii.size == 0 or var.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]
+
+        # `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]
 
     @property
     def Gamma_P1(self):
@@ -6889,56 +6905,44 @@ def _eval_use_robust_location(self) -> bool:
         """
         return (uw.mpi.size > 1) and (self._location_capability() != "none")
 
-    def _get_mesh_sizes(self, verbose=False):
-        """
-        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.
-        """
-
+    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])
-        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)
-            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())
-
-            # 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
+        return radii, centroids
 
     # ==========
 
-    # Deprecated in favour of _get_mesh_sizes (above)
+    # Deprecated in favour of _get_cell_radii (above)
     def _get_mesh_centroids(self):
         """
         Obtain and cache the (local) mesh centroids using underworld swarm technology.
@@ -7028,7 +7032,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)
@@ -7050,7 +7054,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)
@@ -7069,15 +7073,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:

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

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)

src/underworld3/systems/free_surface.py

@@ -121,7 +121,7 @@ def __init__(
         driving_buoyancy=None,
         smooth_length=0.0,
         mass="lumped",
-        max_surface_cfl=0.5,
+        max_surface_cfl=0.437,
         tangent_advect=None,
         tangent_spectral_modes=0,
         surface_mask=None,

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).

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)