From 5e748b73dcf8027cf2e34511e8f1ad2aaf6ec3dc Mon Sep 17 00:00:00 2001 From: lmoresi Date: Mon, 25 May 2026 10:32:44 +1000 Subject: [PATCH 1/2] mesh: add get_mean_radius() + @collective_operation on radii accessors MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Adds a parallel-safe mean-cell-radius accessor mirroring the existing get_min_radius / get_max_radius pattern (MPI allreduce of local sum and count). Together the three methods form the canonical "characteristic mesh length" API. Motivation: user code (e.g. an adaptive mesh harness) reached for mesh._radii.mean() to set a smoothing length default, which is rank-local: different ranks compute different means → different gradient_smoothing_length passed to a screened-Poisson Vector_Projection → different JIT C source generated per rank → the JIT determinism guard correctly raised a "C-source hash differs across ranks" RuntimeError. Symptom looked like a JIT bug but was a rank-local-data bug being correctly caught. Also decorates all three radii methods with @uw.collective_operation so any future use inside selective_ranks() blocks raises a clear deadlock-warning before the run hangs. The methods are documented as parallel-safe and the get_mean_radius docstring explicitly warns against falling back to self._radii.mean() with the JIT-leak explanation, so the same trap doesn't get rebuilt in user scripts. Underworld development team with AI support from Claude Code --- .../discretisation/discretisation_mesh.py | 28 +++++++++++++++++++ 1 file changed, 28 insertions(+) diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py index b2527343..800cea7d 100644 --- a/src/underworld3/discretisation/discretisation_mesh.py +++ b/src/underworld3/discretisation/discretisation_mesh.py @@ -3705,6 +3705,7 @@ def get_min_radius_old(self) -> float: return self._min_radius + @uw.collective_operation def get_min_radius(self) -> float: """ This method returns the global minimum distance from any cell centroid to a face. @@ -3721,6 +3722,7 @@ def get_min_radius(self) -> float: return all_min_radii.min() + @uw.collective_operation def get_max_radius(self) -> float: """ This method returns the global maximum distance from any cell centroid to a face. @@ -3735,6 +3737,32 @@ def get_max_radius(self) -> float: return all_max_radii.max() + @uw.collective_operation + def get_mean_radius(self) -> float: + """ + Global mean distance from any cell centroid to a face — the + characteristic background cell length scale. Parallel-safe via + MPI allreduce of the local sum and count. + + Together with :meth:`get_min_radius` / :meth:`get_max_radius` + 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). + """ + + import numpy as np + + radii = np.asarray(self._radii) + local_sum = float(radii.sum()) + local_n = int(radii.size) + if uw.mpi.size > 1: + local_sum = uw.mpi.comm.allreduce(local_sum) + local_n = uw.mpi.comm.allreduce(local_n) + return local_sum / max(local_n, 1) + # This should be deprecated in favour of using integrals def stats(self, uw_function, uw_meshVariable, basis=None): """ From b7694fac28c89b4e3afad72bad9542d0d6d71a48 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Mon, 25 May 2026 10:47:54 +1000 Subject: [PATCH 2/2] =?UTF-8?q?mesh:=20address=20copilot=20review=20?= =?UTF-8?q?=E2=80=94=20explicit=20MPI.SUM,=20clearer=20docstring,=20tests?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Three review comments from copilot-pull-request-reviewer on PR #206 addressed: * get_mean_radius now passes op=MPI.SUM explicitly to allreduce instead of relying on the default (matches the codebase pattern used elsewhere; less brittle if mpi4py wrapper conventions change) * docstring corrected — clarifies that _radii is the characteristic cell length (volume^(1/dim)) computed by DMPlexComputeGeometryFVM, not literally a centroid-to-face distance. The same correction applies to get_min/max_radius but those existing docstrings are not touched in this PR * test coverage added: - tests/test_0008_mesh_radii_accessors.py: serial unit tests on Annulus and UnstructuredSimplexBox checking min/mean/max ordering and that all three methods carry the @uw.collective_operation decorator - tests/parallel/ptest_0008_mesh_radii_accessors.py: MPI test asserting all ranks see identical min/mean/max to 1e-12 Underworld development team with AI support from Claude Code --- .../discretisation/discretisation_mesh.py | 13 +++-- .../ptest_0008_mesh_radii_accessors.py | 46 +++++++++++++++++ tests/test_0008_mesh_radii_accessors.py | 50 +++++++++++++++++++ 3 files changed, 104 insertions(+), 5 deletions(-) create mode 100644 tests/parallel/ptest_0008_mesh_radii_accessors.py create mode 100644 tests/test_0008_mesh_radii_accessors.py diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py index 800cea7d..bc83ad18 100644 --- a/src/underworld3/discretisation/discretisation_mesh.py +++ b/src/underworld3/discretisation/discretisation_mesh.py @@ -3740,9 +3740,11 @@ def get_max_radius(self) -> float: @uw.collective_operation def get_mean_radius(self) -> float: """ - Global mean distance from any cell centroid to a face — the - characteristic background cell length scale. Parallel-safe via - MPI allreduce of the local sum and count. + Global mean of the characteristic cell length scale + (``volume^(1/dim)``, i.e. the equivalent radius derived from each + cell's volume — the same quantity averaged by ``get_min_radius`` + and ``get_max_radius`` to obtain global min/max). Parallel-safe + via MPI allreduce of the local sum and count. Together with :meth:`get_min_radius` / :meth:`get_max_radius` this is the canonical "mesh length" API. Use this anywhere you @@ -3754,13 +3756,14 @@ def get_mean_radius(self) -> float: """ import numpy as np + from mpi4py import MPI radii = np.asarray(self._radii) local_sum = float(radii.sum()) local_n = int(radii.size) if uw.mpi.size > 1: - local_sum = uw.mpi.comm.allreduce(local_sum) - local_n = uw.mpi.comm.allreduce(local_n) + local_sum = uw.mpi.comm.allreduce(local_sum, op=MPI.SUM) + local_n = uw.mpi.comm.allreduce(local_n, op=MPI.SUM) return local_sum / max(local_n, 1) # This should be deprecated in favour of using integrals diff --git a/tests/parallel/ptest_0008_mesh_radii_accessors.py b/tests/parallel/ptest_0008_mesh_radii_accessors.py new file mode 100644 index 00000000..d56b40db --- /dev/null +++ b/tests/parallel/ptest_0008_mesh_radii_accessors.py @@ -0,0 +1,46 @@ +"""MPI smoke test for the parallel-safe mesh radius accessors. + +Run via mpi_runner.sh (mpirun -np N python ptest_0008_*.py). +Asserts: + - min/mean/max return positive, ordered values on every rank + - every rank sees IDENTICAL min, mean, max (parallel-consistent) + - the methods are tagged as @uw.collective_operation +""" +import underworld3 as uw + +mesh = uw.meshing.Annulus( + radiusOuter=1.0, radiusInner=0.5, cellSize=1.0 / 16, qdegree=3) + +rmin = mesh.get_min_radius() +rmean = mesh.get_mean_radius() +rmax = mesh.get_max_radius() + +assert rmin > 0.0, f"rank {uw.mpi.rank}: min must be positive, got {rmin}" +assert rmin <= rmean <= rmax, ( + f"rank {uw.mpi.rank}: expected min <= mean <= max, got " + f"{rmin} <= {rmean} <= {rmax}") + +# All ranks must agree to within float-round-off (the values come from +# global allreduce / gather, so any disagreement would indicate the +# accessor leaked rank-local state). +all_min = uw.mpi.comm.allgather(rmin) +all_mean = uw.mpi.comm.allgather(rmean) +all_max = uw.mpi.comm.allgather(rmax) + +tol = 1.0e-12 +assert max(all_min) - min(all_min) < tol, ( + f"min_radius differs across ranks: {all_min}") +assert max(all_mean) - min(all_mean) < tol, ( + f"mean_radius differs across ranks: {all_mean}") +assert max(all_max) - min(all_max) < tol, ( + f"max_radius differs across ranks: {all_max}") + +# Decorator must be present so any selective_ranks() misuse fails fast +for name in ("get_min_radius", "get_mean_radius", "get_max_radius"): + fn = getattr(mesh, name) + assert getattr(fn, "_is_collective", False), ( + f"rank {uw.mpi.rank}: {name} must be @uw.collective_operation") + +if uw.mpi.rank == 0: + print(f"OK: min={rmin:.6f} mean={rmean:.6f} max={rmax:.6f} " + f"agree across {uw.mpi.size} ranks", flush=True) diff --git a/tests/test_0008_mesh_radii_accessors.py b/tests/test_0008_mesh_radii_accessors.py new file mode 100644 index 00000000..7c79378c --- /dev/null +++ b/tests/test_0008_mesh_radii_accessors.py @@ -0,0 +1,50 @@ +"""Smoke tests for the parallel-safe mesh radius accessors. + +Verifies that :meth:`Mesh.get_min_radius`, :meth:`Mesh.get_mean_radius`, +and :meth:`Mesh.get_max_radius` return sensible values in serial. +A matching MPI smoke test lives in tests/parallel/. +""" +import pytest + +pytestmark = pytest.mark.level_1 + + +def _check_radii_ordering(mesh): + rmin = mesh.get_min_radius() + rmean = mesh.get_mean_radius() + rmax = mesh.get_max_radius() + assert rmin > 0.0, f"min radius must be positive, got {rmin}" + assert rmin <= rmean <= rmax, ( + f"expected min <= mean <= max, got " + f"min={rmin}, mean={rmean}, max={rmax}") + + +def test_mesh_radii_accessors_annulus(): + import underworld3 as uw + mesh = uw.meshing.Annulus( + radiusOuter=1.0, radiusInner=0.5, + cellSize=1.0 / 16, qdegree=3) + _check_radii_ordering(mesh) + + +def test_mesh_radii_accessors_box(): + import underworld3 as uw + mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), + cellSize=1.0 / 16) + _check_radii_ordering(mesh) + + +def test_mesh_radii_accessors_are_collective(): + """The @uw.collective_operation decorator must be present so a + selective_ranks() misuse fails fast rather than deadlocking.""" + import underworld3 as uw + mesh = uw.meshing.Annulus( + radiusOuter=1.0, radiusInner=0.5, + cellSize=1.0 / 16, qdegree=3) + for name in ("get_min_radius", "get_mean_radius", + "get_max_radius"): + fn = getattr(mesh, name) + assert getattr(fn, "_is_collective", False), ( + f"{name} must be decorated with @uw.collective_operation " + f"so misuse inside selective_ranks() is caught early.")