Cell size: one definition, one name — three are in play, the docstrings describe a fourth, and get_*_radius() is partition-dependent

[lmoresi]

The defect

get_min_radius(), get_max_radius() and get_mean_radius() present as the
canonical global "mesh length" API. Their docstrings say so explicitly:

• get_min_radius: "the smallest cell anywhere in the mesh, not just on this
  rank. Parallel-safe via MPI allreduce."
• get_max_radius: "the largest cell anywhere in the mesh."
• get_mean_radius: "Use this anywhere you need a representative h0 ...
  rather than reaching for the rank-local self._radii array, which gives
  different answers on different MPI ranks and leaks downstream."

All three reduce over self._radii — the array that last docstring warns
against. The allreduce makes the answer identical on every rank. It does not
make it identical for every rank count, because the values being reduced are
themselves partition-dependent: _radii comes from a kd-tree over this rank's
centroids, so a cell near a partition boundary gets a different size depending
on how the mesh was split.

"Same on every rank" and "same for every rank count" are different properties.
These functions deliver the first and their names and docstrings promise the
second.

Measured

UnstructuredSimplexBox, cellSize=0.12, 198 cells, on development:

accessor	np=1	np=2	np=4	np=8
get_min_radius()	0.0502908499358	same	same	same
get_max_radius()	0.0670934714626	same	same	0.070373950134
get_mean_radius()	0.0592989430362	0.0594222426856	0.0594775301273	0.0597593860929

get_max_radius() moves by 4.9% at np=8. get_mean_radius() moves at every
rank count. get_min_radius() is stable here — but incidentally, because the
minimising cell does not happen to sit where the kd-tree lookup differs, not
because anything guarantees it.

Ruling (maintainer, 2026-09-05)

These ought to be rank-independent — that is part of the parallel contract they
advertise. If a quantity is local-only, its name has to say so.

The fix is already half-built

#692 adds _radii_own: each cell's characteristic length from its own
vertices and centroid, which cannot depend on the partition. Measured on that
branch, per-cell values are bit-identical at np=1/2/4 (max abs diff 0.000e+00)
where the legacy field differs by 3.3e-03 and 4.1e-03.

Building the three accessors on _radii_own makes all of them genuinely
rank-count-independent and makes the docstrings true.

Two things to get right when doing it:

1. min and max become exactly reproducible; mean needs care. A
   distributed sum reduces in partition order, so even over identical per-cell
   values the mean can differ in the last bits between rank counts. If
   bit-identity is wanted rather than agreement to ~1e-16, the reduction has to
   be made order-independent.
2. Timesteps are currently insulated, and that should not be relied on.
   Every estimate_dt path reduces to a global minimum, and estimate_dt()
   measures bit-identical at np=1/2/4/8. It is stable for the same incidental
   reason get_min_radius() is.

Also

follow_metric's docstring in src/underworld3/meshing/smoothing/api.py
recommends gradient_smoothing_length=2.0 * mesh._radii.mean() — a rank-local
mean, a different value on every rank. Docstring only, but it teaches exactly the
pattern get_mean_radius's docstring warns against.

Related: #569 and #687 (the same defect in mesh.cell_size(), fixed by #692).

[lmoresi]

Correction to point 1 above — I overstated it.

get_mean_radius() uses a numpy local sum with an MPI.SUM allreduce of both the sum and the count:

local_sum = float(radii.sum())
local_n   = int(radii.size)
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)

That is a true global mean over all cells, not a mean-of-means, and the partition-order last-bit difference is not worth engineering around (maintainer, 2026-09-05). A PETSc vector reduction would reduce in partition order too, so it buys nothing.

So this issue is simpler than I wrote it. All three accessors already reduce correctly; they just reduce over the wrong array. Swapping self._radii for #692's _radii_own fixes all three, and no change to the reduction machinery is needed. The 4.9% jump in get_max_radius() at np=8 is entirely the input.

Point 2 stands: estimate_dt() is currently safe only because every path takes a global min, and that stability is incidental.

[lmoresi]

Naming, raised on #692 and recorded here because this issue inherits it.

own / owned means rank ownership everywhere else in this codebase ("locally-owned points", "owned by a different rank", "non-owned on this rank"), so #692's _radii_own parses as "the cells this rank owns" — the opposite of what it is, which is the one field that is partition-independent by construction.

Neither field is about rank ownership; both are per-cell local arrays. The distinction is provenance:

field	derived from	partition-dependent
_radii	kd-tree over this rank's centroids	yes
_radii_own	this cell's own vertices about its own centroid	no

Suggested: _cell_radii and _neighbour_radii.

This matters here because the fix for this issue is pointing get_min/max/mean_radius() at that field, which puts the name in three more call sites — in the very functions whose contract is about partition independence.

[lmoresi]

The PETSc FV route works, and it is the definition the docstrings already claim

Maintainer's question: is one of these the PETSc FV element-radius routine? It
was tried, errors were found, and it was abandoned — but it is the obvious place
to get this.

Neither field is that routine today. petsc_fvm_get_local_cell_sizes
(cython/petsc_discretisation.pyx) wraps DMPlexComputeGeometryFVM and returns
volume**(1/dim). It is commented out at its call site in
discretisation_mesh.py and marked "Deprecated in favour of _get_mesh_sizes",
still carrying the note:

the petsc4py version of DMPlexComputeGeometryFVM does not compute all cells
and does not obtain the minimum radius for the mesh

That no longer reproduces. Same mesh as the measurements above
(UnstructuredSimplexBox, cellSize=0.12, 198 cells):

np	local cells	values returned	non-finite or zero	global mean
1	198	198	0	0.0708408414261
2	99	99	0	0.0708408414261
4	50	50	0	0.0708408414261

One value per local cell, all finite and positive, and per-cell exactly
equal across rank counts:

max|np1 - np2| = 0.0
max|np1 - np4| = 0.0

A cell's volume is intrinsic, so volume**(1/dim) cannot depend on the split.
It is partition-independent by construction rather than by careful
reimplementation.

This also fixes a documentation bug

The docstrings already describe the PETSc definition, not what the code
computes:

• cell_size(): "the volume**(1/dim) equivalent radius, i.e. self._radii"
• get_mean_radius(): "volume^(1/dim), i.e. the equivalent radius derived
  from each cell's volume"

_radii is not that. _get_mesh_sizes builds a kd-tree over this rank's
centroids, queries it with each cell's vertices, and takes
cell_r = sqrt(distsq.mean()) — an RMS vertex-to-nearest-centroid distance.
Near a partition boundary the true nearest centroid can belong to another rank
and be missing from the tree, which is the drift measured above. So the
docstrings describe the abandoned routine and the code does something else.

What this implies for #692

#692 defines _radii_own as the RMS distance of a cell's vertices from its
own centroid — a third definition, reimplementing by hand a partition-
independent cell size that PETSc already computes correctly.

Worth deciding which is wanted before the naming is settled, because the two are
not the same number: FVM mean here is 0.0708408 against _radii_own's 0.0627204.
Both are defensible characteristic lengths; only one matches what the docstrings
have been promising, and only one avoids a hand-rolled geometric routine that
has to be maintained for triangles, tets, quads and hexes.

If DMPlexComputeGeometryFVM is used, the naming question largely dissolves —
there is one per-cell radius, it is volume**(1/dim), and it is partition
independent. _radii as it stands is then a neighbour-spacing measure, which
is a different quantity that deserves its own name wherever it is still wanted.

Caveat, stated plainly: this is one mesh type, 2-D simplex, at np ≤ 4. Before
acting on it the same check is wanted on quads/hexes, in 3-D, and on a
deformed mesh — the original complaint may have come from a case none of this
touches.

[lmoresi]

Consolidated proposal — one cell size, one name

Superseding the scattered comments above. Everything here is measured on
development unless marked otherwise; the caveats are at the end.

There are three definitions of "cell radius" in play

#	where	how	partition-independent	mean, same mesh
1	_radii (live)	kd-tree over this rank's centroids, queried with each cell's vertices; sqrt(distsq.mean())	no	0.0592989
2	_radii_own (#692)	RMS of the cell's vertices about its own centroid	yes	0.0627204
3	petsc_fvm_get_local_cell_sizes (abandoned)	DMPlexComputeGeometryFVM → volume**(1/dim)	yes	0.0708408

Three different numbers for the same words. Definition 1 is what the code
computes; definition 3 is what the docstrings say it computes:

• cell_size(): "the volume**(1/dim) equivalent radius, i.e. self._radii"
• get_mean_radius(): "volume^(1/dim), i.e. the equivalent radius derived
  from each cell's volume"

Neither is true of _radii. That is the first loose end: the documentation has
been describing the abandoned routine for as long as it has been abandoned.

The abandoned routine works

The note on it reads "does not compute all cells and does not obtain the minimum
radius for the mesh". Measured (2-D simplex, 198 cells):

np	local cells	values returned	non-finite or zero
1	198	198	0
2	99	99	0
4	50	50	0

Per-cell, across rank counts: max|np1-np2| = max|np1-np4| = 0.0. Volume is
intrinsic, so this is partition-independent by construction rather than by
careful reimplementation.

Two of the four arrays are dead

_get_mesh_sizes returns (_min_size, _radii, _centroids, _search_lengths).
_min_size and _search_lengths are assigned and never read anywhere in
src/ or tests/ — a per-cell kd-tree query computes three distance
statistics and two are discarded. _centroids does not need the loop either:
_get_mesh_sizes gets them from _get_coords_for_basis(0, False) before
building the tree.

Proposal

1. One per-cell field, from PETSc: volume**(1/dim) via
   DMPlexComputeGeometryFVM, named _cell_radii.
2. cell_size() and get_min/max/mean_radius() all read it. All four
   become partition-independent, and the docstrings become true without being
   rewritten. The reduction machinery in the three accessors is already correct
   and does not change.
3. Delete the per-cell kd-tree loop, _min_size and _search_lengths.
   Keep _centroids from the basis coords.
4. If a neighbour-spacing measure is wanted anywhere, it is a different
   quantity and gets its own name — _neighbour_spacing — and its own
   justification at the call site. Nothing currently asks for it: cell_size()
   feeds Nitsche penalties and SUPG stabilisation, which want a cell length h,
   and volume**(1/dim) is the standard choice for that.

This makes own unnecessary as a qualifier, which resolves the naming problem
rather than renaming around it — own/owned goes back to meaning rank
ownership everywhere, as it does in the rest of the codebase.

What this costs, stated plainly

It changes answers. cell_size() moves from 0.0592989 to 0.0708408 in the
mean — about 20% — and it feeds add_nitsche_bc(local_h=True), which is the
default, and SUPG stabilisation. That is a physics-affecting change and
needs the Nitsche and SUPG benchmarks run before and after, not just unit tests.
#692's definition 2 has the same problem at about 6%; there is no
partition-independent choice that leaves the numbers alone, because the current
numbers are the partition-dependent ones.

The measurements are narrow. 2-D simplex, UnstructuredSimplexBox, np ≤ 4,
undeformed. Before acting: quads and hexes, 3-D, and a deformed mesh. The
original complaint about DMPlexComputeGeometryFVM may have come from a case
none of this touches — DMPlexGetMinRadius is a different call and the note
conflates the two.

If the FVM route is not taken

#692 stands on its own and should merge: it fixes cell_size(), which is the
user-visible half, and its scope is right. The requests on it are then unchanged
— a guard that runs at more than one rank count, and a name that is not
_radii_own. This proposal is the larger tidy, not a blocker for that.

Loose ends this closes

• mesh.cell_size() is partition-dependent (kd-tree of this rank's centroids) — and add_nitsche_bc(local_h=True) is the DEFAULT, so Nitsche answers move 6.6e-03 between rank counts #569 and mesh.cell_size() is partition-dependent: the per-cell radius comes from a kd-tree over the rank's centroids, so cells at partition boundaries get a different size on each rank count #687 are the same defect, filed twice.
• Docstrings describe a definition the code does not use.
• Three definitions of one quantity.
• own colliding with rank ownership.
• Two dead arrays and a kd-tree loop computing statistics nobody reads.
• follow_metric's docstring recommending 2.0 * mesh._radii.mean(), a
  rank-local mean, in the same breath as get_mean_radius() warns against
  exactly that.

[lmoresi]

Ruling (maintainer, 2026-09-05)

We should not keep broken code for backward compatibility with tests. We have
to own the mistake.

Recorded, and it changes the shape of the proposal above: the ~20% move in
cell_size() is not a cost to be weighed. The current numbers are the
partition-dependent ones — anything validated against them was validated against
a quantity that depends on how many ranks the run was given.

I also framed the cost wrong, and it is smaller than I implied. I said this
"needs the Nitsche and SUPG benchmarks run before and after". Checking what the
guards actually assert:

test	what it pins
test_1065_nitsche_local_h	vrms_ref is computed in-test from an exact essential free-slip solve; Nitsche must match it within 2%. The rest are ratios and inequalities (gfmax/gfmin > 3.0, h_top_min > 2.0*g, leak < 1e-3).
test_1060_nitsche_freeslip	"Nitsche differs from essential by ..." < 0.01 — the same oracle.
test_1064_constrained_spherical_shell_response	published Zhong values within 5%.

None of these stores a cell_size-derived number. They compare against an exact
solution computed alongside, or against the literature. A better cell_size()
should leave them passing, and if it does not, that is a finding rather than a
re-baselining chore.

There is already a workaround in the tree that this fix should delete.
tests/parallel/test_1069_boundary_normal_parallel.py sets local_h=False
specifically to avoid this defect, and records the numbers in its docstring:

local_h=False is deliberate and it is not a workaround for this fix. The
default local_h=True scales the Nitsche penalty by mesh.cell_size()
[...] partition-dependent in its own right (on this mesh the field's sum is
26.0822 at np=1, 26.1211 at np=2, 26.1386 at np=4 [...])

So a partition-independence test is currently forced to disable the default
code path in order to measure anything. That workaround is the debt marker:
removing it, and letting the test run with local_h=True, is part of the fix
and is the best single check that the fix worked.

Revised plan, replacing the "what this costs" section above:

1. One _cell_radii from DMPlexComputeGeometryFVM; cell_size() and all
   three get_*_radius() read it; delete the kd-tree loop, _min_size and
   _search_lengths.
2. Delete the local_h=False workaround in test_1069 and let it exercise the
   default.
3. Run the Nitsche and SUPG tests. They are oracle-based, so they should pass
   unchanged; a failure means the new definition is wrong, not that a baseline
   moved.
4. Widen the FVM check first — quads, hexes, 3-D, deformed — since everything
   measured so far is 2-D simplex at np ≤ 4.

The scope caveat stands: I have audited the three Nitsche tests above and
test_1069. I have not audited SUPG's guards, or anything consuming
get_max_radius() / get_mean_radius() in adaptivity and relaxation.

[lmoresi]

Decision tree and work estimate

Gate: does the PETSc route work? — ANSWERED, yes

The one thing that decided this. DMPlexComputeGeometryFVM, every cell family
we ship, 2-D and 3-D, undeformed and deformed:

case	cells	non-finite/zero	max abs diff np1 vs np2
2-D simplex	66	0	0.000e+00
2-D quad	36	0	0.000e+00
3-D simplex (tet)	184	0	0.000e+00
3-D hex	64	0	0.000e+00
2-D simplex, deformed	66	0	0.000e+00

One finite value per local cell everywhere, bit-identical across rank counts,
and it tracks deformation (mean 0.122822957 → 0.122822102). The "does not
compute all cells" note does not reproduce in any of them.

The tree

Adopt volume**(1/dim) from DMPlexComputeGeometryFVM as THE cell size?
│
├─ YES  ──►  _cell_radii is the single per-cell field.
│            · cell_size() and get_min/max/mean_radius() read it
│            · kd-tree loop, _min_size, _search_lengths deleted
│            · _radii_own (#692) unnecessary
│            · naming problem dissolves — no `own` qualifier needed
│            · docstrings become true without being rewritten
│
└─ NO   ──►  #692's _radii_own is the field.
             · still need: the accessors pointed at it (#694's core)
             · still need: a name that is not `_radii_own`
             · still need: a guard that runs at >1 rank count
             · hand-rolled geometry to maintain for tri/tet/quad/hex
             · docstrings still wrong, or must be rewritten

The NO branch is strictly more work and leaves more debt. The only argument for
it is if the FVM definition is wrong for the purpose — but volume**(1/dim)
is the standard element length for Nitsche and SUPG scaling, and it is what the
docstrings have claimed all along.

What needs updating, and how much

No call site changes. The accessors keep their names and signatures; only
the array behind them changes. That is what makes this small.

#	change	size
1	_cell_radii from petsc_fvm_get_local_cell_sizes, replacing the _get_mesh_sizes kd-tree loop	~40 lines net deleted in discretisation_mesh.py
2	cell_size() fill reads _cell_radii	1 line
3	get_min/max/mean_radius() read _cell_radii	3 lines
4	delete _min_size, _search_lengths (assigned, never read)	2 lines + the unpack
5	keep _centroids from _get_coords_for_basis(0, False)	already how it is obtained
6	delete the local_h=False workaround in test_1069	1 line + its docstring paragraph
7	new partition-independence test in tests/parallel/	~40 lines
8	fix follow_metric docstring (2.0 * mesh._radii.mean())	1 line

Call sites that do not change: cell_size() 11 in src / 6 in tests;
get_min_radius() 19 / 16; get_max_radius() 2 / 4; get_mean_radius() 1 / 3.

Estimate: half a day for 1-8, most of it writing the test rather than the
change. The change itself is a net deletion.

The risk is the answers, not the diff

cell_size() moves ~20% in the mean and feeds add_nitsche_bc(local_h=True),
the default, plus SUPG. Per the maintainer's ruling, that is not a reason to
keep the broken definition — the current numbers are the partition-dependent
ones. But it is where the time goes if something is wrong.

Validation, in order of what would actually catch a mistake:

1. test_1069 with the workaround removed — the single best check. It is a
   partition-independence test currently forced to disable the default path.
2. The three Nitsche tests — all oracle-based (test_1065 and test_1060
   compare against an exact essential free-slip solve computed in-test;
   test_1064 against published Zhong values within 5%). They should pass
   unchanged; a failure is a finding, not a re-baseline.
3. SUPG guards — not yet audited.
4. Adaptivity and relaxation — the consumers of get_max_radius() and
   get_mean_radius(), the two that are demonstrably partition-dependent today.
   Not yet audited.

Items 3 and 4 are the honest unknowns. Everything else above is measured.