Make mesh.cell_size() local and partition independent

[gthyagi]

Summary

Fixes #569 and fixes #687 by making Mesh.cell_size() a per-cell geometric quantity that is independent of MPI partitioning.

The previous field used self._radii, whose nearest-centroid kd-tree contains only centroids held by the current rank. Near partition boundaries, the nearest available centroid can therefore change with rank count. Nitsche penalties and SUPG stabilization terms built from mesh.cell_size() consequently changed with the partition.

Change

• Compute each cell's characteristic length as the RMS distance of its vertices from that cell's own centroid.
• Read vertices from the DMPlex vertex stratum and coordinate-section offsets, covering triangles, tetrahedra, quadrilaterals, and hexahedra.
• Cache the values as _cell_radii; the name describes cell geometry and does not overload own/owned, which elsewhere denotes MPI-rank ownership.
• Use _cell_radii only for mesh.cell_size() and refresh its field after mesh deformation through the existing reinitialization path.
• Preserve legacy kd-tree radii for global timestep estimates, adaptivity, and mesh-relaxation consumers.
• Update the cell_size() and Nitsche test documentation to state the exact definition and the default Nitsche local-h consequence.

Why the scope is narrow

This does not redefine get_min_radius(), get_max_radius(), or get_mean_radius(). It corrects only the local field consumed by spatially varying stabilization and penalty terms. Follow-up issue #694 tracks the global accessor contract.

Regression coverage

• Independent vertex-coordinate oracle for simplex and tensor meshes in 2-D and 3-D.
• Initial and deformed geometries.
• Analytical regular-square control and preservation of legacy global radius values.
• Nitsche local-h magnitude, deformation refresh, and free-slip solve checks.
• Enumerated tests/parallel/ regression that gathers every owned-cell (centroid, cell radius) row and compares the complete sorted table with a fresh np=1 process using the same cached Gmsh mesh.

Validation

Run	Result	Pytest time
Focused serial geometry/Nitsche suite	9 passed	7.93 s
Cell-by-cell partition regression, np=2	passed	9.15 s
Cell-by-cell partition regression, np=4	passed	6.72 s
Cell-by-cell partition regression, np=8	passed	8.68 s

The Stokes case emits the existing expected GAMG fallback warning because its tiny mesh has no multigrid hierarchy; all assertions pass.

This is split from #689 so the generic mesh-size correction can be reviewed and merged independently of Eulerian SUPG transport.

Underworld development team with AI support from Claude Code

[lmoresi]

Reviewed by measurement. The fix works, exactly as claimed, and the defect
is real on current development. One gap worth closing before merge.

The fix does what it says

Per-cell values gathered and compared across rank counts on the same mesh
(198 cells, UnstructuredSimplexBox, cellSize=0.12), reading _radii_own
and _radii directly rather than sampling through evaluate():

field	max abs diff np1 vs np2	np1 vs np4
_radii_own — what cell_size() now uses	0.000e+00	0.000e+00
_radii — the legacy kd-tree field	3.261e-03	4.110e-03

So the new field is bit-identical across partitions and the old one is not.
That also reproduces the original defect independently, which is worth having
on the record.

The gap: the property is never checked at more than one rank count

tests/test_0010_cell_size_geometry.py runs in the serial batch
(scripts/test.sh line 85, tests/test_00[0-4]*py). The new tests are written
with uw.mpi.comm.allgather, so they would be meaningful under mpirun — but
nothing runs them there.

What they assert is a within-rank oracle: each cell's size matches its own
vertices. That is correct, and it is not the claim in the title. Partition
independence is a statement about two different rank counts agreeing, and no
test compares two rank counts.

This is the shape that has bitten this repo twice in the last month. #675 was
two files that matched no glob in either test script and failed the first time
CI ever ran them; #615 fixed the enumeration for tests/parallel/ specifically
so a file could not go unrun. A test that is about parallel behaviour and only
ever executes at np=1 is the same hole in a different place.

Concretely: a tests/parallel/ case that builds the mesh, gathers
(centroid, cell_size), and asserts the sorted set matches a stored np=1
reference — or simply that _radii_own is bit-identical to the serial
reference, which is what the table above measures. tests/parallel/ is
enumerated since #615, so it will actually run.

Two smaller points

#687 duplicates #569 (15 August), which is still open and carries something
this PR description does not: the consequence. add_nitsche_bc(local_h=True) is
the default, so Nitsche answers moved 6.6e-03 between rank counts. Worth
closing both, and worth stating the Nitsche consequence in the changelog — it
is the reason anyone cares.

Scope is right. Leaving get_min_radius() / get_max_radius() /
get_mean_radius() and the timestep, adaptivity and relaxation consumers on the
legacy field is the conservative call, and the PR says so plainly.

Nothing here blocks the approach. The fix is correct and I would merge it with a
parallel guard added.

[lmoresi]

Follow-up on "why the scope is narrow", since the obvious next question is
whether everything else built on the legacy radii moves too. Measured on
development (UnstructuredSimplexBox, cellSize=0.12, 198 cells):

quantity	np=1	np=2	np=4	np=8
per-cell _radii	—	max abs diff 3.3e-03	4.1e-03	—
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
SNES_Stokes.estimate_dt(), prescribed velocity	0.0545712806276149	identical	identical	identical

Timesteps are not currently affected, and the reason matters. Every
estimate_dt path reduces to a global minimum — solvers.py:947, 4745,
swarm.py, free_surface.py via get_min_radius(), and SNES_Stokes via
min(h_i/|v_i|) over the raw per-cell array. The minimum turned out to be
bit-stable across 1/2/4/8 here, so estimate_dt() is bit-identical.

I first measured a ~1e-5 spread in estimate_dt() and was ready to call it
partition dependence. It was not: re-running with the velocity written in
analytically rather than solved gives bit-identical answers at every rank count,
so that spread was the iterative solve, not the radii. Worth recording so nobody
else chases it.

But the insulation is incidental, not structural. The min is stable because
the minimising cell happens not to sit where the kd-tree lookup differs. A mesh
whose smallest cell lies on a partition boundary would move the timestep, and
nothing in the code prevents that.

What is already partition-dependent today, and is not fixed by this PR:

• get_max_radius() — moves at np=8, by 4.9%.
• get_mean_radius() — moves at every rank count.

So the narrow scope is defensible for cell_size(), but "preserve legacy
kd-tree radii for global timestep estimates, adaptivity, and mesh-relaxation
consumers" leaves those two returning partition-dependent numbers. Adaptivity
and relaxation consume them.

One more, minor: the follow_metric docstring in
src/underworld3/meshing/smoothing/api.py shows
gradient_smoothing_length=2.0 * mesh._radii.mean(). That is a rank-local
mean — a different value on every rank — so as a documented recipe it teaches
the partition dependence rather than avoiding it. Docstring only, not live code.

None of this needs to be in this PR. It is the follow-on issue, and #569 is
probably the right place for it.

[lmoresi]

Filed #694 for the other half, with your ruling recorded: get_*_radius() ought to be rank-independent because that is the contract they advertise, and anything local-only needs to say so in its name.

The fix is half-built here — _radii_own is partition-independent by construction, so building the three accessors on it makes all of them rank-count-independent and makes their docstrings true. Measured on this branch: per-cell own values are bit-identical at np=1/2/4 where the legacy field differs by 3.3e-03 and 4.1e-03.

Not asking for it in this PR. The scope here is right; #694 is the follow-on.

[lmoresi]

One naming request before this lands, and it is worth doing now because private
fields get copied.

_radii_own reads as the opposite of what it is. In this codebase own /
owned consistently means rank ownership:

swarm.py:3936             "keeping only locally-owned points"
petsc_compat.h:297        "sum local owned contributions across all ranks"
petsc_generic_snes_solvers.pyx:3032  "ghost ones and owned by a different rank"
petsc_generic_snes_solvers.pyx:3579  "closure DOFs that are non-owned on this rank"

So _radii_own parses as "the radii of the cells this rank owns" — a partition
concept — when it is the one field here that is partition-independent by
construction. That is an unusually costly name for a field whose entire purpose
is the distinction it obscures.

Neither field is about rank ownership. Both are per-cell local arrays. The
real distinction is where the number comes from:

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

A nomenclature that says that, per Charter §3 ("names state what a thing IS"):

• _radii_own → _cell_radii — the cell's own geometry, nothing else in it
• _radii → _neighbour_radii — says where the number comes from, and
  implies why it moves with the partition

_cell_radii / _neighbour_radii reads correctly at every call site and leaves
own meaning what it means everywhere else. Renaming the legacy field is a
larger diff and could reasonably be split out — but at minimum the new one
should not be called _radii_own.

Whatever is chosen, #694 will inherit it: the fix there is pointing
get_min/max/mean_radius() at this field, so the name lands in three more
call sites the moment that is done.

[gthyagi]

Addressed the requested merge gap in 6fde1ac9.

• Renamed the new geometry-derived cache from _radii_own to _cell_radii; own no longer conflicts with the codebase's rank-ownership terminology.
• Added tests/parallel/test_1077_cell_size_partition_independence.py, which gathers owned-cell (centroid, cell radius) rows and compares the complete sorted table against a fresh np=1 run on the same cached Gmsh mesh.
• The test is under the enumerated tests/parallel/ suite, so CI cannot omit it through a filename glob.

Local validation after rebuilding the branch:

Run	Result	Time
Focused serial geometry/Nitsche suite	9 passed	7.93 s
New partition regression, np=2	passed	9.15 s
New partition regression, np=4	passed	6.72 s
New partition regression, np=8	passed	8.68 s

This directly tests the title claim across rank counts; the existing independent vertex-geometry oracle remains as the separate correctness/deformation check. The practical Nitsche consequence remains documented in Mesh.cell_size(): default local-h penalty scaling no longer changes with partition-local neighbour spacing.

[gthyagi]

For sequencing context, this fix is also the direct cause of #689's current sole CI failure (test_error_is_partition_independent, absolute difference 1.823e-08 against a 1e-08 gate). Once #691/#692 land, the Eulerian SUPG target branch can be synced and #689 rerun without carrying either generic mesh fix in its feature diff.

[gthyagi]

Thanks. Consolidated answers and the resulting changes:

1. Cross-rank claim: addressed in 6fde1ac9. The new enumerated tests/parallel/test_1077_cell_size_partition_independence.py compares the complete sorted owned-cell (centroid, cell radius) table with a fresh np=1 process on the same mesh. It passes at np=2, 4 and 8 with atol=1e-14.

2. Naming: agreed. The new field is _cell_radii, not _radii_own; own/owned remains reserved for MPI ownership. I have not renamed legacy _radii in this PR because its remaining users and eventual deletion belong to the coordinated Cell size: one definition, one name — three are in play, the docstrings describe a fourth, and get_*_radius() is partition-dependent #694 cleanup.

3. 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: agreed that they are the same defect. The PR description now says Fixes #569 and fixes #687, so this PR will close both.

4. Default Nitsche consequence: agreed, and now recorded in docs/developer/CHANGELOG.md: the old rank-local centroid kd-tree moved the default local_h=True velocity answer by 6.6e-3 between rank counts. Commit e6eaac28 also removes the local_h=False workaround from test_1069; the test now exercises the public default and matches its fresh serial reference at np=2, 4 and 8 (10.99 s, 7.31 s and 9.60 s).

5. Global radius accessors and timestep: this PR does not claim to fix them. get_min/max/mean_radius() still reduce legacy _radii; that makes each call rank-symmetric but not rank-count independent. The observed get_min_radius()/estimate_dt() stability is incidental, as you measured. The follow_metric example using mesh._radii.mean() is likewise still wrong and belongs to Cell size: one definition, one name — three are in play, the docstrings describe a fourth, and get_*_radius() is partition-dependent #694.

6. PETSc FVM question: neither legacy _radii nor this PR's _cell_radii uses DMPlexComputeGeometryFVM. This PR's value is the RMS distance of a cell's vertices from its own centroid. The evidence on Cell size: one definition, one name — three are in play, the docstrings describe a fourth, and get_*_radius() is partition-dependent #694 shows the PETSc volume**(1/dim) route now returns one finite, partition-identical value for simplex/tensor cells in 2-D/3-D and after deformation. I agree that it is the better final canonical definition: Cell size: one definition, one name — three are in play, the docstrings describe a fourth, and get_*_radius() is partition-dependent #694 should make that the sole _cell_radii, point cell_size() plus all three accessors at it, remove the kd-tree/dead arrays, and fix the follow_metric example. I am keeping that larger semantic cleanup separate because the comments here explicitly treated it as a follow-on rather than a blocker for Make mesh.cell_size() local and partition independent #692.

The branch is pushed at e6eaac28; GitHub CI has restarted for that head.