Empty-rank mesh construction, tensor stats, global metric moments, swarm.points label

src/underworld3/ckdtree.pyx

@@ -91,6 +91,16 @@ cdef class KDTree:
         import underworld3.function.unit_conversion as unit_conv
         self.coord_units = unit_conv.get_units(points_input) if unit_conv.has_units(points_input) else None
 
+        # A legible error beats the memoryview's "Buffer has wrong number
+        # of dimensions" — the common trigger is numpy.array([]) from an
+        # empty point list, which is 1-D (issue #399).
+        if points_input.ndim != 2:
+            raise RuntimeError(
+                f"KDTree points must be a 2-D (n_points, dim) array, "
+                f"got shape {points_input.shape}. (An empty point set must "
+                f"still be shaped (0, dim).)"
+            )
+
         # Extract raw numpy array for C++ interface
         cdef const double[:,::1] points
         if unit_conv.has_units(points_input):
@@ -101,7 +111,13 @@ cdef class KDTree:
         if points.shape[1] not in (2,3):
             raise RuntimeError(f"Provided points array dimensionality must be 2 or 3, not {points.shape[1]}.")
         self.points = points
-        self.index = new KDTree_Interface(<const double *> &points[0][0], points.shape[0], points.shape[1])
+        # An empty point cloud is legitimate — a parallel rank can own zero
+        # cells (issue #399). No C++ index is built (taking &points[0][0]
+        # would be invalid); the query methods return honest empties.
+        if points.shape[0] > 0:
+            self.index = new KDTree_Interface(<const double *> &points[0][0], points.shape[0], points.shape[1])
+        else:
+            self.index = NULL
 
         global _live_instances, _total_constructed
         _live_instances += 1
@@ -190,6 +206,8 @@ cdef class KDTree:
         """
         Build the kd-tree index.
         """
+        if self.index is NULL:
+            return
         self.index.build_index()
 
     def kdtree_points(self):
@@ -235,6 +253,14 @@ cdef class KDTree:
         dist_sqr = np.empty(count, dtype=np.float64, order='C')
         found    = np.empty(count, dtype=np.bool_,   order='C')
 
+        # Empty index (a rank owning zero points, issue #399): nothing can
+        # be found — report that honestly for every query.
+        if self.index is NULL or count == 0:
+            indices[...] = 0
+            dist_sqr[...] = np.inf
+            found[...] = False
+            return indices, dist_sqr, found
+
         cdef long unsigned int[::1]  c_indices = indices
         cdef            double[::1] c_dist_sqr = dist_sqr
         cdef              bool[::1]    c_found = found
@@ -281,6 +307,13 @@ cdef class KDTree:
         n_indices  = np.empty((coords.shape[0], nCount), dtype=np.uint64,  order='C')
         n_dist_sqr = np.empty((coords.shape[0], nCount), dtype=np.float64,  order='C')
 
+        # Empty index (a rank owning zero points, issue #399): no
+        # neighbours exist — infinite distances, sentinel indices.
+        if self.index is NULL or nInput == 0:
+            n_indices[...] = 0
+            n_dist_sqr[...] = np.inf
+            return n_indices, n_dist_sqr
+
         indices  = np.empty(nCount, dtype=np.uint64,  order='C')
         dist_sqr = np.empty(nCount, dtype=np.float64, order='C')
 
@@ -591,7 +624,9 @@ cdef class KDTree:
         # Note: query() returns sqr_dists=True by default, and we use the converted coords
         distance_n, closest_n = self.query(coords, k=nnn)
 
-        if np.any(closest_n > self.n):
+        # valid indices are 0..n-1; the empty-tree sentinel (0 with n=0)
+        # must trip this guard, so the comparison is >= (issue #399).
+        if np.any(closest_n >= self.n):
             raise RuntimeError(
                 "Error in rbf_interpolator_local_from_kdtree - a nearest neighbour wasn't found"
             )

src/underworld3/discretisation/discretisation_mesh.py

@@ -5171,7 +5171,12 @@ def _build_kd_tree_index(self):
             control_points_list.append(cell_centroid)
             control_points_cell_list.append(cell_id)
 
-        self._indexCoords = numpy.array(control_points_list)
+        # A rank can own zero cells (small mesh / imbalanced partition):
+        # shape the point arrays explicitly so an empty list becomes a
+        # well-formed (0, cdim) array rather than a 1-D numpy.array([]),
+        # which crashed the KDTree construction (issue #399).
+        self._indexCoords = numpy.array(
+            control_points_list, dtype=numpy.float64).reshape(-1, self.cdim)
         self._index = uw.kdtree.KDTree(self._indexCoords)
         # self._index.build_index()
         self._indexMap = numpy.array(control_points_cell_list, dtype=numpy.int64)
@@ -5182,7 +5187,8 @@ def _build_kd_tree_index(self):
         # We keep _nav_centroids separate from _centroids (which is
         # the main-DM cell centroids set in __init__) so the FE-side
         # ``_centroids`` semantics are unchanged on manifold meshes.
-        self._nav_centroids = numpy.array(centroids_list)
+        self._nav_centroids = numpy.array(
+            centroids_list, dtype=numpy.float64).reshape(-1, self.cdim)
         self._centroid_index = uw.kdtree.KDTree(self._nav_centroids)
 
         return
@@ -5797,7 +5803,9 @@ def get_closest_cells(self, coords: numpy.ndarray) -> numpy.ndarray:
 
         if len(model_coords) > 0:
             dist, closest_points = self._index.query(model_coords, k=1, sqr_dists=False)
-            if np.any(closest_points > self._index.n):
+            # >= : valid indices are 0..n-1, and the empty-tree sentinel
+            # (0 with n=0) must trip this guard, not index _indexMap (#399).
+            if np.any(closest_points >= self._index.n):
                 raise RuntimeError(
                     "An error was encountered attempting to find the closest cells to the provided coordinates."
                 )
@@ -5870,7 +5878,9 @@ def _get_closest_local_cells_internal(
 
         if len(coords) > 0:
             dist, closest_points = self._index.query(coords, k=1, sqr_dists=False)
-            if np.any(closest_points > self._index.n):
+            # >= : valid indices are 0..n-1, and the empty-tree sentinel
+            # (0 with n=0) must trip this guard, not index _indexMap (#399).
+            if np.any(closest_points >= self._index.n):
                 raise RuntimeError(
                     "An error was encountered attempting to find the closest cells to the provided coordinates."
                 )
@@ -6273,7 +6283,18 @@ def _get_domain_centroids(self):
         import numpy as np
         from underworld3.utilities import gather_data
 
-        domain_centroid = self._centroids.mean(axis=0)
+        # A rank owning zero cells has no centroid; mean() of the empty
+        # array is NaN, and gather_data silently STRIPS NaN rows — the
+        # gathered table's row index then no longer equals rank, and
+        # _route_by_nearest_centroid mis-routes particles to the wrong
+        # rank (issue #399 review). A huge FINITE sentinel keeps the row
+        # (row == rank) while a nearest-centroid search can never select
+        # it, so starved ranks correctly receive no particles. (Finite,
+        # not inf: infinities poison the kd-tree's bounding boxes.)
+        if self._centroids.shape[0] > 0:
+            domain_centroid = self._centroids.mean(axis=0)
+        else:
+            domain_centroid = np.full(self.cdim, 1.0e30)
         all_centroids = gather_data(domain_centroid, bcast=True).reshape(-1, self.cdim)
         return all_centroids
 

src/underworld3/discretisation/discretisation_mesh_variables.py

@@ -3190,17 +3190,29 @@ def _tensor_stats(self):
                 f"Cannot compute tensor stats: a variable named '{temp_name}' "
                 "already exists on this mesh (reserved as a stats temporary)."
             )
-        frobenius_var = uw.discretisation.MeshVariable(
+        # _BaseMeshVariable, matching _vector_stats: the enhanced wrapper's
+        # __getattr__ refuses underscore-name delegation, so the
+        # _scalar_stats call below would raise AttributeError through it
+        # (second latent defect under issue #400).
+        frobenius_var = _BaseMeshVariable(
             temp_name, self.mesh, 1, degree=self.degree
         )
 
         try:
-            # Compute Frobenius norm: ||A||_F = sqrt(sum(A_ij^2))
+            # Compute Frobenius norm: ||A||_F = sqrt(sum_ij(A_ij^2))
+            # Structured (N, d, d) reads (issue #400): correct for full
+            # tensors AND for symmetric storage — the .array view mirrors
+            # the Voigt components, so off-diagonals are counted twice as
+            # the Frobenius sum requires. (Flat component reads counted
+            # each Voigt entry once and under-measured SYM_TENSOR norms;
+            # the original structured read with the FLAT component count
+            # walked off the axis.)
             with uw.synchronised_array_update():
+                arr = np.asarray(self.array)
                 sum_squares = 0.0
-                for i in range(self.num_components):
-                    component = self.array[:, 0, i].flatten()
-                    sum_squares += component**2
+                for i in range(self.shape[0]):
+                    for j in range(self.shape[1]):
+                        sum_squares = sum_squares + arr[:, i, j] ** 2
                 frobenius_var.array[:, 0, 0] = np.sqrt(sum_squares)
 
             # Get scalar stats on Frobenius norm

src/underworld3/meshing/smoothing/metrics.py

@@ -9,7 +9,7 @@
 
 import underworld3 as uw
 
-from .graph import _tri_cells, _signed_areas, _global_sum
+from .graph import _tri_cells, _signed_areas, _global_sum, _global_max
 
 
 # Named adaptation strategies (off / vlow / low / med / high /
@@ -134,45 +134,67 @@ def mesh_metric_mismatch(mesh, metric, resolution_ratio=None):
     """
     import underworld3 as _uw
 
-    coords = np.asarray(mesh.X.coords)
-    if mesh.cdim == 2:
-        tris = _tri_cells(mesh.dm)
-        A_actual = None if tris is None else np.abs(
-            _signed_areas(coords, tris))
-    else:
-        from .graph import _tet_cells, _signed_volumes
-        tris = _tet_cells(mesh.dm)
-        A_actual = None if tris is None else np.abs(
-            _signed_volumes(coords, tris))
-    if tris is None:
+    # Rank-SYMMETRIC refusal for non-simplex meshes: mesh.isSimplex is a
+    # rank-uniform property, so every rank raises together. (The per-rank
+    # tris-is-None check below cannot serve — a starved rank has no cells
+    # to inspect and would sail into the collective reductions while the
+    # populated ranks raise; issue #351 review.)
+    if not mesh.isSimplex:
         raise NotImplementedError(
             "mesh_metric_mismatch: simplex (triangle/tetrahedral) mesh "
             "required")
-    centroids = coords[tris].mean(axis=1)
-    rho = np.asarray(_uw.function.evaluate(
-        metric, centroids)).reshape(-1)
-    rho = np.maximum(rho, 1.0e-12)   # guard
-    inv_rho = 1.0 / rho
-    # TODO(BUG): A_target (and A_mean below) are built from rank-LOCAL
-    # sums, so under MPI the delta moments returned here (rms / max /
-    # median_abs) are partition-dependent diagnostics measured against a
-    # per-rank equidistribution target, not the global one the docstring
-    # describes. The skip signal consumed by smooth_mesh_interior
-    # (alignment / misalignment) is NOT affected — it is computed from
-    # globally reduced moment sums further down. Fix: reduce
-    # A_actual.sum(), inv_rho.sum() and the cell count globally before
-    # forming A_target / A_mean (delta stays a local array; its moments
-    # then need global reductions too).
-    A_target = A_actual.sum() * inv_rho / inv_rho.sum()
+
+    coords = np.asarray(mesh.X.coords)
+    cStart, cEnd = mesh.dm.getHeightStratum(0)
+    if cEnd > cStart:
+        if mesh.cdim == 2:
+            tris = _tri_cells(mesh.dm)
+            A_actual = None if tris is None else np.abs(
+                _signed_areas(coords, tris))
+        else:
+            from .graph import _tet_cells, _signed_volumes
+            tris = _tet_cells(mesh.dm)
+            A_actual = None if tris is None else np.abs(
+                _signed_volumes(coords, tris))
+        if tris is None:
+            raise NotImplementedError(
+                "mesh_metric_mismatch: simplex (triangle/tetrahedral) mesh "
+                "required")
+        centroids = coords[tris].mean(axis=1)
+        rho = np.asarray(_uw.function.evaluate(
+            metric, centroids)).reshape(-1)
+        rho = np.maximum(rho, 1.0e-12)   # guard
+    else:
+        # A rank owning zero cells participates in the global reductions
+        # below with empty contributions — raising or returning early here
+        # would desynchronise the collectives.
+        #
+        # KNOWN LIMIT: this skips uw.function.evaluate, which is itself
+        # collective for metrics containing MESH-VARIABLE data — a starved
+        # rank then deadlocks the populated ranks inside evaluate. Full
+        # starved-rank support for field-valued metrics needs the
+        # evaluate/points_in_domain layer to be empty-rank safe first
+        # (tracked with the #399 follow-up issue). Analytic (pure-sympy)
+        # metrics are fine: their evaluation is rank-local.
+        A_actual = np.empty(0)
+        rho = np.empty(0)
+    inv_rho = 1.0 / rho if rho.size else np.empty(0)
+    # Global equidistribution target (issue #351): the sums and the cell
+    # count are reduced across ranks — rank-local sums measured each
+    # rank's cells against a per-rank target, so the returned moments
+    # were partition-dependent and inconsistent with the docstring.
+    sum_A = _global_sum(A_actual.sum())
+    sum_inv_rho = _global_sum(inv_rho.sum())
+    A_target = sum_A * inv_rho / sum_inv_rho
     if resolution_ratio is not None:
         R = float(resolution_ratio)
-        A_mean = A_actual.sum() / len(A_actual)
+        A_mean = sum_A / _global_sum(A_actual.size)
         # Clip target areas to the mover's achievable band
         # [A_mean/R², A_mean·R²] (h in [h0/R, h0·R] ⇒
         # A in [h0²/R², h0²·R²] = [A_mean/R², A_mean·R²]).
         A_target = np.clip(A_target, A_mean / R ** 2,
                            A_mean * R ** 2)
-    delta = 0.5 * np.log(A_actual / A_target)
+    delta = 0.5 * np.log(A_actual / A_target) if A_actual.size else np.empty(0)
     abs_delta = np.abs(delta)
 
     # Alignment — Pearson r of log(1/A_c) with log(ρ_c).
@@ -210,9 +232,24 @@ def mesh_metric_mismatch(mesh, metric, resolution_ratio=None):
     misalignment = float(
         np.sqrt(max(0.0, 1.0 - max(0.0, alignment) ** 2)))
 
-    return dict(rms=float(np.sqrt(np.mean(delta ** 2))),
-                max=float(abs_delta.max()),
-                median_abs=float(np.median(abs_delta)),
+    # Global moments of |δ| (issue #351): every rank reports the same
+    # diagnostics. The exact median needs the values in one place — this
+    # is a diagnostic at adapt cadence, so gather |δ| to rank 0 and
+    # broadcast the scalar.
+    n_cells = _global_sum(delta.size)
+    rms = float(np.sqrt(_global_sum((delta ** 2).sum()) / n_cells))
+    delta_max = _global_max(abs_delta.max() if abs_delta.size else 0.0)
+    if _uw.mpi.size > 1:
+        gathered = _uw.mpi.comm.gather(abs_delta, root=0)
+        median_abs = (float(np.median(np.concatenate(gathered)))
+                      if _uw.mpi.rank == 0 else None)
+        median_abs = _uw.mpi.comm.bcast(median_abs, root=0)
+    else:
+        median_abs = float(np.median(abs_delta))
+
+    return dict(rms=rms,
+                max=float(delta_max),
+                median_abs=median_abs,
                 alignment=alignment,
                 misalignment=misalignment)
 

src/underworld3/swarm.py

@@ -3247,16 +3247,24 @@ def points(self):
         # Scale model coordinates to physical coordinates
         if hasattr(self.mesh.CoordinateSystem, "_scaled") and self.mesh.CoordinateSystem._scaled:
             coords = model_coords * self.mesh.CoordinateSystem._length_scale
+            scaled_to_si = True
         else:
             coords = model_coords
+            scaled_to_si = False
 
         coords.flags.writeable = False
         coords = coords.view(_ReadOnlyCoordinateSnapshot)
 
         if hasattr(self.mesh, "units") and self.mesh.units is not None:
             from underworld3.utilities.unit_aware_array import UnitAwareArray
 
-            return UnitAwareArray(coords, units=self.mesh.units)
+            # The _length_scale factor converts model coordinates to SI
+            # metres, so scaled values are labelled "meter" — the same
+            # convention as mesh.X.coords. Labelling metre magnitudes with
+            # mesh.units (e.g. kilometres) was a 1000x label/value
+            # mismatch (issue #386).
+            return UnitAwareArray(
+                coords, units="meter" if scaled_to_si else self.mesh.units)
 
         return coords
 

tests/parallel/test_0700_basic_parallel_operations.py

@@ -328,3 +328,18 @@ def test_that_parallel_tests_end():
     import sys
 
     sys.exit(pytest.main([__file__, "-v", "--with-mpi"]))
+
+
+@pytest.mark.mpi(min_size=4)
+def test_empty_rank_mesh_construction():
+    """#399: mesh construction must survive ranks that own zero cells.
+
+    A 2x2 quad box at np4 starves at least one rank; the navigation
+    kd-tree build then received a 1-D empty array and crashed (with the
+    other ranks dying on the resulting collective divergence).
+    """
+    mesh = uw.meshing.StructuredQuadBox(elementRes=(2, 2))
+    uw.discretisation.MeshVariable("T399", mesh, 1, degree=1)
+
+    area = float(uw.maths.Integral(mesh, 1.0).evaluate())
+    assert np.isclose(area, 1.0, rtol=1e-10), f"domain area = {area}"

tests/parallel/test_0750_global_statistics.py

@@ -389,3 +389,25 @@ def test_empty_array_handling():
     import sys
 
     sys.exit(pytest.main([__file__, "-v", "--with-mpi"]))
+
+
+@pytest.mark.mpi(min_size=2)
+def test_metric_mismatch_moments_partition_independent():
+    """#351: mesh_metric_mismatch must return the SAME delta moments on
+    every rank — the sums feeding A_target/A_mean and the moments
+    themselves are globally reduced (rank-local sums made the diagnostics
+    partition-dependent)."""
+    import sympy
+    from underworld3.meshing.smoothing.metrics import mesh_metric_mismatch
+
+    mesh = uw.meshing.UnstructuredSimplexBox(
+        minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=0.1)
+    x, y = mesh.X
+
+    result = mesh_metric_mismatch(mesh, 1.0 + 10.0 * x)
+
+    for key in ("rms", "max", "median_abs", "alignment", "misalignment"):
+        gathered = uw.mpi.comm.allgather(result[key])
+        assert max(gathered) - min(gathered) < 1e-14, (
+            f"{key} differs across ranks: {gathered}"
+        )