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/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: 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) 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/src/underworld3/systems/solvers.py b/src/underworld3/systems/solvers.py index fe997d8ff..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 @@ -223,10 +246,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 +264,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) @@ -2148,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): @@ -4202,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( @@ -5056,7 +5088,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 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() 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 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