From 538fa8a6cf46dc901788be070812473fa4bf1223 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Fri, 14 Aug 2026 10:20:34 +1000 Subject: [PATCH 1/6] gather_data must keep one row per rank: NaN-stripping is now opt-in The gathered table is normally one row per rank, and callers index it by rank. Silently dropping a rank's NaN row shifted every later row up, so the row index no longer matched the rank that contributed it -- the mechanism behind the nearest-centroid particle mis-route on starved ranks (#399), and the trap any empty-rank reduction falls into next (#405). Stripping is still available as strip_nan=True, documented with the reason it is not the default. Underworld development team with AI support from Claude Code --- src/underworld3/utilities/_utils.py | 35 +++++++++++++++++++++-------- 1 file changed, 26 insertions(+), 9 deletions(-) diff --git a/src/underworld3/utilities/_utils.py b/src/underworld3/utilities/_utils.py index 96ae67c48..92fe3e716 100755 --- a/src/underworld3/utilities/_utils.py +++ b/src/underworld3/utilities/_utils.py @@ -193,15 +193,33 @@ def mem_footprint(): return python_process.memory_info().rss // 1000000 -def gather_data(val, bcast=False, dtype="float64"): +def gather_data(val, bcast=False, dtype="float64", strip_nan=False): """ - gather values on root (bcast=False) or all (bcast = True) processors - Parameters: - vals : Values to combine into a single array on the root or all processors - - returns: - val_global : combination of values form all processors + Gather values on root (``bcast=False``) or on all (``bcast=True``) processors. + Parameters + ---------- + val : array-like or scalar + This rank's contribution. Ranks may contribute different lengths + (including zero rows). + bcast : bool, default False + Make the combined array available on every rank rather than root only. + dtype : str, default "float64" + dtype of the combined array. + strip_nan : bool, default False + Drop NaN entries from the combined array. **Off by default and + rarely what you want**: one row per rank is the usual contract, and + silently dropping a rank's NaN row shifts every later row up so the + table index no longer equals the rank that contributed it. That + renumbering was the mechanism behind the nearest-centroid particle + mis-route on starved ranks (issues #399, #405). Pass ``True`` only + where the result is an unordered bag of values and NaN means "no + contribution". + + Returns + ------- + numpy.ndarray + The concatenated contributions from all ranks, in rank order. """ comm = uw.mpi.comm @@ -233,8 +251,7 @@ def gather_data(val, bcast=False, dtype="float64"): comm.barrier() - if uw.mpi.rank == 0: - ### remove rows with NaN + if strip_nan and uw.mpi.rank == 0: val_global = val_global[~np.isnan(val_global)] comm.barrier() From 0718d74d8a95c3a301b6b141faf887baf718d284 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Fri, 14 Aug 2026 10:20:48 +1000 Subject: [PATCH 2/6] A rank with no cells must answer every global mesh query like its peers get_min_radius / get_max_radius took min() / max() of the rank's own _radii array. On a rank owning zero cells that array is empty, so the reduction raised a rank-local ValueError while the populated ranks sat in the matching collective -- the job hung or aborted asymmetrically. This fires before the interpolation path is reached, and get_min_radius feeds estimate_dt and the penalty scaling at some fourteen call sites, so a starved rank took down ordinary time-stepping, not just point evaluation. Both accessors now allreduce the local extremum, with an empty rank contributing the identity element (+inf for a MIN, -inf for a MAX), which is the shape get_mean_radius already used. Every rank gets the same global answer, including the starved ones. The same class of defect, fixed here too: - points_in_domain reached get_max_radius (collective) only after a rank-local short-circuit. It now takes the collective first, then answers False everywhere on a rank that owns no cells -- which is the honest answer, and avoids interrogating an empty cell set. - quality() chose between its simplex and volume-only branches from rank-local evidence. A starved rank collects no triangles, so it ran three reductions where its peers ran eleven: mismatched collective counts, i.e. a hang on any simplex mesh with an empty rank. The branch is now decided collectively and each reduction is empty-safe. - physical_bounds / physical_extent reduced over this rank's nodes only, so every rank reported a different "domain size" and an empty rank raised. Both now share _global_coord_bounds, which reduces across ranks. A suspected double application of the length scale in that path is flagged with a TODO(BUG) rather than changed here. Underworld development team with AI support from Claude Code --- .../discretisation/discretisation_mesh.py | 191 ++++++++++++++---- 1 file changed, 151 insertions(+), 40 deletions(-) diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py index f141d8593..8be2f96c3 100644 --- a/src/underworld3/discretisation/discretisation_mesh.py +++ b/src/underworld3/discretisation/discretisation_mesh.py @@ -1504,6 +1504,21 @@ def _reduce(val, op): val, op=getattr(_MPI, op)) return val + # A rank owning zero cells contributes the identity element of each + # reduction rather than raising on an empty array (issue #405). + def _reduce_min(arr): + return _reduce(float(arr.min()) if arr.size else float("inf"), + "MIN") + + def _reduce_max(arr): + return _reduce(float(arr.max()) if arr.size else float("-inf"), + "MAX") + + def _local_percentile(arr, pct): + # Rank-local estimate (see the docstring); a rank with no cells + # has no local distribution to take a percentile of. + return float(np.percentile(arr, pct)) if arr.size else float("nan") + tri_vertex_lists = [] is_simplex2d = cdim == 2 if is_simplex2d: @@ -1516,17 +1531,32 @@ def _reduce(val, op): break tri_vertex_lists.append(cell_vertices) - if not is_simplex2d or not tri_vertex_lists: + # Choose the branch COLLECTIVELY. A rank owning zero cells collects no + # triangles and would otherwise take the volume-only branch (three + # reductions) while its populated peers took the simplex branch + # (eleven) — mismatched collective counts, i.e. a hang (issue #405). + # A starved rank abstains from the vote instead. + has_cells = cEnd > cStart + n_simplex_ranks = _reduce( + int(has_cells and is_simplex2d and bool(tri_vertex_lists)), "SUM") + n_populated_ranks = _reduce(int(has_cells), "SUM") + is_simplex2d = (n_populated_ranks > 0 + and n_simplex_ranks == n_populated_ranks) + + if not is_simplex2d: try: volume = np.abs(np.array( [dm.computeCellGeometryFVM(cell_id)[0] for cell_id in range(cStart, cEnd)])) except Exception: volume = np.array([1.0]) - if not volume.size: + # A rank that owns cells but cannot compute their geometry keeps + # the unit-volume placeholder; a rank owning NO cells contributes + # nothing at all, so the global cell count stays honest. + if not volume.size and has_cells: volume = np.array([1.0]) n_cells = _reduce(int(volume.size), "SUM") - vol_min = _reduce(float(volume.min()), "MIN") + vol_min = _reduce_min(volume) vol_sum = _reduce(float(volume.sum()), "SUM") metrics = dict( n_cells=n_cells, element="non-2D-simplex", @@ -1538,7 +1568,9 @@ def _reduce(val, op): metrics["per_cell"] = dict(volume=volume) return metrics - tri = np.asarray(tri_vertex_lists, dtype=np.int64) + # reshape(-1, 3) keeps the (0, 3) shape on a rank with no cells, where + # np.asarray([]) would be 1-D and the column indexing below would fail. + tri = np.asarray(tri_vertex_lists, dtype=np.int64).reshape(-1, 3) v0, v1, v2 = (vertex_coords[tri[:, 0]], vertex_coords[tri[:, 1]], vertex_coords[tri[:, 2]]) @@ -1566,7 +1598,9 @@ def _angle_deg(opposite, side1, side2): _angle_deg(edge_c, edge_a, edge_b)]) longest_edge = np.maximum.reduce([edge_a, edge_b, edge_c]) aspect = longest_edge * longest_edge / (2.0 * area) - rel_area = area / area.mean() + # rel_area only ever masks this rank's own cells, so the rank-local + # mean is the right scale — and an empty rank has no cells to mask. + rel_area = area / area.mean() if area.size else area # Neighbour size-jump: map each (undirected) edge to the triangles # sharing it; interior edges (exactly two triangles) contribute the @@ -1585,22 +1619,22 @@ def _angle_deg(opposite, side1, side2): area_sum = _reduce(float(area.sum()), "SUM") metrics = dict( n_cells=n_cells, element="2D-simplex", - q_min=_reduce(float(shape_q.min()), "MIN"), + q_min=_reduce_min(shape_q), q_mean=q_sum / max(n_cells, 1), - q_p01=float(np.percentile(shape_q, 1)), - q_p05=float(np.percentile(shape_q, 5)), + q_p01=_local_percentile(shape_q, 1), + q_p05=_local_percentile(shape_q, 5), n_q_lt_0p3=_reduce(int((shape_q < 0.3).sum()), "SUM"), n_q_lt_0p2=_reduce(int((shape_q < 0.2).sum()), "SUM"), - angle_max_deg=_reduce(float(largest_angle.max()), "MAX"), + angle_max_deg=_reduce_max(largest_angle), n_angle_gt_150=_reduce(int((largest_angle > 150).sum()), "SUM"), n_angle_gt_165=_reduce(int((largest_angle > 165).sum()), "SUM"), - aspect_max=_reduce(float(aspect.max()), "MAX"), - aspect_p99=float(np.percentile(aspect, 99)), + aspect_max=_reduce_max(aspect), + aspect_p99=_local_percentile(aspect, 99), sizejump_max=float(size_jump.max()), - sizejump_p99=float(np.percentile(size_jump, 99)), + sizejump_p99=_local_percentile(size_jump, 99), n_big_thin=_reduce( int(((rel_area > 2.0) & (aspect > 4.0)).sum()), "SUM"), - vol_min_over_mean=(_reduce(float(area.min()), "MIN") + vol_min_over_mean=(_reduce_min(area) / (area_sum / max(n_cells, 1)))) if per_cell: metrics["per_cell"] = dict( @@ -4434,6 +4468,9 @@ def physical_bounds(self): Returns the mesh bounding box scaled to physical units using the model's length scale. + COLLECTIVE: the bounding box spans the whole mesh, so it is reduced + across ranks (see :meth:`_global_coord_bounds`). + Returns ------- tuple of UWQuantity or None @@ -4447,10 +4484,7 @@ def physical_bounds(self): if not hasattr(self, "_model") or self._model is None: return None - import numpy as np - - min_coords = np.min(self.points, axis=0) - max_coords = np.max(self.points, axis=0) + min_coords, max_coords = self._global_coord_bounds() return ( self._model.scale_to_physical(min_coords, dimension="length"), @@ -4464,6 +4498,9 @@ def physical_extent(self): Returns the mesh size (max - min) in each dimension scaled to physical units. + COLLECTIVE: the extent spans the whole mesh, so it is reduced across + ranks (see :meth:`_global_coord_bounds`). + Returns ------- UWQuantity or None @@ -4477,13 +4514,59 @@ def physical_extent(self): if not hasattr(self, "_model") or self._model is None: return None + min_coords, max_coords = self._global_coord_bounds() + + return self._model.scale_to_physical( + max_coords - min_coords, dimension="length") + + def _global_coord_bounds(self): + """Bounding box of the mesh nodes, ``(min_coords, max_coords)``. + + COLLECTIVE. Each rank holds only its own subdomain's nodes, so a + rank-local ``min``/``max`` describes the partition rather than the + mesh — every rank would report a different "domain size". The + reduction makes the answer global and identical everywhere, and lets + a rank owning no cells (hence no nodes) contribute the identity + elements instead of raising on an empty array (issue #405). + + Reads the same node coordinates as the deprecated ``mesh.points`` + (without its warning or unit wrapping), so the physical-bounds / + physical-extent answers are unchanged apart from being global. + + .. TODO(BUG): ``mesh.points`` already multiplies by + ``CoordinateSystem._length_scale`` when the coordinate system is + scaled, and both callers then pass the result through + ``model.scale_to_physical(..., dimension="length")`` — a second + application of the same factor. Reproduced here deliberately so + this fix stays behaviour-neutral; the double scaling is a separate + question for the units owner. + """ import numpy as np + from mpi4py import MPI - min_coords = np.min(self.points, axis=0) - max_coords = np.max(self.points, axis=0) - extent = max_coords - min_coords + coords = np.asarray(self._coords, dtype=np.float64).reshape( + -1, self.cdim) + if getattr(self.CoordinateSystem, "_scaled", False): + coords = coords * self.CoordinateSystem._length_scale + coords = np.ascontiguousarray(coords) + + if coords.shape[0] > 0: + local_min = np.ascontiguousarray(coords.min(axis=0)) + local_max = np.ascontiguousarray(coords.max(axis=0)) + else: + local_min = np.full(self.cdim, np.inf) + local_max = np.full(self.cdim, -np.inf) + + if uw.mpi.size > 1: + # Buffer (uppercase) Allreduce: the pickling `allreduce` applies + # MPI.MIN through Python's `min()`, which is ambiguous for arrays. + global_min = np.empty_like(local_min) + global_max = np.empty_like(local_max) + uw.mpi.comm.Allreduce(local_min, global_min, op=MPI.MIN) + uw.mpi.comm.Allreduce(local_max, global_max, op=MPI.MAX) + local_min, local_max = global_min, global_max - return self._model.scale_to_physical(extent, dimension="length") + return local_min, local_max @timing.routine_timer_decorator def write_timestep( @@ -5876,13 +5959,24 @@ def points_in_domain(self, points, strict_validation=True): # and handles all the complexity of extracting values from unit-aware coordinates model_points = _convert_coords_to_si(points) - self._mark_local_boundary_faces_inside_and_out() - + # get_max_radius() is COLLECTIVE, so it must be reached by every rank + # before any rank takes a short-circuit below — otherwise the starved + # ranks skip the reduction their peers are sitting in (issue #405). max_radius = self.get_max_radius() + self._mark_local_boundary_faces_inside_and_out() + if model_points.shape[0] == 0: return numpy.array([], dtype=bool) + # A rank owning no cells contains no points, so the honest answer is + # False everywhere. Its local boundary skeleton is empty too, and the + # closest-local-cell test below would otherwise have to interrogate a + # cell set that does not exist. + cStart, cEnd = self.dm.getHeightStratum(0) + if cEnd == cStart: + return numpy.zeros(model_points.shape[0], dtype=bool) + # Cd-1 surface mesh: no boundary-face control points exist # (see _mark_local_boundary_faces_inside_and_out). Per the # surface-mesh contract, query points are assumed to lie on @@ -6453,13 +6547,13 @@ def _get_domain_centroids(self): from underworld3.utilities import gather_data # 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.) + # array is NaN. gather_data no longer strips NaN rows (issue #405 + # made that opt-in), but a NaN row would still poison the kd-tree + # this table feeds, and _route_by_nearest_centroid would mis-route + # particles (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: @@ -6499,34 +6593,51 @@ def get_min_radius_old(self) -> float: @uw.collective_operation def get_min_radius(self) -> float: """ - This method returns the global minimum distance from any cell centroid to a face. - It wraps to the PETSc `DMPlexGetMinRadius` routine. The petsc4py equivalent always - returns zero. + Global minimum of the characteristic cell length scale — the smallest + cell anywhere in the mesh, not just on this rank. Parallel-safe via + MPI allreduce of the local minimum. + + A rank owning zero cells contributes the identity element of the + reduction (:math:`+\\infty`) and therefore returns the same global + value as its populated peers. Taking ``min()`` of that rank's empty + ``_radii`` array instead would raise on the starved rank alone, while + its peers waited in the reduction — the rank-asymmetric raise that + deadlocks the job (issue #405). """ ## Note: The petsc4py version of DMPlexComputeGeometryFVM does not compute all cells and ## does not obtain the minimum radius for the mesh. import numpy as np + from mpi4py import MPI - all_min_radii = uw.utilities.gather_data(np.array((self._radii.min(),)), bcast=True) - - return all_min_radii.min() + radii = np.asarray(self._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) + return local_min @uw.collective_operation def get_max_radius(self) -> float: """ - This method returns the global maximum distance from any cell centroid to a face. + Global maximum of the characteristic cell length scale — the largest + cell anywhere in the mesh. Parallel-safe via MPI allreduce of the + local maximum; a rank owning zero cells contributes the identity + element (:math:`-\\infty`) and still returns the global value. + See :meth:`get_min_radius` for why the guard matters. """ ## Note: The petsc4py version of DMPlexComputeGeometryFVM does not compute all cells and ## does not obtain the minimum radius for the mesh. import numpy as np + from mpi4py import MPI - all_max_radii = uw.utilities.gather_data(np.array((self._radii.max(),)), bcast=True) - - return all_max_radii.max() + radii = np.asarray(self._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) + return local_max @uw.collective_operation def get_mean_radius(self) -> float: From 6cdf62790fa94af5c8c10cb5ea95859ccbaec198 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Fri, 14 Aug 2026 10:20:59 +1000 Subject: [PATCH 3/6] Rank-local interpolation must not index an empty source cloud A rank owning no cells owns no DOFs, so there is nothing to interpolate from. MeshVariable.rbf_interpolate built its kd-tree over an empty point cloud and the monotone limiter's kNN bound gathered neighbour values from an empty array (#405). The SwarmVariable rbf path already guarded this way. Both are purely rank-local, so returning early takes no collective with it: rbf_interpolate returns correctly-shaped zeros, and the limiter returns the (empty) value unchanged -- exact, not a fallback. monotone='pick', the only collective mode, is already refused under MPI. Underworld development team with AI support from Claude Code --- .../discretisation/discretisation_mesh_variables.py | 9 +++++++++ src/underworld3/function/functions_unit_system.py | 8 ++++++++ 2 files changed, 17 insertions(+) diff --git a/src/underworld3/discretisation/discretisation_mesh_variables.py b/src/underworld3/discretisation/discretisation_mesh_variables.py index 8c95671f0..e4e563c9a 100644 --- a/src/underworld3/discretisation/discretisation_mesh_variables.py +++ b/src/underworld3/discretisation/discretisation_mesh_variables.py @@ -935,6 +935,15 @@ def rbf_interpolate(self, new_coords, nnn=None, p=1, verbose=False, D = self.data.copy() + # A rank owning no cells owns no DOFs either, so there is nothing to + # interpolate FROM: the kd-tree below would be built over an empty + # point cloud and the stencil gather would index an empty array + # (issue #405). Return the correctly-shaped zeros — this is a purely + # rank-local path, so returning early takes no collective with it. + # (The equivalent SwarmVariable path guards the same way.) + if D.shape[0] == 0: + return np.zeros((np.asarray(new_coords).shape[0], D.shape[1])) + if verbose and uw.mpi.rank == 0: print("Building K-D tree", flush=True) diff --git a/src/underworld3/function/functions_unit_system.py b/src/underworld3/function/functions_unit_system.py index d9289c6e0..ecfea55bf 100644 --- a/src/underworld3/function/functions_unit_system.py +++ b/src/underworld3/function/functions_unit_system.py @@ -698,6 +698,14 @@ def _apply_monotone_limit( # latent mismatch (scaling is inactive in the validated baseline so it # never bites). Do not "fix" without re-validating the trajectory. nnn = mesh.dim + 1 + # A rank owning no cells owns no source DOFs: there is no neighbourhood + # to bound against, and the stencil gather would index an empty array + # (issue #405). Such a rank also has no interior evaluation points, so + # `value` is empty and returning it unchanged is exact, not a fallback. + # Purely rank-local: "pick" (the only collective mode) is refused above + # under MPI, so this early return cannot skip a collective. + if psi_coords_nd.shape[0] == 0 or np.asarray(coords_nd).shape[0] == 0: + return value kdt = uw.kdtree.KDTree(np.ascontiguousarray(psi_coords_nd)) _, idxs = kdt.query( np.ascontiguousarray(coords_nd), k=nnn, sqr_dists=False) From 67a332ba200154ba48e363e17c627d627c359b95 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Fri, 14 Aug 2026 10:20:59 +1000 Subject: [PATCH 4/6] estimate_dt must not raise on a rank that owns no cells Three reductions inside the timestep estimators read this rank's own cells and reduced before the allreduce, so a starved rank raised while its peers waited in the collective (#405): - _global_max_diffusivity took .max() of the diffusivity sampled at this rank's cell centroids. The max is now taken once, in the reduction, where an empty rank contributes the identity element. - SNES_NavierStokes.estimate_dt took .max() of the centroid velocity magnitudes. - SNES_Stokes_SaddlePt.estimate_dt took .max() of the local velocity DOFs. Underworld development team with AI support from Claude Code --- .../cython/petsc_generic_snes_solvers.pyx | 5 ++++- src/underworld3/systems/solvers.py | 17 +++++++++++++---- 2 files changed, 17 insertions(+), 5 deletions(-) diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index 92351010c..d20ad05f2 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -9600,7 +9600,10 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): if self.mesh.dim == 3: magvel_squared += vel[:, 2] ** 2 - max_magvel = math.sqrt(magvel_squared.max()) + # A rank owning no cells owns no velocity DOFs; it contributes the + # identity element of the MAX rather than raising on the empty array + # while its peers wait in the allreduce (issue #405). + max_magvel = math.sqrt(magvel_squared.max()) if magvel_squared.size else 0.0 from mpi4py import MPI diff --git a/src/underworld3/systems/solvers.py b/src/underworld3/systems/solvers.py index fe997d8ff..8b0951760 100644 --- a/src/underworld3/systems/solvers.py +++ b/src/underworld3/systems/solvers.py @@ -223,10 +223,12 @@ def _global_max_diffusivity(constitutive_K, mesh): diffusivity = uw.function.evaluate( K_sym, np.zeros((1, mesh.dim))) else: - # Spatially varying: sample at cell centroids, take the max. + # Spatially varying: sample at cell centroids. The maximum is + # taken once, in the reduction below — a rank owning no cells + # samples nothing, and `.max()` of that empty array would raise + # here while the peers waited in the allreduce (issue #405). diffusivity = uw.function.evaluate( sympy.sympify(K_sym), mesh._centroids, mesh.N) - diffusivity = diffusivity.max() else: diffusivity = K @@ -239,7 +241,10 @@ def _global_max_diffusivity(constitutive_K, mesh): # Plain UWQuantity without units context - use magnitude diffusivity = diffusivity.magnitude - local_max = float(np.asarray(diffusivity).max()) + # A rank with no local samples contributes the identity element of the + # MAX and still receives the correct global value. + local_values = np.asarray(diffusivity) + local_max = float(local_values.max()) if local_values.size else float("-inf") return uw.mpi.comm.allreduce(local_max, op=MPI.MAX) @@ -5056,7 +5061,11 @@ def estimate_dt(self): # the maximum |component| rather than the maximum vector magnitude # (the other estimate_dt implementations squeeze first). Preserved # as-is (Wave D is behaviour-neutral); revisit with a numerical check. - max_magvel = np.linalg.norm(vel, axis=1).max() + # A rank owning no cells has no centroid samples; it contributes the + # identity element of the MAX rather than raising on the empty array + # while its peers wait in the allreduce (issue #405). + magvel = np.linalg.norm(vel, axis=1) + max_magvel = float(magvel.max()) if magvel.size else 0.0 max_magvel_glob = comm.allreduce(max_magvel, op=MPI.MAX) ## get radius From 78dd1b33b7d998d911e2d86ce75606c5ad9f97df Mon Sep 17 00:00:00 2001 From: lmoresi Date: Fri, 14 Aug 2026 10:21:11 +1000 Subject: [PATCH 5/6] Regression test: global reductions on a genuinely zero-cell rank A 1x2 quad box has two cells, so PETSc must leave at least one rank empty at both np=2 and np=4. The tests assert what the fix is for: min/max/mean radius, quality(), estimate_dt and points_in_domain give every rank the same global answer, equal to the serial one, with no hang. Two of them earn their keep explicitly. The premise test fails if the partition is not actually starved, and the negative control shows the cross-rank agreement assertion has teeth: the rank-local minima on this partition do NOT agree, so a rank-local answer would be caught. Fail-before on the same fixture (build without the fix): at np=2 one rank raised "zero-size array to reduction operation minimum which has no identity" from get_min_radius while the other passed, then the run HUNG; at np=4 three ranks raised and the run HUNG. Both needed the mpirun timeout to end. Also records a measured finding in mesh_metric_mismatch: the empty-rank limit on field-valued metrics does NOT lift with this fix, as #405 item 2 expected. Measured at np=4 on a starved mesh, the remaining blocker is below UW3 -- the DMPlex sub-DM clone the mesh-variable path builds fails with MPI_ERR_BUFFER when a rank has no cells. That is #314's territory. Underworld development team with AI support from Claude Code --- src/underworld3/meshing/smoothing/metrics.py | 15 +- .../test_0774_empty_rank_reductions_mpi.py | 182 ++++++++++++++++++ 2 files changed, 192 insertions(+), 5 deletions(-) create mode 100644 tests/parallel/test_0774_empty_rank_reductions_mpi.py diff --git a/src/underworld3/meshing/smoothing/metrics.py b/src/underworld3/meshing/smoothing/metrics.py index 2aae84e9f..976840439 100644 --- a/src/underworld3/meshing/smoothing/metrics.py +++ b/src/underworld3/meshing/smoothing/metrics.py @@ -171,11 +171,16 @@ def mesh_metric_mismatch(mesh, metric, resolution_ratio=None): # # 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. + # rank then deadlocks the populated ranks inside evaluate. Analytic + # (pure-sympy) metrics are fine: their evaluation is rank-local. + # + # Issue #405 made the reduction layer under evaluate empty-rank safe + # (radii, points_in_domain), which was expected to lift this limit on + # its own. Measured at np=4 on a starved mesh: it does not. The + # remaining blocker is below UW3 — the DMPlex sub-DM clone that the + # mesh-variable path builds fails with MPI_ERR_BUFFER when a rank has + # no cells. That is issue #314's territory; revisit this branch when + # #314 closes. 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/parallel/test_0774_empty_rank_reductions_mpi.py b/tests/parallel/test_0774_empty_rank_reductions_mpi.py new file mode 100644 index 000000000..e763ddb8d --- /dev/null +++ b/tests/parallel/test_0774_empty_rank_reductions_mpi.py @@ -0,0 +1,182 @@ +"""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 +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 +element (+inf for a MIN, -inf for a MAX, 0 for a SUM) — so a rank with no +cells returns the SAME global answer as everyone else. + +Run under MPI, e.g.:: + + mpirun -np 2 python -m pytest --with-mpi \ + tests/parallel/test_0774_empty_rank_reductions_mpi.py + mpirun -np 4 python -m pytest --with-mpi \ + tests/parallel/test_0774_empty_rank_reductions_mpi.py + +The fixture is a 1x2 quad box — TWO cells for two-or-more ranks — so PETSc +must leave at least one rank empty at both np=2 and np=4. Every test is +timeout-guarded: the pre-fix failure mode is a hang, not an exception. +""" + +import math + +import numpy as np +import pytest + +import underworld3 as uw + +pytestmark = [ + pytest.mark.mpi(min_size=2), + pytest.mark.timeout(120), + pytest.mark.level_1, + pytest.mark.tier_a, +] + +# Cell is 1.0 wide x 0.5 tall; the characteristic length UW3 reports is the +# centroid-to-corner half-diagonal. Analytic, so this is a real oracle rather +# than a recorded number. +SERIAL_RADIUS = math.sqrt(0.5 ** 2 + 0.25 ** 2) + + +def _starved_box(): + """A 2-cell mesh: at np >= 2 some rank is guaranteed to own no cells.""" + return uw.meshing.StructuredQuadBox(elementRes=(1, 2)) + + +def _local_cell_count(mesh): + cStart, cEnd = mesh.dm.getHeightStratum(0) + return cEnd - cStart + + +def test_premise_some_rank_owns_no_cells(): + """Test premise: without a genuinely empty rank nothing here is tested.""" + mesh = _starved_box() + counts = uw.mpi.comm.allgather(_local_cell_count(mesh)) + + assert sum(counts) >= 2, f"fixture lost cells: {counts}" + assert min(counts) == 0, ( + f"no rank was starved at np={uw.mpi.size}: cells per rank {counts} — " + "this test cannot detect the #405 defect on this partition") + assert max(counts) > 0, f"every rank starved: {counts}" + + +def test_radius_accessors_are_global_on_a_zero_cell_rank(): + """min/max/mean radius: same value on every rank, equal to the serial one.""" + mesh = _starved_box() + + r_min = mesh.get_min_radius() + r_max = mesh.get_max_radius() + r_mean = mesh.get_mean_radius() + + for name, value in (("min", r_min), ("max", r_max), ("mean", r_mean)): + gathered = uw.mpi.comm.allgather(value) + assert max(gathered) - min(gathered) < 1.0e-12, ( + f"get_{name}_radius disagrees across ranks: {gathered} — a " + "starved rank must not return a rank-local answer") + assert np.isclose(value, SERIAL_RADIUS, rtol=1.0e-9), ( + f"get_{name}_radius = {value}, expected the serial " + f"{SERIAL_RADIUS}") + + +def test_negative_control_rank_local_minimum_would_be_caught(): + """Prove the cross-rank agreement assertion above has teeth. + + If ``get_min_radius`` returned this rank's own minimum (the pre-fix + intent, minus the raise), the ranks would NOT agree — so the assertion + in the previous test is not true by construction. + """ + mesh = _starved_box() + + radii = np.asarray(mesh._radii).reshape(-1) + rank_local_min = float(radii.min()) if radii.size else float("inf") + gathered = uw.mpi.comm.allgather(rank_local_min) + + assert max(gathered) - min(gathered) > 1.0e-12, ( + f"rank-local minima {gathered} happen to agree, so this partition " + "cannot distinguish a global answer from a rank-local one") + assert math.isinf(max(gathered)), ( + "expected the starved rank's rank-local minimum to be the identity " + f"element, got {gathered}") + + +def test_points_in_domain_answers_false_on_a_zero_cell_rank(): + """A rank with no cells contains no points — and must not raise.""" + mesh = _starved_box() + + # Interior of the lower cell (avoids the internal face at y = 0.5). + query = np.array([[0.5, 0.25]]) + in_or_not = mesh.points_in_domain(query) + + n_local = _local_cell_count(mesh) + if n_local == 0: + assert not in_or_not.any(), ( + f"rank {uw.mpi.rank} owns no cells but claimed {query}") + + claims = uw.mpi.comm.allreduce(int(in_or_not.any()), op=uw.MPI.SUM) + assert claims >= 1, "no rank claimed an interior point" + + +def test_quality_diagnostic_agrees_across_ranks(): + """quality() reduces globally; the branch choice must be rank-symmetric.""" + mesh = _starved_box() + + n_cells = mesh.quality()["n_cells"] + gathered = uw.mpi.comm.allgather(n_cells) + + assert len(set(gathered)) == 1, ( + f"quality()['n_cells'] disagrees across ranks: {gathered}") + assert n_cells == 2, f"expected the fixture's 2 cells, got {n_cells}" + + +def test_estimate_dt_agrees_across_ranks(): + """estimate_dt is the reason #405 was raised in priority: it feeds every + time-stepping loop through get_min_radius and the diffusivity reduction.""" + mesh = _starved_box() + T = uw.discretisation.MeshVariable("T405", mesh, 1, degree=1) + + solver = uw.systems.Diffusion(mesh, u_Field=T) + solver.constitutive_model = uw.constitutive_models.DiffusionModel + solver.constitutive_model.Parameters.diffusivity = 1.0 + + dt_const = float(solver.estimate_dt()) + gathered = uw.mpi.comm.allgather(dt_const) + assert max(gathered) - min(gathered) < 1.0e-12, ( + f"estimate_dt disagrees across ranks: {gathered}") + # Diffusive CFL with unit diffusivity is exactly h_min**2. + assert np.isclose(dt_const, mesh.get_min_radius() ** 2, rtol=1.0e-12) + + # A spatially varying diffusivity is sampled at THIS rank's cell + # centroids — a starved rank samples nothing, which is where the + # unguarded `.max()` used to raise. + x, _y = mesh.X + solver.constitutive_model.Parameters.diffusivity = 1.0 + x + dt_varying = float(solver.estimate_dt()) + gathered = uw.mpi.comm.allgather(dt_varying) + assert max(gathered) - min(gathered) < 1.0e-12, ( + f"estimate_dt (varying K) disagrees across ranks: {gathered}") + assert dt_varying < dt_const, ( + "a larger diffusivity must shorten the diffusive timestep") + + +def test_gather_data_keeps_one_row_per_rank(): + """#405 item 3: NaN rows must survive, so row index still equals rank.""" + contribution = np.array( + [float("nan") if uw.mpi.rank % 2 else float(uw.mpi.rank)]) + + table = uw.utilities.gather_data(contribution, bcast=True) + + assert table.shape[0] == uw.mpi.size, ( + f"gather_data returned {table.shape[0]} rows for {uw.mpi.size} " + "ranks — a dropped row renumbers every rank after it") + for r in range(uw.mpi.size): + if r % 2: + assert np.isnan(table[r]), f"row {r} should be this rank's NaN" + else: + assert table[r] == r, f"row {r} came from the wrong rank" + + # The old behaviour is still available, explicitly. + stripped = uw.utilities.gather_data(contribution, bcast=True, + strip_nan=True) + assert stripped.shape[0] == (uw.mpi.size + 1) // 2 From 6e97b0766c633c1cfede9205c44efbaca8a83a78 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Fri, 14 Aug 2026 11:34:00 +1000 Subject: [PATCH 6/6] estimate_dt must return a scalar, not a zero-dimensional array CI regression from this branch: tests/test_1005_TransientDarcyCartesian.py failed with "TypeError: Could not convert object to sequence" out of SNES_TransientDarcy.solve. Diagnosis. Every estimate_dt funnels its nondimensional result through np.squeeze. That collapses a numpy scalar to a numpy scalar, but PROMOTES a plain Python float to a 0-d ndarray. solve(timestep=dt) then evaluates `timestep != self.delta_t` against a UWexpression, and ndarray-vs-sympy comparison tries to convert the sympy object to a sequence and raises rather than deferring to sympy's own comparison. The contract only ever held by accident: get_min_radius returned a numpy scalar, so the squeeze was a no-op. Making it return a plain float (the empty-rank work in this branch, where the reduction is now an allreduce of a Python float) turned the squeeze into a promotion and broke the solve. Fix at the contract boundary rather than by restoring the accident: the shared _apply_unit_aware_scaling helper already documents its return as "float or UWQuantity", so it now collapses a 0-d array to a scalar via a small _as_scalar helper. Unit-bearing values are passed through untouched -- .item() would discard their units. The two estimate_dt implementations that return directly instead of via that helper (SNES_Stokes, SNES_AdvectionDiffusion) use the same helper, so all five agree. This makes the contract independent of whichever numeric type a caller passes in. tests/test_1007 pins it: _as_scalar collapses 0-d arrays and only those, and TransientDarcy.estimate_dt returns something solve() actually accepts. Measured on the broken build, estimate_dt returned array(0.01586914) with ndim 0; it now returns the identical value as a float. Underworld development team with AI support from Claude Code --- src/underworld3/systems/solvers.py | 35 ++++++++-- .../test_1007_estimate_dt_scalar_contract.py | 68 +++++++++++++++++++ 2 files changed, 99 insertions(+), 4 deletions(-) create mode 100644 tests/test_1007_estimate_dt_scalar_contract.py diff --git a/src/underworld3/systems/solvers.py b/src/underworld3/systems/solvers.py index 8b0951760..77e1aecdc 100644 --- a/src/underworld3/systems/solvers.py +++ b/src/underworld3/systems/solvers.py @@ -85,6 +85,18 @@ def expression(*args, **kwargs): return public_expression(*args, _unique_name_generation=True, **kwargs) +def _as_scalar(value): + """Collapse a zero-dimensional array to a plain scalar, leave the rest. + + Unit-bearing values (UWQuantity / UnitAwareArray / pint) are returned + untouched: ``.item()`` would discard their units. + """ + if isinstance(value, np.ndarray) and value.ndim == 0 and not hasattr( + value, "units"): + return value.item() + return value + + def _apply_unit_aware_scaling(dt_nondimensional, field, mesh): """ Convert a nondimensional timestep estimate to physical time units. @@ -108,8 +120,19 @@ def _apply_unit_aware_scaling(dt_nondimensional, field, mesh): ------- float or UWQuantity Timestep with physical time units if a time scale is configured, - otherwise the nondimensional input. + otherwise the nondimensional input — always a SCALAR, never a + zero-dimensional array. The callers reach here through + ``np.squeeze(dt)``, which collapses a numpy scalar to a numpy scalar + but promotes a plain Python float to a 0-d ``ndarray``. A 0-d array + is neither of the documented return types, and it poisons the + arithmetic downstream: ``solve(timestep=...)`` compares the estimate + against ``self.delta_t`` (a UWexpression), and ndarray-vs-sympy + comparison raises ``TypeError: Could not convert object to sequence`` + instead of deferring to sympy. Collapsing it here keeps the contract + independent of whichever numeric type a caller happened to pass in. """ + dt_nondimensional = _as_scalar(dt_nondimensional) + try: from ..function.quantities import UWQuantity @@ -2153,7 +2176,9 @@ def estimate_dt(self): try: return uw.dimensionalise(np.squeeze(min_dt_glob), {'[time]': 1}) except Exception: - return np.squeeze(min_dt_glob) + # _as_scalar: np.squeeze promotes a Python float to a 0-d array + # (see _apply_unit_aware_scaling for what that breaks). + return _as_scalar(np.squeeze(min_dt_glob)) class SNES_VE_Stokes(SNES_Stokes): @@ -4207,8 +4232,10 @@ def _reduce_dt(per_elem): try: return uw.dimensionalise(np.squeeze(dt_estimate), {'[time]': 1}) except Exception: - # Fallback: return plain nondimensional number - return np.squeeze(dt_estimate) + # Fallback: return plain nondimensional number. _as_scalar because + # np.squeeze promotes a Python float to a 0-d array, which is not + # a number any caller expects (see _apply_unit_aware_scaling). + return _as_scalar(np.squeeze(dt_estimate)) @timing.routine_timer_decorator def solve( diff --git a/tests/test_1007_estimate_dt_scalar_contract.py b/tests/test_1007_estimate_dt_scalar_contract.py new file mode 100644 index 000000000..84063847b --- /dev/null +++ b/tests/test_1007_estimate_dt_scalar_contract.py @@ -0,0 +1,68 @@ +"""estimate_dt must return a SCALAR, not a zero-dimensional array. + +Every estimate_dt implementation funnels its nondimensional result through +``np.squeeze``. That collapses a numpy scalar to a numpy scalar but promotes a +plain Python float to a 0-d ``ndarray`` -- which is neither of the documented +return types, and which poisons the arithmetic downstream: ``solve(timestep=dt)`` +compares the estimate against ``self.delta_t`` (a UWexpression), and +ndarray-vs-sympy comparison raises ``TypeError: Could not convert object to +sequence`` instead of deferring to sympy. + +The contract used to hold only by accident -- ``get_min_radius`` happened to +return a numpy scalar, so the squeeze was a no-op. When it started returning a +plain float (the #405 empty-rank work), the transient Darcy solve broke. These +tests pin the contract itself so it cannot depend on a caller's numeric type +again. +""" + +import numpy as np +import pytest +import sympy as sp + +import underworld3 as uw +from underworld3.systems.solvers import _as_scalar + +pytestmark = [pytest.mark.level_1, pytest.mark.tier_a] + + +def test_as_scalar_collapses_zero_dim_arrays_only(): + """0-d arrays become scalars; everything else is passed through.""" + collapsed = _as_scalar(np.squeeze(0.5)) + assert not isinstance(collapsed, np.ndarray), ( + "a 0-d array must not survive as an array") + assert collapsed == 0.5 + + # Real arrays and plain scalars are untouched. + arr = np.array([1.0, 2.0]) + assert _as_scalar(arr) is arr + assert _as_scalar(0.5) == 0.5 + assert _as_scalar(np.float64(0.5)) == 0.5 + + +def test_estimate_dt_returns_a_usable_scalar(): + """The value estimate_dt returns must be accepted by solve(timestep=...). + + Regression: a 0-d ndarray reached ``timestep != self.delta_t`` and raised + TypeError against the UWexpression, so the solver could not be stepped with + its own timestep estimate. + """ + mesh = uw.meshing.StructuredQuadBox(elementRes=(4, 8)) + h = uw.discretisation.MeshVariable("h1007", mesh, 1, degree=2) + v = uw.discretisation.MeshVariable("v1007", mesh, mesh.dim, degree=1) + + darcy = uw.systems.TransientDarcy(mesh, h, v, order=1, theta=0.5) + darcy.constitutive_model = uw.constitutive_models.DarcyFlowModel + darcy.constitutive_model.Parameters.permeability = 1.0 + darcy.constitutive_model.Parameters.s = sp.Matrix([0, 0]).T + darcy.storage = 1.0 + darcy.f = 0.0 + + dt = darcy.estimate_dt() + + assert not (isinstance(dt, np.ndarray) and dt.ndim == 0), ( + f"estimate_dt returned a 0-d array ({dt!r}); callers compare this " + "against a UWexpression and numpy raises rather than deferring") + assert np.isfinite(float(dt)) and float(dt) > 0.0 + + # The comparison that actually broke, exercised directly. + assert (dt != darcy.delta_t) is not NotImplemented