From a69211354d833f0c2b8d7b5c735eba08b88b88c9 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Fri, 24 Jul 2026 10:08:50 +1000 Subject: [PATCH 1/5] =?UTF-8?q?fix(kdtree/mesh):=20empty=20ranks=20constru?= =?UTF-8?q?ct=20=E2=80=94=20KDTree=20accepts=20zero-point=20clouds,=20nav?= =?UTF-8?q?=20index=20shapes=20empty=20arrays=20(#399)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Any mesh construction where a rank owned zero cells crashed: the navigation kd-tree build passed numpy.array([]) — a 1-D array — into the KDTree memoryview, and the surviving ranks died on the resulting collective divergence (repro: StructuredQuadBox(elementRes=(2,2)) at np3/np4). Two layers, both fixed: - _build_kd_tree_index shapes its point arrays explicitly to (0, cdim) when the cell lists are empty; - KDTree itself accepts an empty (0, dim) cloud: no C++ index is built (taking &points[0][0] would be invalid), and the query methods return honest empties (found=False, infinite distances). A 1-D input now gets a legible error instead of the memoryview's buffer message. The count==0 query guard also removes the same invalid-address hazard for empty QUERY arrays against a populated tree. Regression tests: empty-cloud + 1-D-rejection in test_0101; starved-rank mesh construction at np4 in tests/parallel/test_0700 (domain integral exact). Direct repro verified at np3 and np4. Closes #399 Underworld development team with AI support from Claude Code --- src/underworld3/ckdtree.pyx | 35 ++++++++++++++++++- .../discretisation/discretisation_mesh.py | 10 ++++-- .../test_0700_basic_parallel_operations.py | 15 ++++++++ tests/test_0101_kdtree.py | 23 ++++++++++++ 4 files changed, 80 insertions(+), 3 deletions(-) diff --git a/src/underworld3/ckdtree.pyx b/src/underworld3/ckdtree.pyx index 2a7164fc6..612a752b9 100644 --- a/src/underworld3/ckdtree.pyx +++ b/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( &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( &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') diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py index 7cc36e9ba..c02b29c3e 100644 --- a/src/underworld3/discretisation/discretisation_mesh.py +++ b/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 diff --git a/tests/parallel/test_0700_basic_parallel_operations.py b/tests/parallel/test_0700_basic_parallel_operations.py index 5bae24889..edf92d6dc 100644 --- a/tests/parallel/test_0700_basic_parallel_operations.py +++ b/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}" diff --git a/tests/test_0101_kdtree.py b/tests/test_0101_kdtree.py index c5de2d42e..2588239f8 100644 --- a/tests/test_0101_kdtree.py +++ b/tests/test_0101_kdtree.py @@ -137,3 +137,26 @@ def test_mesh_centroid(res, dim, fill_param): assert np.any(kdpt > index.n) == False, "Some point weren't found. Error" # `find_closest_point` should return index of pts. assert np.allclose(cellid, kdpt), "Point indices weren't as expected." + + +def test_empty_point_cloud(): + """#399: a KDTree over zero points is legitimate — a parallel rank can + own no cells. Construction succeeds and queries report not-found.""" + kd = uw.kdtree.KDTree(np.empty((0, 2))) + assert kd.n == 0 + + query = np.array([[0.5, 0.5]]) + indices, dist_sqr, found = kd.find_closest_point(query) + assert not found[0] + assert np.isinf(dist_sqr[0]) + + n_idx, n_dist = kd.find_closest_n_points(3, query) + assert n_idx.shape == (1, 3) + assert np.isinf(n_dist).all() + + +def test_one_dimensional_input_rejected(): + """#399: numpy.array([]) is 1-D — the constructor must say so legibly + instead of failing in the memoryview cast.""" + with pytest.raises(RuntimeError, match="2-D"): + uw.kdtree.KDTree(np.empty(0)) From 255e2195bbfce9de93a140d7babb27c768f4053d Mon Sep 17 00:00:00 2001 From: lmoresi Date: Fri, 24 Jul 2026 10:08:50 +1000 Subject: [PATCH 2/5] =?UTF-8?q?fix(stats):=20tensor=20stats=20work=20?= =?UTF-8?q?=E2=80=94=20flat-layout=20component=20reads=20+=20raw=20temp=20?= =?UTF-8?q?variable=20(#400)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Two stacked defects made stats() raise for every genuine TENSOR variable: - the Frobenius accumulation indexed the flat component count (d*d) against the structured .array axis of size d (IndexError); - the temporary was created through the enhanced factory, whose __getattr__ refuses underscore-name delegation, so the _scalar_stats call could never be reached (_vector_stats already used _BaseMeshVariable). Components are now read through the flat .data layout and the temp is a _BaseMeshVariable. Regression: identity-tensor Frobenius = sqrt(d) exact, plus a repeat-call/cleanup lifecycle test (tests/test_0102). Closes #400 Underworld development team with AI support from Claude Code --- .../discretisation_mesh_variables.py | 11 +++- tests/test_0102_meshvariable_stats.py | 51 +++++++++++++++++++ 2 files changed, 60 insertions(+), 2 deletions(-) create mode 100644 tests/test_0102_meshvariable_stats.py diff --git a/src/underworld3/discretisation/discretisation_mesh_variables.py b/src/underworld3/discretisation/discretisation_mesh_variables.py index 188fca2a0..208053ffd 100644 --- a/src/underworld3/discretisation/discretisation_mesh_variables.py +++ b/src/underworld3/discretisation/discretisation_mesh_variables.py @@ -3190,16 +3190,23 @@ 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)) + # Components are read through the FLAT layout: a tensor's + # .array is structured (N, d, d), so indexing it with the flat + # component count d*d walked off the axis (issue #400). with uw.synchronised_array_update(): sum_squares = 0.0 for i in range(self.num_components): - component = self.array[:, 0, i].flatten() + component = np.asarray(self.data[:, i]).flatten() sum_squares += component**2 frobenius_var.array[:, 0, 0] = np.sqrt(sum_squares) diff --git a/tests/test_0102_meshvariable_stats.py b/tests/test_0102_meshvariable_stats.py new file mode 100644 index 000000000..1b3e94950 --- /dev/null +++ b/tests/test_0102_meshvariable_stats.py @@ -0,0 +1,51 @@ +"""MeshVariable.stats() over the variable types. + +Regression coverage for: + +- issue #400: ``_tensor_stats`` indexed the flat component count (d*d) + against the structured ``.array`` axis of size d — ``stats()`` on any + genuine TENSOR variable raised IndexError; +- issue #384 lifecycle: the stats temporaries are created and cleaned up + by name — repeated calls must work and leave no variable behind. +""" + +import numpy as np +import pytest +import underworld3 as uw + +pytestmark = [pytest.mark.level_1, pytest.mark.tier_b] + + +@pytest.fixture(scope="module") +def mesh(): + return uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=0.25) + + +def test_tensor_stats_computes_frobenius(mesh): + T = uw.discretisation.MeshVariable( + "T400", mesh, (mesh.dim, mesh.dim), vtype=uw.VarType.TENSOR, degree=1) + # Identity everywhere: ||I||_F = sqrt(d) + with uw.synchronised_array_update(): + T.array[:, :, :] = np.eye(mesh.dim) + + result = T.stats() + assert result["type"] == "tensor" + expected = np.sqrt(mesh.dim) + assert np.isclose(result["mean"], expected, rtol=1e-12) + assert np.isclose(result["max"], expected, rtol=1e-12) + assert np.isclose(result["min"], expected, rtol=1e-12) + + +def test_vector_stats_repeat_and_cleanup(mesh): + V = uw.discretisation.MeshVariable("V400", mesh, mesh.dim, degree=1) + with uw.synchronised_array_update(): + V.array[:, 0, 0] = 3.0 + V.array[:, 0, 1] = 4.0 + + for _ in range(2): # repeated calls: temp is recreated by name each time + result = V.stats() + assert np.isclose(result["magnitude_mean"], 5.0, rtol=1e-12) + + leftovers = [name for name in mesh.vars if name.startswith("_temp_")] + assert leftovers == [], f"stats temporaries left behind: {leftovers}" From d60c9cdab2ced7058ba61688732b94f9fde3f3a2 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Fri, 24 Jul 2026 10:08:50 +1000 Subject: [PATCH 3/5] fix(metrics): mesh_metric_mismatch moments are globally reduced and empty-rank safe (#351) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit A_target/A_mean were built from rank-local sums, so the returned delta moments (rms/max/median_abs) were partition-dependent diagnostics against a per-rank equidistribution target. The sums, the cell count and the moments are now globally reduced (median via gather-to-root + broadcast — a diagnostic at adapt cadence). A rank owning zero cells — constructible since the #399 fix — participates in every reduction with empty contributions instead of raising through the tris-is-None path and desynchronising the collectives. Regression: moments identical across ranks to 1e-14 at np2 and np4 (tests/parallel/test_0750_global_statistics). Closes #351 Underworld development team with AI support from Claude Code --- src/underworld3/meshing/smoothing/metrics.py | 87 +++++++++++-------- tests/parallel/test_0750_global_statistics.py | 22 +++++ 2 files changed, 75 insertions(+), 34 deletions(-) diff --git a/src/underworld3/meshing/smoothing/metrics.py b/src/underworld3/meshing/smoothing/metrics.py index 93183d8d0..d8ce536db 100644 --- a/src/underworld3/meshing/smoothing/metrics.py +++ b/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 / @@ -135,44 +135,48 @@ 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)) + 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: - 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 - 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() + # A rank owning zero cells still participates in every global + # reduction below with empty contributions — raising or returning + # early here would desynchronise the collectives. + 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 +214,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) diff --git a/tests/parallel/test_0750_global_statistics.py b/tests/parallel/test_0750_global_statistics.py index b2822a6e6..89c0b27c7 100644 --- a/tests/parallel/test_0750_global_statistics.py +++ b/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}" + ) From 005ae58a5c32f741b91a6814acfb994392093fbb Mon Sep 17 00:00:00 2001 From: lmoresi Date: Fri, 24 Jul 2026 10:08:50 +1000 Subject: [PATCH 4/5] fix(swarm): points snapshot labels scaled coordinates as metres (#386) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The deprecated swarm.points scales model coordinates by the SI length scale (metres) but labelled the result with mesh.units — kilometre-unit meshes got metre magnitudes under a km label, a 1000x label/value mismatch. Scaled values are now labelled 'meter', the same convention as mesh.X.coords; the unscaled path keeps its previous label. Regression: label + magnitude consistency under an active km-scale model (tests/test_0114). Closes #386 Underworld development team with AI support from Claude Code --- src/underworld3/swarm.py | 10 +++++++++- .../test_0114_swarm_save_coordinate_units.py | 20 +++++++++++++++++++ 2 files changed, 29 insertions(+), 1 deletion(-) diff --git a/src/underworld3/swarm.py b/src/underworld3/swarm.py index 06c522768..da80015ed 100644 --- a/src/underworld3/swarm.py +++ b/src/underworld3/swarm.py @@ -3247,8 +3247,10 @@ 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) @@ -3256,7 +3258,13 @@ def points(self): 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 diff --git a/tests/test_0114_swarm_save_coordinate_units.py b/tests/test_0114_swarm_save_coordinate_units.py index 2743a78f8..d1acbfd9d 100644 --- a/tests/test_0114_swarm_save_coordinate_units.py +++ b/tests/test_0114_swarm_save_coordinate_units.py @@ -102,3 +102,23 @@ def test_save_read_roundtrip_with_scaling(scaled_model, tmp_path): "unit system (physically scaled, outside the model-unit mesh domain)" ) np.testing.assert_allclose(restored, original, rtol=1e-12) + + +def test_points_label_matches_magnitudes(scaled_model): + """#386: the deprecated swarm.points snapshot scales model coordinates + by the SI length scale (metres) — the unit label must be 'meter', not + mesh.units (kilometres here: a 1000x label/value mismatch).""" + mesh = uw.meshing.StructuredQuadBox( + elementRes=(4, 4), minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), + units="km") + swarm = uw.swarm.Swarm(mesh) + swarm.populate(fill_param=2) + + with pytest.warns(DeprecationWarning): + pts = swarm.points + + if getattr(mesh.CoordinateSystem, "_scaled", False): + assert str(pts.units) in ("meter", "metre"), f"label = {pts.units}" + model_coords = np.asarray(swarm._particle_coordinates.data) + scale = mesh.CoordinateSystem._length_scale + assert np.allclose(np.asarray(pts), model_coords * scale) From ad63bc95599a4d49657d1bf7f27d8109f947a0a1 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Fri, 24 Jul 2026 10:50:20 +1000 Subject: [PATCH 5/5] fix: adversarial-review round 1 responses (SYM_TENSOR Frobenius, symmetric raise, sentinel guards, mis-route) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Responses to the public adversarial review of this branch: - Tensor stats read the structured (N, d, d) array view, which mirrors symmetric storage — SYM_TENSOR off-diagonals now count twice as the Frobenius norm requires (the flat read converted the pre-existing crash into a silently under-measured norm). Validated against numpy for non-symmetric TENSOR and Voigt SYM_TENSOR fields. - mesh_metric_mismatch refuses non-simplex meshes rank-symmetrically via mesh.isSimplex (the per-rank tris check let starved ranks sail into the collective reductions while populated ranks raised). The field-valued metric limitation on starved ranks is documented in place — it needs the evaluate layer to be empty-rank safe first (#405). - The three kd-tree sentinel guards use >= (valid indices are 0..n-1), so the empty-tree sentinel trips the sanity error instead of indexing _indexMap out of bounds. - _get_domain_centroids gives starved ranks a huge finite sentinel centroid instead of NaN: gather_data silently strips NaN rows, which de-aligned the rank/row correspondence and made _route_by_nearest_centroid silently migrate particles to a rank owning zero cells (demonstrated live in the review). Finite, not inf — infinities poison kd-tree bounding boxes. - Test hygiene from the review: the swarm.points test asserts _scaled rather than hiding behind it, and drops the deprecated units= mesh kwarg. Deeper findings outside this diff are tracked as #405 (evaluate / points_in_domain / radii / gather_data NaN-stripping — the next empty-rank layer). Underworld development team with AI support from Claude Code --- src/underworld3/ckdtree.pyx | 4 ++- .../discretisation/discretisation_mesh.py | 21 ++++++++++++-- .../discretisation_mesh_variables.py | 19 +++++++----- src/underworld3/meshing/smoothing/metrics.py | 24 +++++++++++++-- tests/test_0102_meshvariable_stats.py | 29 ++++++++++++++++--- .../test_0114_swarm_save_coordinate_units.py | 16 +++++----- 6 files changed, 88 insertions(+), 25 deletions(-) diff --git a/src/underworld3/ckdtree.pyx b/src/underworld3/ckdtree.pyx index 612a752b9..d8eb90b90 100644 --- a/src/underworld3/ckdtree.pyx +++ b/src/underworld3/ckdtree.pyx @@ -624,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" ) diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py index c02b29c3e..30472ea7c 100644 --- a/src/underworld3/discretisation/discretisation_mesh.py +++ b/src/underworld3/discretisation/discretisation_mesh.py @@ -5803,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." ) @@ -5876,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." ) @@ -6279,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 diff --git a/src/underworld3/discretisation/discretisation_mesh_variables.py b/src/underworld3/discretisation/discretisation_mesh_variables.py index 208053ffd..1a71d2c45 100644 --- a/src/underworld3/discretisation/discretisation_mesh_variables.py +++ b/src/underworld3/discretisation/discretisation_mesh_variables.py @@ -3199,15 +3199,20 @@ def _tensor_stats(self): ) try: - # Compute Frobenius norm: ||A||_F = sqrt(sum(A_ij^2)) - # Components are read through the FLAT layout: a tensor's - # .array is structured (N, d, d), so indexing it with the flat - # component count d*d walked off the axis (issue #400). + # 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 = np.asarray(self.data[:, 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 diff --git a/src/underworld3/meshing/smoothing/metrics.py b/src/underworld3/meshing/smoothing/metrics.py index d8ce536db..2aae84e9f 100644 --- a/src/underworld3/meshing/smoothing/metrics.py +++ b/src/underworld3/meshing/smoothing/metrics.py @@ -134,6 +134,16 @@ def mesh_metric_mismatch(mesh, metric, resolution_ratio=None): """ import underworld3 as _uw + # 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") + coords = np.asarray(mesh.X.coords) cStart, cEnd = mesh.dm.getHeightStratum(0) if cEnd > cStart: @@ -155,9 +165,17 @@ def mesh_metric_mismatch(mesh, metric, resolution_ratio=None): metric, centroids)).reshape(-1) rho = np.maximum(rho, 1.0e-12) # guard else: - # A rank owning zero cells still participates in every global - # reduction below with empty contributions — raising or returning - # early here would desynchronise the collectives. + # 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) diff --git a/tests/test_0102_meshvariable_stats.py b/tests/test_0102_meshvariable_stats.py index 1b3e94950..c21cc066f 100644 --- a/tests/test_0102_meshvariable_stats.py +++ b/tests/test_0102_meshvariable_stats.py @@ -25,16 +25,37 @@ def mesh(): def test_tensor_stats_computes_frobenius(mesh): T = uw.discretisation.MeshVariable( "T400", mesh, (mesh.dim, mesh.dim), vtype=uw.VarType.TENSOR, degree=1) - # Identity everywhere: ||I||_F = sqrt(d) + # Non-symmetric, spatially varying — compare against numpy directly. with uw.synchronised_array_update(): - T.array[:, :, :] = np.eye(mesh.dim) + x = np.asarray(T.coords)[:, 0] + T.array[:, 0, 0] = 1.0 + x + T.array[:, 0, 1] = 2.0 + T.array[:, 1, 0] = -3.0 * x + T.array[:, 1, 1] = 0.5 + expected = np.sqrt((np.asarray(T.array) ** 2).sum(axis=(1, 2))) result = T.stats() assert result["type"] == "tensor" - expected = np.sqrt(mesh.dim) + assert np.isclose(result["mean"], expected.mean(), rtol=1e-12) + assert np.isclose(result["max"], expected.max(), rtol=1e-12) + assert np.isclose(result["min"], expected.min(), rtol=1e-12) + + +def test_sym_tensor_stats_counts_off_diagonals_twice(mesh): + # SYM_TENSOR stores Voigt components; the Frobenius norm must count + # the mirrored off-diagonals twice (adversarial-review finding: a flat + # single-count read under-measured ||A||_F silently). + S = uw.discretisation.MeshVariable( + "S400", mesh, (mesh.dim, mesh.dim), vtype=uw.VarType.SYM_TENSOR, degree=1) + with uw.synchronised_array_update(): + S.data[:, 0] = 1.0 # xx + S.data[:, 1] = 2.0 # yy + S.data[:, 2] = 3.0 # xy (mirrored) + + expected = np.sqrt(1.0 + 4.0 + 2.0 * 9.0) + result = S.stats() assert np.isclose(result["mean"], expected, rtol=1e-12) assert np.isclose(result["max"], expected, rtol=1e-12) - assert np.isclose(result["min"], expected, rtol=1e-12) def test_vector_stats_repeat_and_cleanup(mesh): diff --git a/tests/test_0114_swarm_save_coordinate_units.py b/tests/test_0114_swarm_save_coordinate_units.py index d1acbfd9d..9f4cbfac3 100644 --- a/tests/test_0114_swarm_save_coordinate_units.py +++ b/tests/test_0114_swarm_save_coordinate_units.py @@ -109,16 +109,18 @@ def test_points_label_matches_magnitudes(scaled_model): by the SI length scale (metres) — the unit label must be 'meter', not mesh.units (kilometres here: a 1000x label/value mismatch).""" mesh = uw.meshing.StructuredQuadBox( - elementRes=(4, 4), minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), - units="km") + elementRes=(4, 4), minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0)) swarm = uw.swarm.Swarm(mesh) swarm.populate(fill_param=2) with pytest.warns(DeprecationWarning): pts = swarm.points - if getattr(mesh.CoordinateSystem, "_scaled", False): - assert str(pts.units) in ("meter", "metre"), f"label = {pts.units}" - model_coords = np.asarray(swarm._particle_coordinates.data) - scale = mesh.CoordinateSystem._length_scale - assert np.allclose(np.asarray(pts), model_coords * scale) + # The scaled_model fixture activates coordinate scaling — assert it, + # so a fixture change cannot silently vacuate this test. + assert mesh.CoordinateSystem._scaled + + assert str(pts.units) in ("meter", "metre"), f"label = {pts.units}" + model_coords = np.asarray(swarm._particle_coordinates.data) + scale = mesh.CoordinateSystem._length_scale + assert np.allclose(np.asarray(pts), model_coords * scale)