Take the cell radius from PETSc, and fix the radius accessors with it (#694) by lmoresi · Pull Request #697 · underworldcode/underworld3

lmoresi · 2026-09-06T19:40:12Z

Summary

Takes the per-cell characteristic length from PETSc's
DMPlexComputeGeometryFVM — volume**(1/dim) — instead of a kd-tree over this
rank's centroids. A cell's volume is a property of that cell, so this cannot
depend on how the mesh was split.

Builds on #692, which fixed cell_size() with a hand-rolled own-centroid RMS
and established the _cell_radii name.

Why PETSc rather than the own-vertex RMS

	#692	this
`cell_size()` partition independent	yes	yes
`get_min/max/mean_radius()`	left on the kd-tree field	fixed
docstrings	still describe a third definition	become true
geometry to maintain for tri/tet/quad/hex	hand-rolled	PETSc's
kd-tree loop, `_min_size`, `_search_lengths`	still present	deleted

get_max_radius() moved 4.9% at np=8 and get_mean_radius() at every rank
count, while their docstrings advertise "the smallest cell anywhere in the mesh"
and "use this rather than the rank-local array". An allreduce makes an answer
agree on every rank; it does not make it agree at every rank count when the
values being reduced are themselves partition-dependent.

And both cell_size() and get_mean_radius() have always documented "the
volume**(1/dim) equivalent radius" — which is what this restores.

The abandoned routine works

It carried a note that it "does not compute all cells". That does not reproduce:

case	cells	max abs diff np1 vs np2
2-D simplex	66	0.000e+00
2-D quad	36	0.000e+00
3-D tetrahedra	184	0.000e+00
3-D hexahedra	64	0.000e+00
2-D simplex, deformed	66	0.000e+00

The note also named DMPlexGetMinRadius, a different call that is not used here.

Accessors after the change, identical to 12 digits at np=1/2/4/8:

min=0.0575495804241  max=0.0868948034925  mean=0.0708408414261

Two constants move with the length

The definition of h changed, so constants calibrated against it change too, or
the methods are silently different. Both were found by tests failing, not by
inspection.

Nitsche gamma: 10.0 → 12.5 (maintainer's call). The penalty is
gamma*mu/h and h is ~19% larger in the mean:

development, gamma=10,   old h    8.954e-05
this branch, gamma=10,   new h    1.312e-04     <- past test_1060's 1e-4
this branch, gamma=12.5, new h    7.851e-05

test_1060 passed gamma=10.0 explicitly, so the default alone did not reach
it; it now omits gamma and follows the default, because what it asserts is that
Nitsche enforces v.n = 0 at the recommended penalty.

max_surface_cfl: 0.5 → 0.437. The free-surface step is
max_surface_cfl * get_min_radius() / speed, and get_min_radius() grew by
1.1443, so the surface took ~14% larger steps and test_1070's strong-constraint
volume leak reached 1.063e-03 against a 1e-03 threshold. 0.5 / 1.1443 = 0.437
preserves the step.

Tests

tests/parallel/test_1078 guards the accessors — the half #692 does not reach —
by comparing against this module's own np=1 answer through serial_reference,
with the mesh fingerprint asserted so a host that triangulates differently is
reported as that rather than as partition dependence. #692's test_1077 keeps
covering the field cell by cell; the two are complementary and both stay.

Negative control: restoring the kd-tree field makes it fail with
max radius is partition dependent — np=1 0.0670934714626447 vs np=8 0.07037395013400187.

test_0010's oracle is rewritten, since it asserted the own-vertex RMS this
supersedes. It now checks volume**(1/dim) against a determinant volume from
the cell's vertices (simplices) and against the analytic cell volume scaled by
the affine map's determinant (structured boxes) — computed from coordinates, not
from PETSc, so it still checks the implementation rather than restating it.

test_1069's local_h=False workaround is removed. That test measures partition
independence and had to disable the default code path to avoid measuring this
defect at the same time; it now runs the default.


test_1060, 1065, 0010, 0008, 1070	26 passed, 1 skipped
parallel test_1077 + test_1078	3 passed at np=2 and np=4
test_10, test_0008, test_005* (before the CFL fix)	782 passed, 1 failed — that failure is test_1070, fixed here

Also

follow_metric's docstring recommended 2.0 * mesh._radii.mean(), a
rank-local mean, in the same breath as get_mean_radius() warns against
exactly that. It now recommends get_mean_radius().
The TODO(BUG) on _assemble_cell_size describing this defect is removed.

Underworld development team with AI support from Claude Code

🤖 Generated with Claude Code

https://claude.ai/code/session_01E87Q7KrpapxeQiLD1RiNXv

…tition (#694) mesh.cell_size() and get_min/max/mean_radius() were built on a kd-tree over THIS RANK's centroids, queried with each cell's vertices. Near a partition boundary the true nearest centroid can belong to a cell owned by another rank and be absent from the tree, so the answer moved with the rank count: per-cell by 3.3e-03 at np=2, get_max_radius() by 4.9% at np=8, get_mean_radius() at every rank count. It reached users through the DEFAULT add_nitsche_bc( local_h=True), which scales the penalty by cell_size(). An allreduce made those accessors agree on every rank without making them agree at every rank count -- the values being reduced were themselves partition dependent -- while their docstrings advertise "the smallest cell anywhere in the mesh" and "use this rather than the rank-local array". The radius is now PETSc's volume**(1/dim) from DMPlexComputeGeometryFVM. A cell's volume is a property of that cell, so this cannot depend on the split. That routine had been abandoned with a note that it "does not compute all cells"; that does not reproduce. Measured on 2-D simplex, 2-D quad, 3-D tetrahedra, 3-D hexahedra and a deformed mesh: one finite positive value per local cell every time, and bit-identical across rank counts in all five (max abs diff 0.000e+00). The old note also named DMPlexGetMinRadius, which is a different call and is not used here. It is also the definition the docstrings have claimed all along -- cell_size() and get_mean_radius() both describe "the volume**(1/dim) equivalent radius", which the kd-tree field was not. The documentation becomes true rather than being rewritten. Consequences beyond the swap: - `_min_size` and `_search_lengths` are gone. They were unpacked from _get_mesh_sizes and never read anywhere in src/ or tests/ -- the loop computed three distance statistics per cell and discarded two. - test_1069's `local_h=False` workaround is removed. That test measures partition independence and had to disable the DEFAULT code path to avoid measuring this defect at the same time; it now runs the default and passes (7 passed, 1 skipped at np=2). - follow_metric's docstring recommended `2.0 * mesh._radii.mean()`, a RANK-LOCAL mean, in the same breath as get_mean_radius() warns against exactly that. It now recommends get_mean_radius(). tests/parallel/test_0798 is the guard, and it compares two RANK COUNTS rather than asserting a within-rank property: the reference is this module's own np=1 answer via serial_reference, with the mesh fingerprint asserted so a host that triangulates differently is reported as that. Negative control: restoring the kd-tree field makes it fail with "max radius is partition dependent -- np=1 0.0670934714626447 vs np=8 0.07037395013400187". Verified: accessors identical to 12 digits at np=1/2/4/8 (max and mean both moved before); test_0798 passes at np=2/4/8; test_1069 passes on the default path. KNOWN, UNRESOLVED: the Nitsche leak in test_1060 grows from 8.954e-05 to 1.312e-04, crossing that test's 1e-4 threshold. This is not a defect in the new radius -- the penalty is gamma*mu/h and h is ~19% larger in the mean, so the same gamma enforces less. gamma's default of 10.0 was calibrated against the old definition and has to move with the quantity it scales, or the method is silently weaker. That is a solver-behaviour decision and is left for the maintainer rather than resolved by loosening a threshold. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01E87Q7KrpapxeQiLD1RiNXv

…of h (#694) The penalty is gamma*mu/h, so gamma is calibrated against whatever h means. mesh.cell_size() is now PETSc's volume**(1/dim) rather than a kd-tree distance to a neighbouring centroid, about 19% larger in the mean, so gamma=10.0 enforced correspondingly less: the free-slip leak in test_1060 grew from 8.954e-05 to 1.312e-04, past that test's 1e-4 threshold. A calibration constant has to move with the quantity it scales, or the method is silently weaker. gamma=12.5 (maintainer's call) restores it and slightly improves on the original: development, gamma=10, old h 8.954e-05 this branch, gamma=10, new h 1.312e-04 this branch, gamma=12.5, new h 7.851e-05 test_1060 passed gamma=10.0 explicitly, so the default alone did not reach it. It now omits gamma and follows the default, because what it asserts is that Nitsche enforces v.n = 0 at the RECOMMENDED penalty -- not that 10.0 in particular does. Pinning the old number there would have kept the test green while leaving every user on a weaker constraint, which is the failure mode this whole change exists to remove. test_1065's helper default moves with it. tests/test_1060 + test_1065: 9 passed (was 4 failed, 5 passed). Underworld development team with AI support from Claude Code Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01E87Q7KrpapxeQiLD1RiNXv

…own-vertex RMS #692 fixed cell_size() by computing each cell's RMS vertex distance about its own centroid. This takes the same quantity from DMPlexComputeGeometryFVM instead -- volume**(1/dim), which PETSc already computes and which cannot depend on the partition because a cell's volume is a property of that cell. What that adds over #692: - get_min/max/mean_radius() are fixed too. #692 left them on the kd-tree field, where get_max_radius() moved 4.9% at np=8 and get_mean_radius() at every rank count, while their docstrings advertise "the smallest cell anywhere in the mesh". - The docstrings become true rather than being rewritten: cell_size() and get_mean_radius() have always described "the volume**(1/dim) equivalent radius", which neither the kd-tree field nor the own-vertex RMS was. - The kd-tree loop, `_min_size` and `_search_lengths` are gone. The last two were unpacked and never read anywhere in src/ or tests/. - No hand-rolled geometry to maintain across triangles, tetrahedra, quads and hexahedra. Conflicts resolved toward the FVM definition, except where development's side was better on its own terms: test_1065 now uses the Charter's `.array[:, 0, 0]` access rather than the deprecated `.data`, which came from #692. #692's rename to `_cell_radii` is kept -- it answers the review request that `_radii_own` collided with `own` meaning rank ownership everywhere else. The two parallel guards are complementary and both stay: test_1077 (#692) compares the field cell by cell against its own serial answer; test_1078 (here, renamed from test_0798) does the same for the three accessors, which is the half #692 does not reach. test_0010's oracle is rewritten. It asserted the own-vertex RMS, which this supersedes; it now checks volume**(1/dim) against a determinant volume computed from the cell's vertices for simplices, and against the analytic cell volume -- scaled by the determinant of the affine map when deformed -- for the structured boxes. Both are computed from coordinates rather than from PETSc, so the test still checks the implementation instead of restating it. max_surface_cfl moves 0.5 -> 0.437. The free-surface timestep is max_surface_cfl * get_min_radius() / speed, and get_min_radius() grew by a factor of 1.1443 (0.0502908 -> 0.0575496), so the surface took ~14% larger steps and test_1070's strong-constraint volume leak grew to 1.063e-03 against its 1e-03 threshold. Same principle as the gamma change: a constant calibrated against a length has to move with the length, or the scheme is silently different. 0.5 / 1.1443 = 0.437 preserves the step. tests/test_1060, test_1065, test_0010, test_0008, test_1070: 26 passed, 1 skipped. tests/parallel/test_1077 + test_1078: 3 passed at np=2 and np=4. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01E87Q7KrpapxeQiLD1RiNXv

Copilot

🟡 Changes recommended

The updated simplex-volume oracle in tests/test_0010_cell_size_geometry.py computes closures using incorrect cell point IDs and unstable vertex selection, and there are misleading/outdated comments in the updated radius accessor implementations.

Once you've addressed the issues Copilot identified, you can request another Copilot review.

Pull request overview

This PR makes Underworld3’s “cell radius / mesh length” quantities partition-independent by switching the per-cell characteristic length to PETSc’s DMPlexComputeGeometryFVM (volume**(1/dim)), and wiring the public radius accessors and dependent defaults/tests to that definition.

Changes:

Replace the legacy rank-local kd-tree-derived radii with PETSc-derived per-cell volume**(1/dim) radii and use these for cell_size() and get_{min,max,mean}_radius().
Recalibrate dependent defaults (Nitsche gamma, max_surface_cfl) and update tests to assert partition independence (including a new parallel regression for the accessors).
Update downstream call sites / documentation to prefer global accessor APIs (get_mean_radius) over rank-local arrays.

File summaries

File	Description
tests/test_1065_nitsche_local_h.py	Updates local-h Nitsche regression and aligns expectations with new radius definition / gamma.
tests/test_1060_nitsche_freeslip.py	Switches Nitsche test to rely on default gamma and updates assertion comment.
tests/test_0010_cell_size_geometry.py	Rewrites the cell-size geometry oracle to validate `volume**(1/dim)` and deformation scaling.
tests/parallel/test_1078_radius_accessors_partition_independence.py	Adds parallel regression ensuring radius accessors are rank-count independent.
tests/parallel/test_1069_boundary_normal_parallel.py	Removes `local_h=False` workaround now that `cell_size()` is partition independent.
tests/parallel/test_0774_empty_rank_reductions_mpi.py	Updates empty-rank reduction negative control to use `_cell_radii`.
src/underworld3/systems/solvers.py	Ensures timestep estimation uses `_cell_radii` terminology/field consistently.
src/underworld3/systems/free_surface.py	Adjusts default `max_surface_cfl` to preserve prior effective surface step size.
src/underworld3/swarm.py	Updates comments to reference `_cell_radii` consistency for nondimensional velocity.
src/underworld3/meshing/smoothing/api.py	Updates docstring example to use `mesh.get_mean_radius()` instead of rank-local mean.
src/underworld3/discretisation/discretisation_mesh.py	Implements PETSc-based `_get_cell_radii()`, rewires rebuild paths, and updates accessors / docs.
src/underworld3/cython/petsc_generic_snes_solvers.pyx	Updates default Nitsche `gamma` and documents calibration against the new `h` definition.

Review details

Files reviewed: 12/12 changed files
Comments generated: 3
Review effort level: Lite

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

+    cell_start, cell_end = mesh.dm.getHeightStratum(0)
+    point_start, _point_end = mesh.dm.getDepthStratum(0)
+
+    volumes = np.empty(cell_end - cell_start)
+    factorial = 2.0 if dim == 2 else 6.0
+    for cell in range(cell_end - cell_start):
+        corners = mesh.dm.getTransitiveClosure(cell)[0][-(dim + 1):]
+        coords = mesh._coords[corners - point_start]
+        edges = coords[1:] - coords[0]
+        volumes[cell] = abs(np.linalg.det(edges)) / factorial
+    return volumes


        ## 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


        ## 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


lmoresi · 2026-09-06T20:01:01Z

Adversarial review of this PR

Copilot cannot be requested here (not a repo collaborator), so this is ours,
by execution.

What was attacked

The change reads radii from PETSc's cellgeom vector while centroids still
come from _get_coords_for_basis(0, False). Everything downstream indexes both
by cell, so the risks are ordering, starved ranks, and deformation.

check	result
radii aligned with the cell ordering everything else uses	verified — `test_0010` compares each cell's value against a determinant volume from that cell's vertices, max error < 1e-12
simplex `deform()` refreshes the radii	verified — doubling x scales lengths by sqrt(2), max error 5.6e-17
starved rank, build + accessors	verified at np=8 with cells per rank `[2,2,2,2,0,2,2,2]`
accessors agree across ranks and rank counts	verified np=1/2/4/8
both parallel guards	3 passed at np=2, np=4 and np=8
`test_1069` on the restored default `local_h=True`	7 passed at np=2 and np=8

One thing found, and it is NOT this PR's

At np=8 on a deliberately tiny mesh (14 cells, so one rank owns none),
mesh.cell_size() diverges: ranks 0-3 and 5-7 sit inside _assemble_cell_size
while rank 4 has already moved on to the next statement.

It is not introduced here:

_assemble_cell_size is documented as deliberately rank-local, and its
empty-partition early return dates from Nitsche BC: scale penalty by a local per-cell mesh size #275 (a7f0b11f), well before this
work or Make mesh.cell_size() local and partition independent #692.
The same 14-cell mesh at np=8 hangs on development earlier still —
during mesh construction, before cell_size() is reached — so that
configuration is pathological independently of anything here. That also means
a clean before/after attribution is not available at that size.

At realistic sizes np=8 is green throughout, as the table above shows. Worth its
own issue in the empty-rank family alongside #691, rather than holding this.

Corrections to my own earlier claims in this thread

I proposed the mean needed an order-independent reduction. It does not —
get_mean_radius() allreduces both the sum and the count, which is a true
global mean, and partition-order last-bit differences are out of scope
(maintainer, 2026-09-05).
I said this "needs the Nitsche and SUPG benchmarks re-run". The guards are
oracles, not baselines: test_1065 and test_1060 compare against an exact
essential free-slip solve computed in-test, test_1064 against published
Zhong values. Nothing stored a cell-size-derived number.
I first reported estimate_dt() as partition-dependent at ~1e-5. It is not —
re-running with the velocity written in analytically rather than solved gives
bit-identical answers at np=1/2/4/8. That spread was the iterative solve.

Residual risk, stated

The two constants that moved — gamma 10.0 -> 12.5 and max_surface_cfl
0.5 -> 0.437 — restore the measured behaviour of the tests that caught them.
They are not derived from first principles, and no benchmark outside the test
suite has been run. SUPG's guards and the adaptivity/relaxation consumers of
get_max_radius()/get_mean_radius() have not been audited.

…#698) _assemble_cell_size short-circuited on `radii.size == 0` and returned without ever touching `var.data`. That access is not rank-local: it lazily reaches MeshVariable._set_vec, which calls dm.createSubDM and createGlobalVector, both collective on the DM. So a rank owning no cells skipped two collectives its populated peers made, and the job diverged -- observed at np=8 with cells per rank [2,2,2,2,0,2,2,2]: ranks 0-3 and 5-7 inside _assemble_cell_size, rank 4 already past it. `var.data` is now read on every rank before the branch. The docstring claimed the whole routine was "purely RANK-LOCAL ... no collective", which was the false premise behind the early return; it now says which access is collective and why it has to come first. A zero-cell rank is routine on a region SUBMESH, which keeps only the cells carrying a label -- a partition whose share of the parent lies outside that region owns none. That is the case this matters for. HONESTLY REPORTED: this fix is argued from the code, not demonstrated by a reproducer. Neither route to a zero-cell rank could be made to run: * a 14-cell mesh at np=8 hangs during mesh CONSTRUCTION, before cell_size() is reached, on this branch and on development alike (intermittently -- one earlier run did get through and showed the divergence above); * `extract_region` on an annulus hangs at np=2 before returning, so the submesh route never reaches cell_size() either. Filed separately; it is not this. What is verified is that the change regresses nothing: parallel guards 3 passed at np=2 and np=4, serial cell-size and Nitsche tests 14 passed. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01E87Q7KrpapxeQiLD1RiNXv

… are collective (#698) The previous commit hoisted `var.data` above the empty-partition test, which fixed one collective and left the other. With a real reproducer the job still hung, one step later, with the populated ranks in `pack_raw_data_to_petsc`: assigning into `var.data` fires the array's write-back callback, and that is collective too. There is now no early return. A rank owning no cells walks the same path and writes a zero-length slice, so both collectives happen everywhere: data = var.data # createSubDM + createGlobalVector n = min(data.shape[0], radii.shape[0]) data[:n, 0] = radii[:n] # pack_raw_data_to_petsc; n == 0 is fine The reproducer is a region SUBMESH, which is the case that matters and the one the maintainer pointed at: a submesh keeps only the cells carrying its label, so a partition whose share of the parent lies outside that region owns none. A box split at z=0.85 gives an outer slab with cells per rank [12, 11, 0, 19, 0, 0, 0, 0] at np=8 -- five starved ranks out of eight, which is routine rather than contrived. My earlier attempts used an over-decomposed 14-cell full mesh instead. That was a poor choice: it is degenerate, it hits an unrelated intermittent hang in mesh construction, and it obscured the defect rather than exposing it. The previous commit's "argued from the code, not demonstrated" caveat is now discharged. tests/parallel/test_0779 guards it. The test asserts that every rank RETURNS -- a regression hangs rather than fails -- plus a second case that the field is still filled where cells exist, so the fix cannot be satisfied by never filling it. The starved count is reported and not asserted: how the parent splits is the partitioner's business, and at np=2 the slab reaches both ranks. np=8: cells [12, 11, 0, 19, 0, 0, 0, 0], starved 5, 2 passed. np=2: 2 passed. No regressions: parallel guards 3 passed at np=4, serial cell-size and Nitsche 14 passed. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01E87Q7KrpapxeQiLD1RiNXv

lmoresi and others added 3 commits September 5, 2026 22:24

Copilot AI lite review requested due to automatic review settings September 6, 2026 19:40

Copilot started reviewing on behalf of lmoresi September 6, 2026 19:40 View session

Copilot AI reviewed Sep 6, 2026

View reviewed changes

lmoresi mentioned this pull request Sep 6, 2026

cell_size() diverges on a rank that owns no cells (np=8, tiny mesh) — and that mesh already hangs in construction on development #698

Open

lmoresi and others added 2 commits September 6, 2026 13:26

Navigation Menu

Uh oh!

Take the cell radius from PETSc, and fix the radius accessors with it (#694) - #697

Take the cell radius from PETSc, and fix the radius accessors with it (#694)#697
lmoresi wants to merge 5 commits into
developmentfrom
bugfix/cell-radii-from-petsc-fvm

lmoresi commented Sep 6, 2026

Uh oh!

Copilot AI left a comment

Uh oh!

lmoresi commented Sep 6, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

Uh oh!

Conversation

lmoresi commented Sep 6, 2026

Summary

Why PETSc rather than the own-vertex RMS

The abandoned routine works

Two constants move with the length

Tests

Also

Uh oh!

Copilot AI left a comment

Choose a reason for hiding this comment

🟡 Changes recommended

Uh oh!

lmoresi commented Sep 6, 2026

Adversarial review of this PR

What was attacked

One thing found, and it is NOT this PR's

Corrections to my own earlier claims in this thread

Residual risk, stated

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants