Complete the boundary-normal sum across ranks, and make the partition-independence tests self-referential (#564)

[lmoresi]

Fixes the constrained half of #564: a boundary normal built from one rank's facets

Two independent things are in here, and they should be read separately.

1. A solver fix.mesh.boundary_normal() — the default constraint direction for
   add_constraint_bc and add_nitsche_bc — was assembled from each rank's own
   boundary facets. On a curved boundary that puts up to 3.3 degrees of error into the
   constraint direction at partition-seam nodes and moves a constrained free-slip
   answer by 3.4 % between np=1 and np=2. Fixed by completing the sum across ranks.
2. A test fix. The four remaining Constrained free-slip is partition-dependent (3.4%): boundary normals accumulated rank-locally — and the rotated tests compare against a golden recorded on another host #564 xfails are not a solver defect at all. The
   rotated path is partition-independent to the 13th significant figure; the tests were
   comparing against constants recorded on a different host. They are now
   self-referential.

---

1. The solver fix

The mechanism

The nodal boundary normal is Σ_f |f| n̂_f over every facet of the boundary that meets
the node, normalised at the end. A boundary facet is labelled on exactly one rank,
so a node on a partition seam sees only some of its facets locally. Normalising that
partial sum does not give a slightly noisy normal — it gives the wrong one.

Mesh._assemble_boundary_normal accumulated rank-locally, under a TODO(parallel)
that prescribed exactly the missing reduction and then asserted the code was
"parallel-SAFE as-is (only the handful of partition-seam surface vertices get a
slightly-rotated unit normal)". That sentence was the bug. It is replaced by the
measurement.

This is the same defect class as #560 (fixed for rotated_bc._boundary_velocity_nodes
in #561), in a different copy of the same accumulation. There were three copies.

Measured, before

Worst nodal normal against the exact radial one, Annulus(cellSize=0.12):

boundary	np=1	np=2	np=3	np=4
Upper (54 nodes)	3.0e-10	5.8e-02	5.8e-02	5.8e-02
Lower (28 nodes)	5.5e-10	1.1e-01	1.1e-01	1.1e-01

5.8e-02 is 3.3 degrees. It does not shrink with more ranks because it is not round-off:
it is what you get by taking ONE of a vertex's two facet normals instead of the average
of both, so its size is set by the facet's angular span. (54 nodes over 2π is 0.116 rad
per facet; half of that is 0.058.)

Measured, after

Identical to serial to every digit at every rank count — 3.005063e-10 (Upper) and
5.488346e-10 (Lower) at np=1, 2, 3 and 4. Serial results are bit-identical to before.

Same for a 3-D SphericalShell: bit-identical at np=1, 2, 4. (Its absolute error
against the analytic radial normal is 8.8e-02 on the outer boundary at cellSize=0.35
— that is honest faceting on a coarse triangulated sphere, not a defect, which is why
the 3-D test uses a global-facet-sum oracle instead of the analytic one.)

End to end — test_1063, run at four rank counts

Velocity L2 of the constrained free-slip annulus, this branch:

np	[iso]	[ti]
1	6.194547793939e-01	3.925981603048e-01
2	6.194547793947e-01	3.925981604290e-01
3	6.194547793476e-01	3.925981603382e-01
4	6.194547793462e-01	3.925981604662e-01

Spread 1.2e-10 and 4.1e-10.

With only the cross-rank sum disabled — a bisection, not an inference:

np	[iso]	[ti]
1	6.194547793939e-01	3.925981603048e-01
2	5.982807168537e-01 (−3.4 %)	3.939330452800e-01 (+0.34 %)
4	6.187878061677e-01 (−0.11 %)	4.036468804580e-01 (+2.8 %)

5.982807168537e-01 is #495's number to every digit. The mean-stripped topography
moved 2.4 % at np=2 on the same runs.

The 3.4 % closes completely. There is no residual, so there is no second mechanism
behind that row.

The fix

• The weighted contributions are summed through the variable's own sub-DM
  local↔global scatter before normalising (ADD into the global vector, scatter back),
  so every rank ends with a bit-identical normal. Work vectors are created from the
  sub-DM rather than borrowed from the variable — the TODO records that a first attempt
  SEGV'd on the lazily-built global vec, and that is what it hit.
• Collective: a rank owning no facet of the boundary still takes part (Empty-rank support, layer 2: evaluate/points_in_domain, radii/centroid reductions, gather_data NaN-stripping #405).
• No de-duplication needed, and this was verified rather than assumed: no boundary
  facet is labelled on two ranks and none is labelled away from its owner. Measured at
  np=2, 3, 4 on the annulus and the spherical shell, for both label sources — the
  per-boundary label this routine uses and the consolidated UW_Boundaries label.
• Orientation comes from the facet's own support cell with the getSupportSize == 1
  guard. An internal boundary has two support cells and support[0] is arbitrary;
  those facets are skipped (so boundary_normal("Internal") returns zero — a
  pre-existing limitation, now documented rather than silent).
• The DOF a facet contributes to now comes from the variable's own section on that
  sub-DM instead of a kd-tree lookup of the DOFs nearest the facet centroid. The
  kd-tree is a heuristic that can mis-assign on a graded mesh, and — decisively — it
  cannot be made to agree across ranks, because each rank's tree is built from its own
  local coordinates. The cross-rank sum would have landed on different points on
  different ranks.

Yes, a shared helper — one of three, not all five

utilities/facet_normals.facet_measure_and_normal now owns the orientation rule, the
measure and the internal-boundary guard, and all three copies of the facet loop use it:
rotated_bc._boundary_velocity_nodes, Mesh._assemble_boundary_normal, and
boundary_flux._node_normals — the last of which was still on the pre-#560 rule
(orient against the mean of this rank's coordinates, which points into the domain on
a concave boundary) under a TODO(BUG) from #561. Three copies drifting apart is
exactly how #564 happened twelve weeks after #560 fixed the first one.

The cross-rank reduction is deliberately not shared. The two are genuinely
different shapes: rotated_bc reduces a dict keyed by DMPlex point on the solver DM and
needs _local_boundary_candidates plus a constrained-DOF fallback; this one reduces a
dense array over the variable's local DOFs on a mesh sub-DM, where a node whose labelled
facets all live on a neighbour is simply a row that stays zero until the reduction fills
it, and where the mesh DM carries no essential-BC constraints. Forcing one function to
serve both would mean parameterising over the DM, the section, the field id, the node
enumeration and the fallback — more coupling than the ten lines it would save.

The other two flagged in #561 are left alone, with reasons:

• FreeSurface._normal_direction is unconditionally +radial. It is not a facet loop
  and shares nothing with this; it is a sign-convention question against the post-Rotated free-slip: the vertex normal is a bisector but assembly is facet-wise, so a kinked boundary with unequal facets leaves the pressure gauge unpinned #560
  σ_nn, and belongs with the free-surface work.
• fault_contact._fault_pair_nodes walks an internal surface where both support
  cells exist by design and the orientation is the ±side split. It is a different
  rule, not a stale copy of this one.

A second, unrelated defect found while testing — NOT fixed here

mesh.cell_size() is partition-dependent in its own right. _get_mesh_sizes
measures a cell by the distance from its vertices to the nearest centroid in a kd-tree
built from this rank's centroids, so near a seam the nearest centroid can simply be
absent. On Annulus(cellSize=0.12) the field's sum is 26.0822 / 26.1211 / 26.1386 at
np=1/2/4, and its max moves at np=4.

It reaches users through add_nitsche_bc(local_h=True), the default: the Nitsche
free-slip annulus moves 6.6e-03 in velocity between rank counts, and that does not
shrink with solver tolerance (stable from 1e-9 to 1e-12). With local_h=False the same
solve agrees to 3.6e-10 at np=1…4, which is how we know the normal is clean and the
local h is what is left.

Marked TODO(BUG) where it lives, with the numbers. Not fixed here: _get_mesh_sizes
also feeds get_min_radius, the adaptivity metrics and the free-surface relaxation, and
it needs its own benchmarking. It should get its own issue.

---

2. The test fix — the rotated rows of #564 are a TEST-DESIGN problem

This is a separate matter and this PR does not claim to fix a rotated solver defect,
because there is not one.

Running test_1064's own _annulus_diagnostics() on one CI host at both rank
counts: np=1 1.897329151623790e-02, np=2 1.897329151623740e-02 — agreement to the
13th significant figure, with leak_lo and leak_up identical to every digit. Both
differ from the recorded GOLDEN_ANNULUS by exactly the same +1.676e-04.

The whole discrepancy is that the goldens are constants recorded on macOS/arm64 while
gmsh builds a different triangulation on the Linux runner, and the files are
mpi(min_size=2) so CI has never run np=1 to notice. The assertion message said
"differs serial vs np=N" when what it measured was "differs from a number recorded
elsewhere" — and that phrasing is what read as seven instances of one solver bug and
sent the investigation the wrong way for a day.

Every partition assertion in test_1063, test_1064 and test_1066 now compares
against a np=1 run of the same file in the same environment
(tests/parallel/serial_reference.py spawns a single-rank child, scrubbing the
launcher's OMPI_*/PMIX_* so it runs as a singleton instead of trying to join the
parent job). Whatever mesh the host generates, both sides of the comparison use it.
Each failure message carries both runs' mesh fingerprints — owned cell count and
∫1 dV, both partition-independent — so a future host difference is legible.

Checks that are genuinely absolute are kept and relabelled as accuracy checks, with
a comment saying they are not partition checks: the analytic SolCx velocity error,
σ_nn's relL2 against the exact solution, the gathered top-node count, the Schur
iteration ceiling, the converged reason, and the datum residual. Conflating the two is
what this change exists to stop.

One tolerance moved on evidence: test_1063's topography goes 1e-6 → 1e-5. It is read
off the multiplier, whose [p,λ] Schur sub-block grinds into its 200-iteration cap on
this problem, and the [ti] case sits at 6.4e-07 against the old gate — a 1.6× margin
is not a tolerance. 1e-5 still leaves three orders below the 2.4 % it has to catch.

xfails removed — all seven

test	was
test_1063::test_constrained_freeslip_partition_independent[iso]	1.4 %
test_1063::…[ti]	0.3 %
test_1063::test_constrained_raw_gauge_partition_independent	2.3e-05
test_1064::test_rotated_freeslip_annulus_partition_independent	1.7e-04
test_1064::test_rotated_freeslip_spherical3d_partition_independent	1.1e-03
test_1064::…spherical3d_topography_partition_independent	5.9e-03
test_1066::test_rotated_datum_prescribed_normal_partition_independent	9.4e-06

The first three close because of the fix. The last four close because they now measure
what they claim to. grep -c xfail tests/parallel/test_106*.py is 0 across the board.

---

Tests

New: tests/parallel/test_1069_boundary_normal_parallel.py (7 tests), and
tests/parallel/serial_reference.py as the shared np=1-reference helper.

The new tests carry their own oracles rather than a recorded constant:

• analytic oracle — on a circular arc the measure-weighted average of a vertex's two
  chord normals is exactly radial by symmetry, so max‖n − ±r̂‖ on the annulus needs
  no golden and no serial run. Gate 1e-6: four orders below the defect, four above the
  discretisation floor;
• global-facet-sum oracle, 2-D and 3-D — the assembled field against the sum over
  the GLOBAL facet set, gathered and done in numpy. This never touches the section, the
  scatter or the DM after the gather, so it cannot reproduce a plumbing bug in the
  reduction. It is what covers the 3-D shell, where faceting error is honest;
• corner preservation on a box — Top keeps (0,1) and Right keeps (1,0) at the shared
  vertex. This is the regression the fix could most easily introduce;
• add_nitsche_bc end to end on a curved boundary with the default normal — the
  other exposed consumer, previously untested. Pre-fix it moved 6.6e-03 in velocity and
  more than doubled the wall-normal leakage across rank counts.

Negative controls, both run. With _sum_local_dofs_across_ranks reverted to a no-op
(i.e. the pre-fix code), 7 of 10 assertions across test_1063 + test_1069 fail at
np=2 — the three that pass are the normal=unit_r gauge test (correctly unaffected, and
the same bisection the investigation used), the box corner test (flat walls are immune)
and the control itself. With a 1e-3 relative partition dependence injected into every
diagnostic, all 14 self-referential assertions across the four files fire. Neither
suite passes vacuously.

Runs

what	result
test_1069 new	7 passed at np=2, 7 at np=4
test_1063	3 passed at np=2, 3 at np=4 (was 3 xfail)
test_1064 + test_1066	11 passed at np=2, 11 at np=4 (was 4 xfail)
parallel batch tests/parallel/test_10*py	34 passed at np=2 (120 s), 34 passed at np=4 (120 s), 0 xfail
serial stakeholders — test_1018_rotated_freeslip, test_0056_projected_normals_deform, test_1060_nitsche_freeslip, test_1061_constrained_freeslip, test_1062_constrained_solcx, test_1065_nitsche_local_h, test_1065_rotation_gauge_freeslip, test_1024_multiplier_schur_pc	56 passed, 0 failed
scripts/test.sh --p 2	see the note below

tests/test_1064_constrained_spherical_shell_response.py is the one stakeholder not
run to completion: it is level_3 / slow / tier_c, is deliberately excluded from
scripts/test.sh's serial batches, and its 3-D LU solves were still running after 35
minutes. It exercises add_constraint_bc on a shell, so it is worth a run before merge.

---

Also

• docs/developer/subsystems/rotated-freeslip.md gains a section on the other two
  free-slip paths, the shared helper, the internal-boundary zero and the cell_size
  caveat, so the doc that already explains this rule for the rotated path now says where
  else it applies.
• The TODO(parallel) at discretisation_mesh.py is gone, replaced by the measurement,
  as is the TODO(BUG) on boundary_flux._node_normals' orientation (its remaining gap
  — no cross-rank sum on an unreachable branch — is stated precisely instead).

Commits

sha	subject
4d660447	A boundary normal built from one rank's facets is a rotated normal: complete the sum across ranks
91384dd3	The constrained partition tests now agree at every rank count: drop their xfails
5575f886	The rotated partition tests were comparing against another host: make them self-referential
a4dc58c0	Keep the new collective out of reach of an existing swallow-and-continue

The last one is worth reading on its own: completing the sum made
_assemble_boundary_normalcollective, and deform() already called it inside a
try: … except: pass. A rank-local failure there would have taken one rank out of the
reduction and hung the rest — a swallowed error turned into a deadlock. The rank-local
guard now sits inside the collective, so the outer one can only ever see a symmetric
failure.

Fixes the constrained half of #564, and #495. The rotated rows of #564 are addressed as
the test-design problem the CI experiment showed them to be; no rotated solver change is
claimed.

Underworld development team with AI support from Claude Code

[Copilot]

Pull request overview

This PR fixes a parallel correctness issue in mesh.boundary_normal() by completing the measure-weighted facet-normal accumulation across MPI ranks before normalization, eliminating partition-dependent constraint directions on curved boundaries. It also updates several parallel solver tests to be self-referential by computing their own np=1 reference in the same environment, avoiding false failures caused by host-dependent meshing differences.

Changes:

• Complete the per-node boundary-facet normal sum across ranks (collective reduction via the variable’s sub-DM local↔global scatter) before normalization.
• Centralize facet measure + outward-orientation logic in a shared helper (facet_measure_and_normal) and reuse it across the three facet-walking accumulators to prevent drift.
• Replace hardcoded “golden” serial constants in parallel partition-independence tests with an in-environment np=1 child-run reference (serial_reference.py), and add a dedicated parallel regression suite for mesh.boundary_normal().

Reviewed changes

Copilot reviewed 11 out of 11 changed files in this pull request and generated no comments.

Show a summary per file

File	Description
src/underworld3/discretisation/discretisation_mesh.py	Fix boundary_normal() assembly by summing facet contributions across ranks before normalizing; add collective reduction helper.
src/underworld3/utilities/facet_normals.py	New shared utility defining a single outward-orientation + measure rule for boundary facets.
src/underworld3/utilities/rotated_bc.py	Switch rotated free-slip facet loop to use the shared facet normal/measure helper.
src/underworld3/utilities/boundary_flux.py	Align geometric normal accumulation with the shared facet rule; document remaining parallel reduction TODO for an unreachable branch.
tests/parallel/serial_reference.py	New test helper to compute per-test-module np=1 diagnostics via a scrubbed-environment singleton child process (self-referential references).
tests/parallel/test_1063_constrained_freeslip_parallel.py	Make partition-independence assertions compare against in-environment np=1 references rather than stored constants; adjust tolerances with documented rationale.
tests/parallel/test_1064_rotated_freeslip_parallel.py	Convert rotated partition-independence checks to use serial_reference and keep absolute “accuracy” gates explicitly separate.
tests/parallel/test_1066_rotated_datum_parallel.py	Remove xfail + hardcoded serial constant; compare solve-energy against in-environment np=1 reference while keeping an absolute datum-imposition gate.
tests/parallel/test_1069_boundary_normal_parallel.py	New parallel regression tests for mesh.boundary_normal() including analytic and global-facet-sum oracles plus a negative control and Nitsche end-to-end coverage.
scripts/test.sh	Enable/clarify execution of the tests/parallel/test_10*py solver-parallel suite in CI.
docs/developer/subsystems/rotated-freeslip.md	Document that the same facet-normal rule applies to constrained/Nitsche paths and note remaining caveats (internal boundaries, cell_size() partition dependence).

---

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

[lmoresi]

Adversarial review — PR #568

"Complete the boundary-normal sum across ranks, and make the partition-independence
tests self-referential", branch bugfix/issue-564-boundary-normal,
commits 4d660447 / 91384dd3 / 5575f886 / a4dc58c0, base bdb56713 (#561).

Reviewed independently of the implementer. Worktree r568-review, own amr-dev
env, built from pr568-head (uw.__file__ verified in the worktree's
site-packages). All runs sequential.

---

Verdict

Request changes. The diagnosis is right, the mechanism is right, the fix is right,
and the evidence behind it is unusually good — the no-de-duplication claim is guarded
by a test rather than asserted, and the negative control fires. We reproduced every
headline number and extended them to np=8. Three things block:

1. mesh.boundary_normal() changes by up to 0.089 in the unit normal (≈5.1°) in
   SERIAL in 3-D, and the PR says serial is unchanged. The new value is the correct
   one — the old one was 0.103 (≈5.9°) from the true facet sum, so this is a second,
   previously unknown defect that the PR fixes by accident — but it is undisclosed,
   undocumented, and untested (no test in the tree reaches the default normal on a
   3-D curved boundary at np=1).
2. The new except Exception: accum[...] = 0.0 inside _assemble_boundary_normal
   silently reproduces Constrained free-slip is partition-dependent (3.4%): boundary normals accumulated rank-locally — and the rotated tests compare against a golden recorded on another host #564 on the error path. Measured: one rank failing the facet
   walk leaves that rank's owned boundary DOFs with a zero normal, no message,
   converged solve.
3. deform()'s guard is not collective on every path. Measured: a rank-local
   failure raised before dm.createSubDM hangs the job (mpirun --timeout 180, exit
   241). The PR's own comment claims the outer guard "only ever sees a symmetric
   failure"; that holds only for the subset of exceptions raised inside the inner try.

None of these is a design objection. (1) is disclosure plus one serial test, (2) is a
collective flag, (3) is moving one line inside the existing try.

---

MERGE-BLOCKERS

B1 — boundary_normal() moves by 0.10 in serial in 3-D, and the PR says it does not

The PR body: "Serial results are bit-identical to before. Same for a 3-D
SphericalShell: bit-identical at np=1, 2, 4." We re-implemented the pre-PR kd-tree
assembly verbatim from bdb56713 and ran both, np=1, in one process:

mesh / boundary	max‖n_new − n_old‖	new vs global-facet-sum oracle	old vs oracle
Annulus(0.12) Upper	1.11e-16	1.24e-16	1.24e-16
Annulus(0.12) Lower	1.11e-16	0.00	1.57e-16
StructuredQuadBox(8,8) Top / Right	0 (bit-identical)	—	—
SphericalShell(0.55,1.0,cs=0.35) Upper	4.32e-02	1.92e-16	4.71e-02
SphericalShell Lower	8.94e-02	2.23e-16	1.03e-01

The oracle is the same construction test_1069 uses, run in serial where there is no
partition effect at all, so the only variable is the DOF-row change. Read it plainly:
the old kd-tree assignment was wrong by up to 0.103 in serial in 3-D and the new
section-based one is exact. On a tetrahedral boundary the three DOFs nearest a face
centroid are not always that face's three vertices, and k=nverts silently grabs a
neighbour. The PR argues the kd-tree "can mis-assign on a graded mesh" — it mis-assigns
on a uniform one, and the PR does not appear to know it.

Why this blocks rather than being a bonus:

• It is a numerical change to mesh.boundary_normal(), i.e. to the default constraint
  direction of add_constraint_bc and add_nitsche_bc, and to
  systems/free_surface.py:635,645 (which the PR's stakeholder list omits) — in
  serial, for every 3-D curved boundary. Any existing 3-D result with normal=None
  moves. The repo's own standing instruction is that solver-facing numerics do not
  change silently.
• The radial oracle cannot see it: at this resolution the faceting error is
  8.8e-02 (Upper) / 2.26e-01 (Lower) and the worst node is the same for old and new, so
  max‖n − r̂‖ is identical before and after. Only the facet-sum oracle catches it.
  This is exactly why the PR's own 3-D justification ("that is honest faceting … which
  is why the 3-D test uses a global-facet-sum oracle") is right and why the number
  needed to be measured rather than inferred.
• No test covers it.test_1069 is pytest.mark.mpi(min_size=2), so it never runs
  at np=1. Sweeping the tree: every serial add_nitsche_bc / add_constraint_bc call
  without normal= is on a box (test_1060, test_1065_nitsche_local_h,
  test_0641, test_1066) — flat walls, immune. test_1064_constrained_spherical_shell_response.py
  passes explicit normal=±unit_r on all four BCs (lines 104-117). So the 3-D serial
  path had, and still has, zero coverage.

Asked for: correct the claim in the PR body and the docs (it is a second fix, and
a good one — say so); and cover the 3-D serial path. The remedy is one line: we deleted
pytest.mark.mpi(min_size=2) from test_1069's pytestmark and ran it at np=1 —

6 passed, 1 skipped in 8.68s
SKIPPED tests/parallel/test_1069_boundary_normal_parallel.py:213:
a rank-local stencil is the complete stencil in serial

— every oracle, the corner test and the Nitsche end-to-end pass serially in 8.7 s, and
the only test that cannot run at np=1 skips itself, by a guard the file already carries.
Had that mark not been there, the 3-D facet-sum oracle would have caught the kd-tree
defect on its own. (The file is under tests/parallel/, so scripts/test.sh still only
reaches it via the --p N line; a serial companion under tests/ would be needed for
the serial job proper.)

B2 — the new swallow reproduces #564, silently, as its error path

discretisation_mesh.py, inside _assemble_boundary_normal:

exceptException:
accum[...] =0.0

Keeping the guard inside the collective is the right call and the reasoning in the
comment is correct. Making it silent is not. The except is unqualified and wraps the
whole facet walk including the new indexing (ssec.getOffset(q) // ncomp), so an
IndexError from the very code this PR introduces is caught and turned into "this rank
contributed nothing".

Measured (np=2, Annulus(cs=0.2), facet_measure_and_normal patched to raise on rank 1
only, then deform()):

| | Σ|n| over local rows, rank 0 | rank 1 |
|---|---:|---:|
| control | 21.306341 | 21.306341 |
| rank-1 facet walk raises | 21.492745 | 2.186404 |

deform() returns on both ranks, nothing is printed, and rank 1's owned boundary DOFs
come back with a zero normal — so the Nitsche/constraint direction is the zero
vector over that part of the boundary. That is strictly worse than the 3.3° the PR
exists to remove, and it is indistinguishable from success.

Asked for: set a local failed flag in the except, comm.allreduce(…, MAX) it
alongside the reduction that is already collective, and raise (or at minimum
uw.pprint a warning) on every rank when any rank failed. The collective is already
there; this costs one extra reduction on the same communicator.

B3 — the deform() guard is not collective on every path

The PR moved the guard inside the collective and wrote:

"this outer guard therefore only ever sees a symmetric failure"

That is true only for exceptions raised inside the inner try. Everything before
dm.createSubDM(var.field_id) — var.num_components, numpy.zeros_like(numpy.asarray(var.data))
— and the one statement between createSubDM and the inner try,
ssec = subdm.getLocalSection(), is outside it, and createSubDM / localToGlobal
are collective.

Measured, np=2, _assemble_boundary_normal patched to raise on rank 1 only, under
mpirun --timeout 180:

[rank 1] entering deform (MODE=outer)
[rank 0] entering deform (MODE=outer)
[rank 1] deform returned
<rank 0 never returns; job killed by the launcher timeout, exit 241>

Rank 1 walks out through deform()'s except Exception: pass and rank 0 blocks in the
sub-DM collective. Before this PR the routine was rank-local and this was benign; the
PR created the hazard and closed part of it. deform() runs every free-surface
timestep, so a hang here is a wall-clock loss on a cluster, not a test failure.

Asked for: move ssec = subdm.getLocalSection() inside the inner try (one line —
it closes the only realistic remaining hole), and replace deform()'s per-boundary
except Exception: pass with a collective decision: allreduce the failure flag from B2
so every rank takes the same branch and an asymmetric failure aborts loudly instead of
hanging.

---

Findings — high

H1 — test_1063's gates were loosened 10× and the PR body does not say so

assertion	before	after
test_constrained_freeslip_* velocity	rtol=1e-9	rtol=1e-8
test_constrained_freeslip_* topography	rtol=1e-6	rtol=1e-5
test_constrained_raw_gauge_* velocity	rtol=1e-8	rtol=1e-8 (unchanged)
test_constrained_raw_gauge_* topography	rtol=1e-6	rtol=1e-5

The PR body reports the residual spread as "1.2e-10 and 4.1e-10" and concludes "the
3.4 % closes completely". Those are absolute; the gates are relative. 4.1e-10 on a
[ti] velocity of 0.3926 is 1.0e-9 relative — i.e. the [ti] case sits exactly on
the old 1e-9 gate, and the loosening was necessary, not cosmetic. That is a defensible
engineering call (one order of headroom on a solve whose [p,λ] Schur block caps at 200
iterations), but a reader of the PR body will assume the old assertions now pass, and
they do not. The topography loosening at least carries its reason in a comment.

Asked for: one line in the PR body stating both loosenings and the 1.0e-9 relative
[ti] spread that forces the velocity one.

H2 — eleven absolute anchors deleted; the design note's own recommendation was half-implemented

The self-referential rewrite is the right diagnosis and the right shape. But
GOLDEN_BOX, GOLDEN_ANNULUS, GOLDEN_ANNULUS_FMG, GOLDEN_SPHERICAL3D,
GOLDEN_SPHERICAL3D_TOPO, GOLDEN_BOX_NONLINEAR, GOLDEN_BOX_SIGMA,
GOLDEN_TOPO_BDL2, GOLDEN, GOLDEN_GAUGE and _INT_VV_REF are all removed, and only
some are replaced by an accuracy check. Nothing now pins:

• the annulus radial leakage (4.6e-05 / 9.3e-06) — a rotated constraint that stopped
  constraining at every rank count equally would pass;
• the custom-FMG annulus answer — no absolute anchor at all, and it is not
  cross-checked against the non-FMG annulus either (1.907e-02 vs 1.897e-02);
• the Zhong l=2 topography coefficients — physics benchmark numbers, now free;
• test_1066's ∫v·v; test_1063's velocity and topography.

The investigation note that this PR implements says (§"Recommended fix shape", and again
in the test-plan recipe):

"Where a stored constant is unavoidable, store a mesh fingerprint beside it … and
assert the fingerprint first, so a mesh change reports itself instead of masquerading
as a physics regression."

The PR built mesh_fingerprint() — and then deleted the constants instead of gating
them. That is the more expensive half (a full serial solve per test, see M3) and the
less protective one. Keeping each golden at a loose rtol=1e-2 as an explicitly-labelled
accuracy gate, alongside the self-referential tight gate, costs nothing and restores
what was lost.

Related and cheap: compare()reports both fingerprints in the failure message but
never asserts they match. The docstring's promise ("the day that stops being true it
says so") is only kept if the test fails for another reason first. Assert it.

---

Findings — medium

M1 — serial_reference deadlocks if the child launch fails for anything but a timeout

_run_child catches subprocess.TimeoutExpired only. Any other exception on rank 0 —
OSError/PermissionError from subprocess.run, a json.JSONDecodeError on a
truncated SERIALREF line, a MemoryError on capture_output — propagates out of
serial_referencebeforeuw.mpi.comm.bcast, leaving every other rank blocked in
the broadcast. pytest-timeout eventually fires (900 s), and Open MPI busy-polls, so those
ranks burn a core each in the meantime. One try/except Exception returning the message
string (the function already has that contract for the timeout case) closes it.

Also: the child's default timeout=1800 is longer than the pytest.mark.timeout(900)
that wraps it, so it can never fire. Make it a fraction of the mark, or drop it.

M2 — CI cost: every partition test now runs a full serial solve in a child

test_1064 spawns eight children, including the 3-D spherical shell twice
(spherical3d, spherical3d_topo). Measured here at np=2 (warm gmsh cache):

file	np=2 wall	tests
test_1063	30 s	3
test_1066	10 s	2
test_1064	64 s	9
test_1069	9 s	7

The timeouts were raised 180 → 600/900 s to accommodate it, which is honest. Two things
to weigh: ranks 1..N-1 sit in MPI_Bcast while the child runs, and Open MPI busy-waits
by default, so on a 2-vCPU GitHub runner the child competes with a spinning rank for the
same core; and _CACHE is per-process, so the eight children are re-run for every fresh
pytest invocation. Neither blocks. If it bites, OMPI_MCA_mpi_yield_when_idle=1 in the
scrubbed child env, or an on-disk cache keyed by the mesh fingerprint, are the levers.

M3 — the dropped fallback in _sum_local_dofs_across_ranks tests the wrong predicate

dropped= (~summed.any(axis=1)) &values.any(axis=1)
summed[dropped] =values[dropped]

The question is "was this DOF constrained out of the global vector?" The code asks "did
it come back all-zero?". Those differ on precisely the pathological case: a node whose
global contributions cancel exactly (opposed facets on a degenerate or
zero-thickness boundary) has its rank-local partial value restored — reintroducing
the partition dependence on the one input where it matters. The author states it is dead
today, which also means it is untested and the negative control cannot reach it.
ssec.getConstraintDof(q) is the actual predicate; or delete the branch and let a future
constrained section fail loudly.

M4 — "a different rule, not a stale copy" overstates the fault_contact exclusion

The PR excludes fault_contact._fault_pair_nodes on the grounds that it "walks an
internal surface where both support cells exist by design and the orientation is
the ±side split. It is a different rule, not a stale copy of this one." Reading
fault_contact.py:561-571, the loop is the same six lines as the block just
extracted from rotated_bc — same computeCellGeometryFVM, same unit-normalise with
the identical +1e-30, same flip against getSupport(f)[0]'s centroid, same
float(vol) measure weight, same full-length dot — with exactly one difference: the
getSupportSize == 1 guard is dropped so the flip is unconditional.

That is a different guard, not a different rule, and it is a fourth copy of the six
lines whose duplication the PR itself calls out ("three copies drifting apart is exactly
how #564 happened twelve weeks after #560 fixed the first one"). facet_measure_and_normal
already returns exterior precisely so a caller can decide what to do about it, so
fault_contact can consume the helper and apply its own flip, keeping the
normalise/epsilon/measure conventions shared. Not a blocker — the code is correct today —
but the stated reason for leaving it out does not survive reading it.

M5 — _assemble_boundary_normal assumes one node per DMPlex point

accum[ssec.getOffset(q) // ncomp] is the row index only when every point carries
exactly ncomp DOFs. True for the degree-1 variable boundary_normal() builds — but
that path is existing = self.vars.get(f"_n_bd_{name}"), so a pre-existing variable of
another degree (checkpoint restore, user code) would be accepted and, for a P3 edge (two
nodes per point in 2-D), would accumulate onto the first node only and leave the second
at zero. One assert var.degree == 1 at the top, or take the row from
var.coords-independent section arithmetic that handles multi-node points as
boundary_flux._boundary_field_nodes already does.

---

Findings — low

• L1 — boundary_flux._node_normals is genuinely dead; the concave sign change does
  NOT propagate. Verified: the sole call site is boundary_flux.py:638, guarded by
  if normal is not None, and boundary_flux_field (:691) routes through the same
  boundary_flux(). With normal=None a vector solver takes the full-traction branch
  and never calls _node_normals. boundary_normal_traction / dynamic_topography
  live in rotated_bc and read the rotated constraint reaction, not this. So the
  Rotated free-slip: weight the nodal normal by the facet measure the assembly integrates over (#560) #561 review's concave σ_nn sign flip is not re-applied a second time. The PR's
  claim here is correct.
• L2 — facet_measure_and_normal dots full-length vectors in the flip test
  where the old _assemble_boundary_normal sliced to [:cdim]. Identical whenever
  mesh.cdim == dm.getCoordinateDim(); a 2-D manifold embedded in 3-D would differ.
  One line in the helper docstring.
• L3 — "bit-identical" is used for the 2-D serial result, which is actually
  1.11e-16. Immaterial numerically, but the PR leans on the word.
• L4 — the test files assert absolute gates before the collective
  serial_reference() call (test_1066 relL2, test_1064nnodes/its/reason).
  Every one of those quantities is collective or bcast today, so it is symmetric; but
  the pattern is a latent hang generator and nothing enforces the invariant. A comment
  at least.
• L5 — test_1069 is mpi(min_size=2), so the corner test — which the file itself
  calls "the regression the fix is most likely to introduce" — and the negative control
  and the analytic oracle never run in the serial job. See B1.
• L6 — the stakeholder list omits systems/free_surface.py:635,645, which reads
  mesh.boundary_normal(self.surface) when self.normal is None.
• L7 — _spherical3d_topography_diagnostics's cell_size argument is no longer
  reachable from __main__ (_DIAGNOSTICS[_kind]() takes no argv). Minor loss of a
  debugging affordance the old block had.

---

What we verified that holds up

Measured in this worktree, sequentially.

Parallel test batch — all green, including odd rank counts.

file	np=2	np=3	np=4
test_1063	3 passed	3 passed	3 passed
test_1064	9 passed	9 passed	9 passed
test_1066	2 passed	2 passed	2 passed
test_1069	7 passed	7 passed	7 passed

No xfails, no skips, no strict-xpass. The seven xfails are genuinely gone rather than
weakened: test_1064's and test_1066's rtols are unchanged from the pre-PR goldens
(H1 is confined to test_1063), the nonlinear iteration count is still exact equality
(rtol=0.0, atol=0), and the absolute checks that were doing real work (verr < 1e-3,
relL2 < 0.10, nnodes == 49, reason > 0, its <= 25, datum relL2 < 1e-3) are all
retained and correctly relabelled.

Empty ranks and np=8.StructuredQuadBox(elementRes=(1,32)) gives Top exactly one
facet, so at np=8 seven of eight ranks own no facet of the boundary being normalised.
No hang; max‖n − (0,1)‖ = 0 at np=1, 2, 4 and 8. The curved case is stable to the digit
at every rank count:

	np=1	np=2	np=4	np=8
annulus Upper max‖n − r̂‖	3.005e-10	3.005e-10	3.005e-10	3.005e-10
annulus Lower max‖n − (−r̂)‖	6.255e-10	6.255e-10	6.255e-10	6.255e-10

No facet is double-labelled — independently confirmed. Per-rank labelled-facet counts
sum to exactly 54 (Upper) and 28 (Lower) at np=1, 2, 4 and 8, matching the serial totals.
The PR's "no de-duplication needed" is not an assumption: test_1069's facet-sum oracle
de-duplicates by centroid while the assembly does not, so a facet labelled twice would
double its weight and break the 1e-12 gate. That is a properly load-bearing guard.

The negative control fires (test_boundary_normal_oracle_fires_without_the_reduction
monkeypatches _sum_local_dofs_across_ranks to the identity and both oracles must exceed
1e-3), and it correctly skips at np=1.

The new sub-DM coupling survives deform(). The design note flagged deformed/adapted
meshes as a risk, and the new code resolves DOF rows through
dm.createSubDM(var.field_id) on the current DM where the old code used var.coords
and never touched the DM's fields — a genuinely new coupling to the remesh. Radial
scaling by 1.37 on Annulus(cs=0.15), worst‖n − r̂‖:

	np=1	np=2	np=4
before deform() (Upper / Lower)	4.0503e-10 / 7.6439e-10	same	same
after deform()'s eager refresh	4.0503e-10 / 7.6439e-10	same	same
after an explicit re-request	4.0503e-10 / 7.6439e-10	same	same

Identical to the digit in all nine cells. The refresh path is sound.

The rotated_bc extraction is behaviour-preserving. Line-by-line, the new
facet_measure_and_normal reproduces the inlined block from bdb56713 exactly:
same computeCellGeometryFVM, same unit-normalise with the +1e-30, same
getSupportSize == 1 guard, same flip against the support-cell centroid, same
wgt = float(vol). This is a pure extraction. (See L2 for the one edge case.)

The gap the PR admits is not a gap.tests/test_1064_constrained_spherical_shell_response.py
passes explicit normal=unit_r / normal=-unit_r to every add_nitsche_bc and
add_constraint_bc (lines 104-117), so it never reaches mesh.boundary_normal() and
cannot be affected by this change — provable by reading, no 35-minute LU run needed.
It is also excluded from CI by construction: scripts/test.sh deliberately does not
batch test_106*py. Risk from not running it: nil. We recommend the PR say that
instead of listing it as unverified.

The serial stakeholders — and a note on where they live. The consumers of the
default (assembled) normal are test_1060_nitsche_freeslip,
test_1061_constrained_freeslip, test_1062_constrained_solcx,
test_1065_nitsche_local_h, test_1065_rotation_gauge_freeslip,
test_1066_stokes_jacobian_layout, plus test_1024_multiplier_schur_pc and
test_0056_projected_normals_deform. 41 passed, 0 failed in 123 s.

Worth flagging separately: most of those are test_106*py, which scripts/test.sh
deliberately does not batch (the comment at scripts/test.sh:82-84 excludes
test_106*py / test_107*py pending an issue-#504 triage). So a green scripts/test.sh
does not exercise the default-normal path in serial at all — they have to be run by
hand, as above. That is a pre-existing CI gap, not this PR's, but it compounds B1: the
one batch that would have noticed a 3-D serial normal change is the one CI does not run.

scripts/test.sh (serial). Run end to end from the review worktree. Every batch that
had completed at the time of writing is clean — zero FAILED, zero ERROR lines:

batch	result
test_000*…test_02*	129 + 43 + 168 + 27 passed
test_05* + test_07*	532 passed, 5 skipped, 9 xfailed, 1 xpassed
test_08*	603 passed, 13 skipped, 11 xfailed
test_101* + test_102*	143 passed
test_105* onward (VEP/VE, diffusion, advection, named)	still running at time of writing

The tail that had not finished (test_105*, test_11*, test_1450*, the named files)
contains no consumer of mesh.boundary_normal — we grepped the whole tree for
add_nitsche_bc / add_constraint_bc / boundary_normal( and every hit is in the
files listed above or in the stakeholder set run separately. We would still want the
tail green before merge, but it carries no exposure to this change.

---

Reproduction

./uw worktree create r568-review
cd .claude/worktrees/r568-review
git fetch origin pull/568/head:pr568-head && git reset --hard pr568-head
./uw build
export PATH="$PWD/.pixi/envs/amr-dev/bin:$PATH"
# parallel batch
for N in 2 3 4; do for F in tests/parallel/test_106{3,6,4,9}*.py; do
mpirun -n $N python -u -m pytest --with-mpi -q -p no:randomly $F; done; done
# probes (scratchpad)
python p1_serial_bitid.py # old kd-tree vs new, np=1
python p2_serial_oracle.py # both vs the facet-sum oracle, np=1
mpirun --oversubscribe -n 8 python p3_emptyrank.py
mpirun --timeout 180 -n 2 python p4_deform_collective.py {none,inner,outer}

Probe sources:
/private/tmp/claude-501/-Users-lmoresi--Underworld-underworld3-pixi/64e34dd9-4af2-4774-9d55-4238a998eb6e/scratchpad/p{1,2,3,4}_*.py

---

Underworld development team with AI support from Claude Code

[lmoresi]

Response commits bb3c2819, 9270aca4, a2c7a5cd, 79c348b2 — all three blockers, both high findings and the minors.

B1 — the serial 3-D defect is now the PR's headline, not an accident. Reproduced independently beside the pre-PR kd-tree assembly in one serial process, matching the review digit for digit:

np=1, vs the global facet sum	old	new
SphericalShell(0.55,1.0,cs=0.35) Upper	4.712169e-02	1.92e-16
SphericalShell Lower	1.029976e-01 (5.9°)	2.23e-16
Annulus(0.12) Upper / Lower	1.24e-16 / 1.57e-16	1.24e-16 / 0.0

mpi(min_size=2) deleted (test_1069: 7 passed, 1 skipped at np=1), a dedicated serial 3-D assertion added so it cannot silently become parallel-only again, and the file wired into scripts/test.sh's SERIAL job — the gap that let this hide. Documented in the boundary_normal docstring, the test module docstring and the PR body, which now leads with it: 2-D is unaffected; every 3-D normal=None curved-boundary result moves, toward the right answer.

B2 — the rank-local flag is all-reduced and every rank raises, carrying the failing ranks' messages. Control at np=2: a rank-1 facet-walk failure raises on both ranks and both reach the next barrier; Σ|n| never collapses (was 2.186 vs a 21.306 control, silently).

B3 — getLocalSection and the pre-createSubDM region are inside the guard, and deform() all-reduces on the exception path and warns rather than skipping silently. Control: a pre-createSubDM failure now raises on both ranks and both reach the barrier (was exit 241 on the launcher timeout). Stated rather than hidden: replacing the entire routine still hangs — rank 0 is inside createSubDM before any flag can be exchanged. What remains unguarded is createSubDM, two destroys, the raise-path allgather and the final pack — all collective, hence symmetric.

H1 — both loosenings stated with the measurement that forces them: 4.1e-10 absolute on a [ti] velocity of 0.3926 is 1.0e-9 relative, exactly on the old gate.

H2 — all eleven anchors restored behind accuracy_anchor(), gated on the mesh fingerprint (skip, don't fail, on a different mesh) at rtol=1e-2, alongside the tight self-referential gate. Every value is the pre-PR golden, each verified to reproduce. compare() now asserts the fingerprints match instead of only printing them. Partition independence and accuracy are separate claims and now have separate tests.

M4 — the reviewer was right: it was the same six lines with one guard dropped, so fault_contact now consumes the helper via orient="support0", killing the fourth copy (48 fault tests pass). M1 _run_child catches everything and returns the message (timeout 1800→600 s, shorter than the mark wrapping it). M3 branch deleted — it restored a rank-local partial value on exactly the input the reduction exists to correct. M5 refuses a non-degree-1 variable. M2 accepted: 114 s for the whole test_10*py batch, and the alternative — dropping the np=1 side — is what produced #564's four false rows.

Verification: np=1 test_1069 7 passed/1 skipped · np=2/3/4 tests/parallel/test_10*py34 passed, 1 skipped each · fault suite 48 · serial stakeholders 107 · free surface 8 · scripts/test.sh --p 2 exit 0, 1923 serial + 38 + 34 parallel, zero FAILED/ERROR, batch counts identical to the pre-review run.

The CI gap this exposed is filed separately as #570: four coverage gaps in two days, each hiding a real defect, and test_106* — where the serial default-normal consumers live — is the one that hid B1.

Underworld development team with AI support from Claude Code

[lmoresi]

Re-review addendum — APPROVE

All three blockers cleared on measurement.

B1 — the new serial 3-D guard genuinely discriminates. Restoring the pre-PR kd-tree routing and re-running the actual test functions at np=1: the dedicated serial 3-D shell test and the parametrised facet-sum [shell] case both FAIL at 4.712e-02, while [annulus] and the exact-radial test still pass. That is the right discrimination — 3-D catches it, 2-D does not move, and the analytic radial oracle stays blind, which is precisely why the facet-sum oracle exists. CI wiring confirmed live: the serial batch appears in the test.sh log as its own line (7 passed, 1 skipped) and the workflow runs ./scripts/test.sh --p 2.

B2 — the attack does not land. A rank-1-only failure injected at each reachable point (facet walk, the guarded set-up before createSubDM, degree-2, and through deform()) raises on both ranks every time, with no hang; through deform() the normals are left at their pre-deform value instead of the earlier silent-zero corruption. The reduce cannot be skipped — the only statement between the guarded set-up and the allreduce is createSubDM itself.

B3 — the disclosure is honest, not optimistic. Replacing the whole routine gives exit 241 with PETSc error 98 from the mismatched collective, and triggering it required replacing the bound method, which no caller can do. What remains unguarded is genuinely all collective. Recording that plainly so it does not resurface as a phantom blocker.

H2 — the anchors fire here (fingerprint matches, zero skips at np=2 and np=4), and a same-mesh regression fails rather than skips: +3%, +3.4% (the #564 velocity move) and +2.4% (its topography move) all raise, while +0.5% passes inside the deliberate gate. The fingerprint discriminates round-off (1e-14, 1e-13 → same) from a re-triangulation (1e-3 → different), and compare() precedes accuracy_anchor() at all 12 call sites (AST-checked), so a skip can never take the partition assertion with it.

Also confirmed bb3c2819 is behaviour-preserving — orient="support0" reproduces fault_contact's unconditional flip exactly; fault suites 603 passed, identical to before.

Three non-blocking items now fixed in a follow-up commit: the test.sh figures were mislabelled (1851 serial + 38 + 34 parallel = 1923 grand total; the body said "1923 serial", double-counting the 72 parallel tests — 1851 is the pre-response 1844 plus the 7 newly-wired test_1069 serial tests, which closes exactly); H1 disclosed only one of the two tolerances that moved (the velocity gate also went 1e-9→1e-8, equally justified since the [ti] spread is 1.05e-9 relative and would fail the old gate); and in test_rotated_freeslip_box_partition_independent the values[1] < 1e-3 assertion sat after the skippable anchor, so a different-host skip would have dropped it too.

Underworld development team with AI support from Claude Code

[lmoresi]

Follow-up commit eee9781 takes the three non-blocking items: the SolCx accuracy assertion now sits BEFORE the fingerprint-gated anchor (it would otherwise be skipped on a host whose mesh differs — the anchor skips by design, the accuracy check must not), and the body now states that TWO tolerances moved in test_1063, not one, with the velocity move justified by this PR's own measurement. test_1064 at np=2: 9 passed.

Underworld development team with AI support from Claude Code