A rank with no cells must answer every global mesh query like its peers (#405)

[lmoresi]

Empty-rank support, layer 2: a rank with no cells gets the same global answer as its peers

Fixes #405.

The defect

A rank that owns zero cells raised a rank-local ValueError from an unguarded
local reduction — self._radii.min() on an empty array — while its populated peers
sat in the matching collective. The job then hung, or aborted asymmetrically.

This is distinct from #314's collective-count divergence and fires earlier, before
the interpolation path is reached. It matters because get_min_radius feeds
estimate_dt and the penalty scaling at ~14 call sites and sits under Swarm.migrate
as well as under evaluate: a starved rank took down ordinary time-stepping, not just
point evaluation.

The correctness argument

Every quantity here is global — the smallest cell anywhere, the domain bounding
box, the largest diffusivity. It happens to be computed from rank-local data, so it
must be reduced across ranks, and a rank with no cells must contribute the identity
element of that reduction (+inf for a MIN, -inf for a MAX, 0 for a SUM) rather
than raise. Every rank then calls the collective and every rank — starved or not —
returns the same value. That is the shape get_mean_radius already used; the fix
copies it.

Two supporting rules the diff follows throughout:

• No rank-asymmetric raise around a collective. Where a short-circuit is right
  (points_in_domain on an empty rank is honestly all-False), the collective is
  taken first, by every rank, and the short-circuit comes after.
• A genuinely rank-local path may return early, because it takes no collective with
  it. The two interpolation guards are of this kind, and say so.

Sites fixed

From the issue's own list (all verified still present at development @ e475246f):

Site	Was	Now
discretisation_mesh.py — get_min_radius / get_max_radius	_radii.min() / .max() on an empty array	MPI allreduce MIN/MAX; empty rank contributes ±inf
discretisation_mesh.py — points_in_domain	rank-local short-circuit placed before the collective get_max_radius; then interrogated an empty cell set	collective first; then all-False on a zero-cell rank
utilities/_utils.py — gather_data	silently stripped NaN rows, so the table's row index no longer equalled rank	stripping is now strip_nan=True, off by default and documented
function/functions_unit_system.py — monotone kNN limiter	kd-tree over an empty source cloud, stencil gather on an empty array	returns the (empty) value unchanged; rank-local, no collective skipped
discretisation/discretisation_mesh_variables.py — rbf_interpolate	same, for the MeshVariable RBF rung	correctly-shaped zeros (matches the already-guarded SwarmVariable path)
meshing/smoothing/metrics.py — mesh_metric_mismatch	KNOWN LIMIT note said this would lift once (1) was fixed	measured: it does not — see below

Found beyond the list, same class (a global quantity computed from rank-local data):

• Mesh.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
  every reduction is empty-safe. (The quad-box fixtures escaped this by luck: both
  branches coincide there.)
• Mesh.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 a new _global_coord_bounds() which allreduces the box. Kept behaviour-neutral
  otherwise: a suspected double application of the length scale in that path is flagged
  with a TODO(BUG) rather than changed here.
• estimate_dt, three implementations: _global_max_diffusivity (.max() of the
  diffusivity sampled at this rank's centroids), SNES_NavierStokes.estimate_dt
  (.max() of centroid velocity magnitudes), and
  SNES_Stokes_SaddlePt.estimate_dt in the .pyx (.max() of the local velocity
  DOFs). All three reduced before the allreduce, so a starved rank raised while its
  peers waited. These are the priority path named in the issue.

Assessed and deliberately left:

• SNES_AdvectionDiffusion.estimate_dt and SNES_Stokes.estimate_dt per-element
  reductions — already len()-guarded, verified empty-safe.
• Swarm.estimate_dt — already has a sanctioned except (ValueError, IndexError)
  contributing 0.0 to the MAX.
• _get_domain_centroids — already guarded (Mesh construction crashes on empty ranks: navigation kd-tree build passes a 1-D empty array (StructuredQuadBox(elementRes=(2,2)) at np4) #399) with a finite sentinel. Only its
  comment changed, to record that gather_data no longer strips NaN but the sentinel
  is still required (a NaN row would poison the kd-tree).
• quality()'s percentiles and neighbour size-jump — rank-local estimates by
  documented design; they now return NaN on a rank with no local distribution rather
  than raising.
• mesh._radii per-cell loop in _get_mesh_sizes, stats() (PETSc Vec.max() is
  already global), the XDMF connectivity check (rank-0, file data).

Item 2 does not lift — measured, not assumed

The issue expected mesh_metric_mismatch's starved-rank limit on field-valued metrics
to lift automatically once the reduction layer was safe. It does not. Measured at np=4
on a starved mesh with this fix in place: get_min_radius, get_max_radius,
get_mean_radius, points_in_domain, quality() and estimate_dt all succeed and
agree, but any path that builds a mesh-variable sub-DM fails below UW3 with
MPI_ERR_BUFFER out of DMCreateSubDM_Plex → DMClone when a rank has no cells. The
KNOWN LIMIT comment now records that, names #314 as the remaining blocker, and says
what was measured.

Fail-before evidence

Same test file, same fixture, on a build with the fix stashed out:

• np=2 — test_radius_accessors_are_global_on_a_zero_cell_rank FAILED on one
  rank only with ValueError: zero-size array to reduction operation minimum which has no identity (raised at discretisation_mesh.py:6512, get_min_radius), the
  other rank did not fail, and the run then HUNG at the next test. Ended by the
  mpirun --timeout guard.
• np=4 — three ranks FAILED, one did not, then HUNG. Ended by the timeout
  guard.

That asymmetry is the bug: one rank raises, the rest are still in the collective.

With the fix: 7 passed at np=2 and 7 passed at np=4.

Tests

tests/parallel/test_0774_empty_rank_reductions_mpi.py (7 tests, level_1 +
tier_a, pytest.mark.timeout(120) — the pre-fix failure mode is a hang, so every
test is timeout-guarded).

The fixture is StructuredQuadBox(elementRes=(1, 2)): two cells, so PETSc must leave a
rank empty at both np=2 ([2, 0]) and np=4 ([0, 0, 0, 2]).

Two tests exist to keep the others honest:

• test_premise_some_rank_owns_no_cells fails if the partition is not actually
  starved, so the suite cannot quietly stop testing the thing it is named for.
• test_negative_control_rank_local_minimum_would_be_caught (house negative-control
  rule) computes the rank-local minimum on this partition and asserts the ranks do
  not agree — proving the cross-rank agreement assertion is not true by
  construction, and would fire if get_min_radius ever returned a rank-local answer.

Oracle for the value: the reported characteristic length is the cell's
centroid-to-corner half-diagonal, sqrt(0.5² + 0.25²) for this fixture — analytic, not
a recorded number. Matches the np=1 value to 1e-9.

Verification

Run	Result
test_0774_empty_rank_reductions_mpi.py np=2	7 passed
test_0774_empty_rank_reductions_mpi.py np=4	7 passed
ptest_0008_mesh_radii_accessors.py np=2 / np=4	OK, values unchanged from before the fix
test_0700_basic_parallel_operations.py np=2 / np=4	12 passed 2 skipped / 14 passed
test_0750_global_statistics.py np=2 / np=4	13 passed / 13 passed
test_0755_swarm_global_stats.py (a gather_data consumer) np=2 / np=4	11 passed / 11 passed
tests/test_1018_rotated_freeslip.py	22 passed
pytest tests -m "level_1 and tier_a" -q --ignore=tests/test_0050_utils.py	604 passed, 24 skipped, 1 xfailed, 0 failed

Unblocks #314

#314 (global_evaluate np=4 hang on empty-interior ranks) is the same family seen from
the DMInterpolation side, and could not be worked while the layer beneath it raised
first: any starved-rank reproducer died in get_min_radius before reaching
DMInterpolation. That layer is now clean, and the measurement above localises what is
left — DMCreateSubDM_Plex → DMClone returning MPI_ERR_BUFFER on a zero-cell rank —
which is a concrete starting point for #314 rather than a symptom. #512's item 1 should
fall out of the same work.

Underworld development team with AI support from Claude Code

[Copilot]

Pull request overview

This PR fixes MPI deadlocks and asymmetric failures when a rank owns zero cells by ensuring “global” mesh queries and reduction-based diagnostics use collective reductions with identity elements (rather than rank-local .min()/.max() on empty arrays), and by adding a targeted parallel regression suite for issue #405.

Changes:

• Make global reductions (min/max radius, max diffusivity, max velocity magnitude, global bounds) empty-rank safe by contributing identity elements and still participating in collectives.
• Reorder/guard evaluation helpers (points_in_domain, quality(), interpolation/limiter helpers) to avoid rank-asymmetric control flow around collectives.
• Add tests/parallel/test_0774_empty_rank_reductions_mpi.py to prevent regressions on starved partitions (np=2 and np=4).

Reviewed changes

Copilot reviewed 8 out of 8 changed files in this pull request and generated 2 comments.

Show a summary per file

File	Description
tests/parallel/test_0774_empty_rank_reductions_mpi.py	New MPI regression tests covering radii, points_in_domain, quality, estimate_dt, and gather_data rank-row preservation on empty ranks.
src/underworld3/utilities/_utils.py	Changes gather_data semantics to keep NaN rows by default; adds strip_nan opt-in and expands documentation.
src/underworld3/systems/solvers.py	Makes _global_max_diffusivity and Navier–Stokes estimate_dt reductions empty-safe.
src/underworld3/meshing/smoothing/metrics.py	Updates KNOWN LIMIT note with measured behavior and points remaining blocker to #314.
src/underworld3/function/functions_unit_system.py	Adds an empty-input early return in the monotone limiter to avoid indexing empty source clouds.
src/underworld3/discretisation/discretisation_mesh.py	Fixes empty-safe min/max radii; enforces collective ordering in points_in_domain; makes quality() branch choice collective; introduces _global_coord_bounds() for global physical bounds/extent.
src/underworld3/discretisation/discretisation_mesh_variables.py	Guards rbf_interpolate against empty DOF data by returning correctly shaped zeros.
src/underworld3/cython/petsc_generic_snes_solvers.pyx	Makes Stokes saddle-point estimate_dt velocity-max reduction empty-safe.

---

💡 Add a code-review agent skill or configure MCP servers for context-aware, tailored reviews. Learn more in the docs.

[lmoresi]

CI regression diagnosed and fixed in 6e97b076 — it was ours, and the mechanism is worth recording because it is a type contract that had only ever held by accident.

np.squeeze(np.float64(0.5)) returns np.float64(0.5); np.squeeze(0.5) returns array(0.5) — a 0-d ndarray. Every estimate_dt funnels its result through np.squeeze. Making get_min_radius return a plain float (what the allreduce naturally produces, and what its -> float annotation already promised) therefore promoted the timestep to a 0-d array, and solve(timestep=dt) evaluating timestep != self.delta_t against a UWexpression made numpy try to convert the sympy object to a sequence instead of deferring to it. Measured, merge base vs head: same value, different type (np.float64(0.01586914…) vs array(0.01586914…), ndim 0).

Fixed at the contract boundary rather than by restoring the accident: _apply_unit_aware_scaling already documents its return as "float or UWQuantity", so it now collapses a 0-d array to a scalar (unit-bearing values pass through untouched — .item() would discard units), and the two estimate_dt implementations that return directly use the same helper. Returned numbers are bit-identical to the merge base. tests/test_1007_estimate_dt_scalar_contract.py pins the contract so it no longer depends on a caller's numeric type — and it sits in the test_100[0-9] batch CI actually runs.

Incidental catch: SNES_AdvectionDiffusion.estimate_dt has been returning a 0-d array since before this branch — a live trap, now closed by the same helper.

Re-ran the batches this diff can reach that the tier-A gate does not (they are the exposure via solvers.py and the .pyx): test_100x 28 passed (was 1 failed), 101x/102x 137, 105x 83, 1100 5+1 xpassed, 1110/1120 4, 1450 3; plus test_0774 7 at np=2 and 7 at np=4, and the full tier-A gate 604 passed / 0 failed.

For a reviewer: the part of this PR to look at hardest is the radius accessors' return type, consumed at ~14 sites — that is what had the blast radius here.

Underworld development team with AI support from Claude Code