Rotated free-slip: weight the nodal normal by the facet measure the assembly integrates over (#560)

[lmoresi]

Weight the rotated free-slip nodal normal by the facet measure

Fixes #560.

The mechanism, in two sentences

A boundary node that sits on more than one facet gets one nodal normal, while the
assembler integrates the boundary term facet by facet — and the normal we built was
the plain bisector of the adjoining facet normals, which is only right when those facets
are the same size. On a kinked wall with unequal facets the node's free tangential row
therefore keeps a residual sin(Δ/2)·(|f₁|−|f₂|)/6, the exact constant-pressure vector
stops being a null vector of the constrained rotated operator, PETSc projects the wrong
direction out of the right-hand side, and the pressure gauge is left unpinned — a silent
mean(p) ≈ 1e4 offset and a velocity whose last 13% is set by round-off, with
converged = True.

The weight, derived from the assembly

For a constant pressure p the boundary term the assembler produces in the row of
velocity basis function φ_i is p ∮ φ_i n ds, summed facet by facet. In the rotated
frame the free tangential row therefore carries

r_i = p · Σ_f ( ∫_f φ_i ds ) · ( t̂_i · n_f )

which vanishes for every tangent only when the nodal normal is parallel to
Σ_f ( ∫_f φ_i ds ) n_f. So the weight is the boundary basis integral, not the facet
count. For simplicial P1/P2 velocity that integral is the facet measure times a
constant that does not depend on which facet it is, so the measure is the correct
weight and the same argument settles 2-D and 3-D:

element	node type	∫_f φ_i ds	weight
2-D P2, edge facet	vertex (sees the kink)	|f|/6	length
2-D P2, edge facet	edge-midpoint	2|f|/3, one facet only	irrelevant
3-D P1, triangle facet	vertex	|f|/3	area
3-D P2, triangle facet	edge-midpoint (sees the kink)	|f|/3	area
3-D P2, triangle facet	vertex	0 — the row is consistent whatever normal it gets	any

Two things worth stating because they are not obvious. In 3-D the node that sees the kink
is the edge-midpoint, not the vertex: a boundary edge is shared by exactly two faces,
while the P2 vertex basis function integrates to exactly zero over a triangle, so the
vertex row is consistent regardless. And the constant in front of |f| cancels out of the
parallelism condition, which is why the measure works while the basis integral itself
does not (it would give a zero normal at every 3-D P2 vertex).

computeCellGeometryFVM already returns that measure — the change is to stop discarding
it. On a node whose facets all carry the same normal the weights are one common
positive factor, so the code skips the weighting there and keeps the historical unweighted
sum bit-for-bit.

Before / after

10×10 StructuredQuadBox, penalty=0, tolerance=1e-9,
petsc_use_pressure_nullspace=True, rotated free-slip on all four walls; the annulus rows
are radiusOuter=1.0, radiusInner=0.55, cellSize=0.1, both arcs rotated. Measured with
the issue's own probes (probe1_taper.py, probe2_nullspace.py, probe4_sweep.py,
probe7_areaweight.py, probe8b_annulus.py), unchanged, before and after.

config	|A z|/|A| before	after	principal angle before	after	mean(p) before	after
A straight box	2.733e-18	2.733e-18	5.24e-10 deg	5.24e-10 deg	6.504e-14	6.504e-14
B deformed top, tilted corner	5.902e-10	2.653e-18	1.11e-05 deg	3.96e-10 deg	+1.860e+04	4.055e-06
T tapered corner	2.139e-09	2.798e-18	3.39e-05 deg	4.46e-10 deg	+1.247e+04	6.476e-06
T2 C∞ bump, flat at the corner	1.887e-08	2.752e-18	2.90e-04 deg	4.52e-10 deg	+1.425e+03	1.011e-05
E annulus, uniform, geometric normal	2.843e-14	7.037e-19	3.54e-07 deg	6.52e-07 deg	−4.982e-04	−4.982e-04
E annulus, skew 0.01, geometric normal	1.503e-07	7.005e-19	2.08e-02 deg	6.21e-08 deg	−4.996e-04	−4.979e-04
E annulus, skew 0.04, geometric normal	5.650e-07	6.498e-19	3.69e-02 deg	3.18e-07 deg	−1.208e-02	−4.931e-04
E annulus, skew 0.04, analytic normal	1.130e-06	1.130e-06	6.50e-02 deg	6.50e-02 deg	−4.956e-04	−4.956e-04

(The uniform-annulus principal angle moves from 3.5e-07 to 6.5e-07 deg. Both are the
noise floor of a dense SVD on a 2-D near-null subspace whose second vector is a rigid
rotation; the column to read on that row is |A z|, which drops by five orders.)

The skewed annulus is the row that matters most: it keeps both boundaries exactly
circular and only moves the nodes to non-uniform angles, so nothing changes but the facet
lengths. Before the fix it acquires a misalignment thirty times worse than the deformed
box (6.5e-02 deg) and simply does not express it, because the skew also lifts the second
singular value out of machine zero and weakly pins the mode. After the fix, on the
geometric path, that misalignment collapses to 3e-07 deg — the fix addresses the
mechanism, not the symptom.

3-D, from the same argument, measured on a tet box (cellSize=0.25, rotated free-slip on
all six faces, top and side walls curved by z → z + amp·z·sin(πx), N=2474):

amplitude	|A z|/|A| before	after	mean(p) before	after
0.00 (flat control)	2.361e-18	2.361e-18	−3.403e-02	−3.403e-02
0.02	1.942e-06	2.439e-18	−1.912e+01	−3.439e-02
0.05	4.813e-06	2.559e-18	−7.840e+00	−3.481e-02

3-D is not "unverified by analogy": the face area is the weight the derivation asks
for, and it is measured. The 3-D defect is three orders of magnitude larger than the 2-D
one at the same amplitude, which is what the table above shows and what the P2
edge-midpoint argument predicts.

Amplitude sweep (shape B, res=10) — the convergence character, which is the point of the
|A z| metric:

amp	0.0005	0.002	0.005	0.01	0.02	0.05
before	9.37e-15	5.99e-13	9.34e-12	7.44e-11	5.90e-10	8.92e-09
after	2.72e-18	2.65e-18	2.72e-18	2.75e-18	2.65e-18	2.72e-18

Before, |A z| grew as amplitude³ (kink angle × facet-length difference, exactly the law).
After, it does not respond to the deformation at all. The run-to-run (2-eps) velocity move
over the same sweep falls from as much as 2.6e-01 to 4e-15 everywhere.

Do not judge this on mean(p)

mean(p) and the velocity move are the round-off amplitude of an unpinned direction:
in the pre-fix sweep, amplitude 0.05 moved by 7e-09 while amplitude 0.01 moved by 2.6e-01.
They carry the presence of the defect and not its size, and they can report a fix that is
luck. |A z|/|A| — the constrained operator applied to the attached null vector, normalised
by σ_max — scales as amplitude³ and as h³·³ and is what the regression tests assert.

What is bit-identical

• Straight walls. Every facet at such a node carries the same normal, the weights are
  a common factor, and the code takes the unweighted branch. Configuration A above is
  identical in every digit, including mean(p) and the 2-eps move.
• The analytic-normal override (normal= as a sympy matrix or a constant array). Its
  normal is a function of the node coordinate, not of the facet, so it also takes the
  unweighted branch. probe8b_annulus.py, which runs the analytic radial normal on the
  uniform and skewed annulus, is byte-identical before and after.
• normal="trace" and the fault-contact pair normals never reach this function
  (utilities/fault_contact.py has its own node walk) and are untouched.

The analytic path is therefore left with the consistency error the geometric path no longer
has — the last row of the table. That is deliberate: normal= is a user override that says
"use exactly this normal", and it is the right call on a genuinely curved boundary. It is
now stated in docs/developer/subsystems/rotated-freeslip.md, and the general guard the
issue's §8 proposes (attach the operator's actual null vector, or assert that whatever is
attached is one) is left as follow-up.

What changes for existing results

boundary_normal_traction / dynamic_topography project the nodal reaction onto this same
normal, so on a kinked geometric boundary the recovered σ_nn now uses the
measure-weighted normal. That is the assembly's own conjugate direction, so it is the more
correct choice; on flat walls and on analytic normals nothing moves. The free-surface
benchmarks are configuration C (the deforming boundary carries a traction, the rotated
free-slip lives on the straight side walls) and the spherical/annulus relaxation work is
configuration E with an analytic normal, so neither is touched — confirmed by running them.

Tests

New tests/test_1018_rotated_nodal_normal.py (level_1 / tier_a, one level_2 sweep):

• the deformed box keeps the constant pressure as a null vector, with the straight box as
  control;
• a negative control that re-runs the same measurement with the pre-fix bisector normal
  (frozen in the test file) and asserts it is ≥ 1e-11 — without it the assertions could pass
  on a metric that cannot see the defect;
• the skewed annulus on the geometric path, with its own negative control;
• bit-identity of the flat walls and of the analytic normal, with a control asserting that
  the two rules do differ on the kinked top (otherwise the comparison is vacuous);
• the amplitude sweep: no growth with deformation.

Fail-before, run against the unmodified branch: 4 failed, 1 passed. The four that fail
are the four that assert the fix; the one that passes is the negative-control test, which
patches the bisector normal in explicitly and therefore passes on both sides — that is the
point of it. After the fix: 5 passed.

Verified

Everything below on the fixed build, in the worktree's own environment, run sequentially.

what	result
tests/test_1018_rotated_nodal_normal.py (new)	5 passed
tests/test_1018_rotated_freeslip.py	22 passed
free surface: test_1070, test_1071, test_1072	6+1 skipped, 1, 2 passed
spherical / annulus: test_1011_stokesSph, test_1110_advDiffAnnulus	3, 1 passed
faults: test_0845 / 0846 / 0847 / 0848	6, 8, 6, 10 passed
parallel rotated (mpirun -n 2 --with-mpi): test_1064_rotated_freeslip_parallel, test_1066_rotated_datum_parallel	11 passed
scripts/test.sh — every CI serial batch	1845 passed, 17 skipped, 20 xfailed, 2 xpassed, 0 failed

The CI batches specifically: test_00[0-4]* 129; test_005[1-9]*+test_006[0-1]* 43;
test_01* 168; test_02* 27; test_05*+test_07* 532; test_08* (which is where the
fault suite lives) 605; test_100[0-9]* 30; test_101*+test_102* (the rotated free-slip
suite and the new file) 142; test_105* 83; test_1100* 5; test_1110*+test_1120* 4;
test_1450* 3; named tests 74. scripts/test.sh finished on "All tests passed".

One caveat, stated because it is a gap rather than a result. The whole-tree marker
invocation pytest -m "level_1 and tier_a" never reaches its run phase on this machine:
it sits in collection and then exits 1 having printed only the session header. That is
not caused by this change — the identical command, run in a separate worktree on
unmodified development (ffbcf98b), produces byte-identical behaviour. It is a
pre-existing property of collecting the whole tests/ tree in one process and deserves
its own issue; every test that gate would have selected is covered by the batches above,
which is how CI runs them.

Underworld development team with AI support from Claude Code

[Copilot]

Pull request overview

This PR fixes a rotated free-slip correctness issue on kinked/curved boundaries by changing how geometric nodal normals are accumulated so they match the facet-wise boundary assembly, restoring the constant-pressure nullspace compatibility and preventing an unpinned pressure gauge.

Changes:

• Update rotated free-slip geometric nodal-normal accumulation to be facet-measure weighted (edge length in 2D, face area in 3D), while preserving bit-identical behavior on coplanar facets.
• Add regression tests that capture the constrained rotated operator and assert the attached constant-pressure mode remains a null vector (with explicit negative controls using the pre-fix bisector rule).
• Document the measure-weighted normal rule and its implications for analytic normal= overrides.

Reviewed changes

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

File	Description
tests/test_1018_rotated_nodal_normal.py	New regression suite validating constant-pressure nullspace compatibility; includes negative controls.
src/underworld3/utilities/rotated_bc.py	Implements measure-weighted accumulation for geometric facet normals in rotated BC construction.
docs/developer/subsystems/rotated-freeslip.md	Documents the measure-weighted nodal-normal rule and analytic-normal behavior.

---

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

[lmoresi]

Adversarial review — PR #561, rotated free-slip nodal normal (fixes #560)

Reviewed at 2a8515f0 (3884eec1 / d735e18b / 2a8515f0), in an independent worktree
(.claude/worktrees/r561-review, own pixi amr-dev env, own build; uw.__file__ verified
in site-packages). Every number below is our own measurement in that environment, serial
unless a rank count is given. Our pre-fix control is an independent re-implementation of
the a25dd1e4 accumulation — we did not import the PR's test helper.

Verdict

The code change is right and we found no regression anywhere we looked. The derivation
is correct, it reproduces to four digits, and it strictly improves every configuration we
measured — including two the PR says it does not touch. What we are blocking on is the
claims shipped alongside it: the PR asserts a machine-exact geometric normal without the
word "serial", it adds documentation advice that our measurements show is wrong in
parallel, and its blast-radius statement about the free surface is factually inverted —
the free surface is in the broken configuration, its answers move by 24%, and the PR says
it is untouched. None of that requires re-deriving anything; it requires four edits and one
re-validation.

Request changes. Merge once MB-1..MB-3 are addressed.

---

MERGE-BLOCKERS

MB-1 — The fix does not survive partitioning, and nothing in the PR says so

_boundary_velocity_nodes accumulates over _boundary_stratum_is(...), which is the
rank-local facet list. There is no SF exchange and no reduction: a seam vertex is
summed from whatever its own rank happens to see. The derivation demands Σ over all
facets at the node; the code sums a subset.

Skewed annulus (skew=0.04, cellSize=0.15), geometric normal, |A z|/|A|_F from a PETSc
matvec — same formula at every rank count:

np	before	after
1	1.541e-07	6.619e-20
2	1.611e-06	1.604e-06
4	1.392e-05	1.336e-05

13 orders in serial; 4% at np=2 and np=4. On the deformed box at np=4 the fix is a
bit-for-bit no-op — |A z|/|A|_F = 8.204871e-08 identical to seven digits before and
after — and that value is 2000x worse than the serial pre-fix 3.9e-11.

Mechanism, pinned rather than inferred. At a seam vertex both adjacent exterior edges
are present in the rank's local mesh, but only one carries the boundary label:

annulus np=2  rank=1 wall=Upper coord=(0.71142568,-0.70276134)
              labelled-in-stratum=1  local-exterior-edges=2  of-which-labelled=1

So this is a label-distribution gap, not an overlap gap. Consequences we measured on the
deformed box at np=4: one owned node's normal differs from serial by 9.79e-03, and one
node's normal disagrees between two ranks by 1.958e-02. Affected seam vertices scale
with rank count — annulus: 4 at np=2, 8 at np=4, 16 at np=8.

Two things follow that the PR must not ship without.

• The geometric path is now rank-count dependent, while the analytic path is
  partition-independent and byte-identical (below). At np=4 on the skewed annulus the
  geometric normal is 43x worse than the analytic one (1.336e-05 vs 3.082e-07).
• The new test file is serial-only. test_skewed_annulus_constant_pressure_stays_a_null_vector
  asserts < 1e-15; the same measurement at np=2 reads 1.6e-06. We ran
  test_1064_rotated_freeslip_parallel + test_1066_rotated_datum_parallel at np=2 and
  np=4 — 11 passed each — and they pass because every one of them is either a straight box
  or an analytic normal. The PR's parallel evidence does not touch the changed code path on
  a curved boundary.

This is not a regression (before ≈ after at np≥2), so we are not blocking the code. We are
blocking the unqualified claim. Required: say "serial" where the PR says machine-exact; add
a parallel case recording the measured np=2/4 values (xfail is fine); open the seam-label
issue.

MB-2 — The PR adds documentation advice that is wrong in parallel

docs/developer/subsystems/rotated-freeslip.md (new text, lines 44-47):

Keep the facets near-uniform on an analytic-normal boundary, or use the geometric
normal there.

Measured, skewed annulus, |A z|/|A|_F:

np	ANALYTIC normal	geometric normal
1	3.082204e-07	6.619e-20
2	3.082204e-07	1.604e-06
4	3.082204e-07	1.336e-05

The analytic column is identical to seven digits at every rank count and byte-identical
before and after the fix. The advice is right in serial and inverted in parallel. It needs
the rank-count caveat.

Worse, the same file now contradicts itself 25 lines apart: line 20-21 still says an
analytic normal is "the preferred choice" on curved boundaries, and the same is true inside
rotated_bc.py — the new note at lines 119-146 argues one way while the surviving bullet at
lines 112-114 ("Best for exact curved/planar faces") argues the other. Fix both in place.

Six further shipped locations still carry the old advice and are untouched by the PR. The
two that matter most are user-facing:

• src/underworld3/cython/petsc_generic_snes_solvers.pyx:6141-6145 — the
  add_rotated_freeslip_bc docstring, "preferred on curved boundaries".
• .claude/skills/adapt-on-top-faults/SKILL.md:193-195 and CLAUDE.md:333-334 — the
  guidance every annulus/fault session is told to follow.

Given MB-1 the old advice is not actually wrong, so these are cheap to reconcile — but they
cannot be left contradicting the PR's own new paragraph.

MB-3 — The free surface IS configuration B, and its answers move by 24%

The PR states:

The free-surface benchmarks are configuration C (the deforming boundary carries a
traction, the rotated free-slip lives on the straight side walls) ... so neither is
touched — confirmed by running them.

That is wrong for the shipped FreeSurface class. With the Cartesian default normal=None
(free_surface.py:118,138), _build_held puts rotated free-slip on the deforming
surface with the geometric normal plus petsc_use_pressure_nullspace = True
(:562), _apply_walls replays the side walls as rotated free-slip (:316), and
_build_consistent does the same (:614). Its own docstring says so: "rotated u.n = 0
on every wall and the surface". That is configuration B. tests/test_1070 builds exactly
this at lines 62 and 132 — the Crameri Case-2 driver.

Measured directly (res=20 simplex box, essential Bottom, rotated Left/Right,
FreeSurface(normal=None), four steps at dt=0.02):

step	surf_amp before	surf_amp after	mean(held.p) before	after	|held.u| before	after
0	3.8133133449527734e-03	identical	2.347890e-02	identical	1.8180280009366026	identical
1	6.6608291626436000e-03	7.5751516748401304e-03	-1.546938e+07	2.166149e-02	1.4670825104735921	1.8160573877910593
2	1.0898250414321353e-02	1.1289493028559860e-02	-1.424484e+04	1.966209e-02	1.7872988565970607	1.8139258021073050
3	1.5136564353666282e-02	1.4954782899239905e-02	-1.587650e+03	1.760255e-02	1.7873763741703950	1.8115955965299231

Step 0 (surface still flat) is bit-identical — the PR's flat-wall claim holds in a real
workflow. From step 1 the pre-fix held lid loses the gauge (mean(p) runs 1e7 → 1e4 → 1e3,
chaotic, exactly the #560 signature) and the held velocity is 23.8% different. The
nodal normal on Top moves by at most 1.45e-05 (mean 2.09e-06) — a 1e-5 normal change
buying a 24% velocity change is the unpinned-mode amplification, and it is the reason
"we ran the tests and they passed" cannot settle this.

So this PR fixes the free surface, which is a better result than it claims. But
h_∞ is the recovered reaction of that held solve, so every free-surface and
dynamic-topography number changes. The PR must say so, and the Crameri numbers need
re-validating rather than re-running for pass/fail — test_1070's assertions are soft bands
(:94 factor-5, :97 2%, :165 a decade) and will not notice.

---

HIGH

H-1 — The scope claim is right for simplices and for all of 2-D; over-stated for 3-D hexes

The docstring says "the measure IS the right weight and the same argument covers 2D and
3D". It covers simplicial facets. A 3-D hex boundary facet is a bilinear quad, and once it
is not a parallelogram ∫_f φ_i ds is no longer measure × a facet-independent constant.
A deformed hex box is an ordinary UW3 case. Measured (z → z + amp·z·sin(πx), rotated
free-slip on all six faces):

case	before	after	flat control
3-D TET, amp 0.02	5.972e-07	8.329e-19	8.078e-19
3-D TET, amp 0.05	1.480e-06	8.285e-19	8.078e-19
3-D HEX, amp 0.02	1.066e-08	3.445e-17	5.080e-18
3-D HEX, amp 0.05	1.610e-07	8.345e-17	5.080e-18

Tets return exactly to their flat control. Hexes do not: the remainder is 7x and 16x the
flat control and grows with amplitude, which is the non-affine facet showing through. It is
not a regression — still 3.1e8 and 1.9e9 better than the bisector, and the crossover below
which the bisector would win is amplitude ≈ 1e-6. Reword the docstring to "simplicial
facets, and any 2-D facet; on non-affine 3-D quad facets the measure is the leading-order
weight". No code change needed.

For completeness on the other two scope axes we were asked about:

• Degree. In 2-D a boundary facet is a 1-D Lagrange element at any degree, so
  ∫_f φ_i ds = |f| × (node-type constant) and the parallelism condition is degree-free —
  the derivation extends. We could not confirm it by measurement, because degree-3
  velocity gives |A z|/|A| = 3.2222e-04 on a straight box, identically at qdegree=3
  and qdegree=6, and identically before and after the fix. That is a separate, pre-existing
  failure of the rotated constant-pressure nullspace at P3 (14 orders off the P2 straight-box
  value of 2.7e-18); it is not this PR's doing and not fixed by it. Worth its own issue.
• Curved / isoparametric facets. Not reachable here: computeCellGeometryFVM returns a
  straight-facet measure and normal regardless, so the weight is the chord measure. This is
  the same approximation the assembler makes for the P2 boundary integral, so it is
  self-consistent — but it is the reason the analytic normal keeps a consistency error
  (measured 3.082e-07 on the skewed annulus, unchanged), which the PR correctly identifies
  and defers.

H-2 — fault_contact carries a second, unfixed copy of the same bisector

src/underworld3/utilities/fault_contact.py:551-566 accumulates unit facet normals and
normalises — verbatim the pre-#560 code this PR replaced, in a second copy the PR neither
fixes nor mentions. Its own docstring (:83-93) describes the symptom ("on a SAMPLED CURVE
zig-zags at node scale") and docs/advanced/split-node-faults.md:238-266 tells users to work
around it with an analytic normal. tests/test_0847_fault_api.py:225 puts add_fault_bc
on a curved "Arc" with no normal. Out of scope for this PR, but it should be named in the
PR body and get a follow-up issue rather than being left for someone to rediscover.

---

MEDIUM

M-1 — The existing deform test flips branch; the promised cleanup is not done

tests/test_1018_rotated_freeslip.py::test_rotated_workspace_deform_invalidates gates on
if unpinned_size > 1.0e-3:. We measured it directly: unpinned_size = 9.926848e-15, so
it now takes the else branch and asserts err < 1e-6 on the raw difference — the test is
back at full strength, which is a real win and should be in the PR body. But its own else
message says "Delete the projection and the _deformed_box_rotated_solve probe" and that
was not done, so it still runs an extra Stokes solve for nothing, and the file still carries
pre-fix prose as current fact (:534-546, "Measured: ... moves the velocity by 1.33e-01").

M-2 — The undeformed annulus on the geometric path changes; benign, and nothing depends on it

We confirmed the PR's own table row rather than taking it on trust: uniform annulus,
geometric normal, 9.5737e-14 → 1.1343e-18. Not bit-identical, and correctly so — an
"unstructured uniform" annulus does not have exactly equal facets. We audited the suite for
recorded values that could ride on it: every annulus / spherical / cubed-sphere golden
(test_1064 GOLDEN_ANNULUS, GOLDEN_ANNULUS_FMG, GOLDEN_SPHERICAL3D,
GOLDEN_SPHERICAL3D_TOPO; test_1066 _INT_VV_REF; test_1018_rotated_freeslip:761-766
0.41920 / 0.77060) is on the analytic path, and every geometric-path golden
(GOLDEN_BOX, GOLDEN_BOX_SIGMA, GOLDEN_TOPO_BDL2, GOLDEN_BOX_NONLINEAR) is on a
straight box. So the change is benign and we found nothing that needs re-baselining.
Two shipped scripts do move: scripts/aniso_convection_demo.py:87 and
scripts/adapt_vs_ref24_cost.py:57 take the geometric normal on an Annulus.

M-3 — Test-metric discipline: controls fire, but one assertion can pass vacuously

We checked this rather than assuming it.

• The negative controls are real. Our independent bisector control reproduces the pre-fix
  value to four digits (5.9015e-10 against the PR's 5.902e-10), and the existing deform
  test's own injection control (poisoned_err > 1e-4) still fires.
• The metric cannot silently read machine-zero. If the capture never fires,
  self.nsp is None and nsp[0] raises TypeError; if _rotated_nullspace returns
  None, kw["nsp"].getVecs() raises AttributeError; a zero null vector gives nan,
  and nan < 1e-15 is False. All three fail loudly.
• One genuine hole. test_flat_walls_and_analytic_normals_are_bit_identical:238-243
  (the annulus/analytic block) has no assert fixed and fixed.keys() == old.keys() guard,
  unlike the straight-box block above it at :228. An empty stratum — a relabelled mesh, a
  typo'd arc, a rank owning no arc — makes both iterations pass with zero comparisons,
  and this is the only test in the PR that exercises the analytic path at all. Add the guard.
• _MACHINE = 1e-15 is an absolute constant gating three configurations with a control
  value recorded for only one (the straight box, 3e-18). We measured the other two:
  annulus uniform 1.134e-18, skew 0.04 1.049e-18 — 1000x headroom, fine today. A threshold
  relative to the measured flat control would not need re-checking on the next mesh.

---

LOW / NITS

• L-1. "Flat walls are unchanged to the last bit" is false for a tilted flat wall. The
  shortcut is np.array_equal on computeCellGeometryFVM normals — exact float equality.
  Measured node counts, before vs after: axis-aligned structured quad box 68/68 identical;
  2-D unstructured simplex box 68/68; 3-D unstructured simplex box 390/390; 2-D box
  rigidly rotated 30°: 20 of 68 nodes move, max |dn| = 1.570e-16 (one ulp). Harmless in
  magnitude, but the claim should say "coplanar facets whose computed normals agree bitwise"
  — which is not every flat wall.
• L-2. Latent RecursionError in the frozen control.
  test_1018_rotated_nodal_normal.py:38-39 delegates the analytic path to
  rbc._boundary_velocity_nodes, which under _BisectorNormals is _bisector_nodes. We
  hit this for real in our own copy of the helper before binding the original at import. No
  current test triggers it, purely by ordering accident — and the obvious next test (an
  analytic-normal solve under the control, given the PR's own thesis) detonates it. One-line
  fix: def _bisector_nodes(..., _real=rbc._boundary_velocity_nodes).
• L-3. nsp[0] is the constant pressure only because it is appended first
  (rotated_bc.py:1757 before the rotation loop at :1766). On the annulus nsp has two
  vectors. Load-bearing and uncommented in the test.
• L-4. Cost for the tier. The new file is level_1 / tier_a but runs ~14 Stokes solves
  and dense O(n³) SVDs; test_1018_rotated_nodal_normal.py + test_1018_rotated_freeslip.py
  took 143 s together. Three of those solves (:223, :234) exist only to populate a DM and
  their results are discarded.
• L-5. The pre-PR docstring already claimed the accumulation was "area-weighted"
  (preserved in docs/docstrings/inventory.json:65694) when it was not. Worth one line in
  the PR body — the wrong docstring is part of why this survived.
• L-6. No 3-D and no simplex coverage in the new test file. The docstring makes three
  distinct 3-D claims (P1 vertex |f|/3, P2 edge-midpoint |f|/3, P2 vertex integral
  identically zero); none is tested, and the test's own box is a StructuredQuadBox while
  the justification is stated "for simplicial P1 and P2". Our H-1 table is the missing
  coverage — a 3-D tet case is cheap and would have caught the hex caveat.

---

What we verified of the PR's own claims

Independently reproduced, in our own worktree and environment:

claim	PR	ours	verdict
2-D deformed box B, before	5.902e-10	5.9015e-10	✓
2-D deformed box B, after	2.653e-18	2.6241e-18	✓
straight box 2-D bit-identical	2.733e-18 both	2.7437e-18 both, all digits	✓
straight box 3-D bit-identical	2.361e-18 both	8.0783e-19 both; normals 390/390 identical	✓
analytic override byte-identical	byte-identical	1.9147e-13 and 4.7248e-06, identical before/after and at np=1,2,4	✓
uniform annulus geometric changes	2.843e-14 → 7.037e-19	9.5737e-14 → 1.1343e-18	✓ benign (M-2)
3-D tet, amp 0.02 / 0.05	1.942e-06 / 4.813e-06 → 2.4e-18	5.972e-07 / 1.480e-06 → 8.3e-19	✓ same character
test_1018_rotated_nodal_normal + test_1018_rotated_freeslip	5 + 22	27 passed (143 s)	✓
free surface + spherical (1070/1071/1072/1011/1110)	6+1s, 1, 2, 3, 1	13 passed, 1 skipped (222 s)	✓
CI batch test_101* + test_102*	142	142 passed (319 s)	✓
parallel rotated (1064 + 1066)	11 at np=2	11 passed at np=2 and 11 at np=4	✓ (but see MB-1)
CI batch test_08* (the fault suite)	605	605 passed, 11 skipped, 11 xfailed (763 s)	✓

Four CI batches run end-to-end here (test_101*+test_102*, test_08*, the free-surface
and spherical files, both rotated files) plus the parallel rotated suite at two rank counts.
Every count the PR quotes for those matched ours exactly. We did not re-run the other nine
batches; nothing we found suggests they would move.

Marker-collection gap: reproduced. pytest -m "level_1 and tier_a" over the whole
tests/ tree exits 1 having produced zero bytes of output — and it does so even with
--co -q, i.e. it dies in collection, not in the run phase. That matches the implementer's
report; we did not chase it, and we agree it deserves its own issue. Separately, another
session was running pytest -m level_1 -n 8 on this machine throughout, which is worth
knowing before anyone reads timings off this review.

---

Artifacts

• Review: /private/tmp/claude-501/-Users-lmoresi--Underworld-underworld3-pixi/64e34dd9-4af2-4774-9d55-4238a998eb6e/scratchpad/r561_review.md
• Probes and logs: /private/tmp/claude-501/-Users-lmoresi--Underworld-underworld3-pixi/64e34dd9-4af2-4774-9d55-4238a998eb6e/scratchpad/r561/
  — common.py (independent pre-fix control + |A z| capture, serial and PETSc-matvec),
  probeA_repro_2d.py, probeB_scope_3d.py (tet vs hex),
  probeC_parallel.py + analyse_par.py (cross-rank normals),
  probeD_seam.py (label vs mesh), probeE_annulus_par.py (np sweep, both normal sources),
  probeF_freesurface.py (held lid), probeG_bitident.py (tilted / simplex flat walls),
  branch_check.py.
• Worktree: /Users/lmoresi/+Underworld/underworld3-pixi/.claude/worktrees/r561-review
  (branch bugfix/r561-review, reset to pr561-head).

Nothing was posted, pushed or merged.

Underworld development team with AI support from Claude Code

[lmoresi]

Response commits fc81ba32, 91dc8fda, 7eb92dd5, f6e27044.

MB-1 fixed — and it was hiding two more defects. The weighted sum is now completed across ranks through the DM's own local↔global scatter (ADD, then scatter back), collective on every rank including those owning no facet. Fixing it exposed: (a) the outward test was rank-local — it oriented facets away from the mean of that rank's coordinates, so two facets at a seam node could be oriented oppositely and cancel in the sum; now oriented away from the facet's own support cell, which also fixes non-convex domains where the coordinate mean sits in the hole; (b) a rank can own a node whose labelled facets all live on neighbours and got no constraint row at all — the node list now comes from the local mesh's exterior facets. The co-planar shortcut is gone: it decided co-planarity from the facets one rank could see, which is exactly what a seam splits.

|A z|/|A|_F before → after:

boundary	np=1	np=2	np=4
uniform annulus	5.71e-15 → 6.62e-20	1.49e-06 → 6.78e-20	1.40e-05 → 6.72e-20
skewed annulus 0.04	1.54e-07 → 6.62e-20	1.61e-06 → 6.62e-20	1.39e-05 → 6.69e-20
deformed box	3.92e-11 → 1.73e-19	3.92e-11 → 1.73e-19	8.20e-08 → 1.73e-19

Read the uniform annulus first: the configuration everyone treats as clean was the worst one in parallel. The seam defect does not need a kink, it needs a partition. Shared-node normals now agree bitwise across ranks (was max |Δn| 1.30e-01 at np=2, 1.63e-01 at np=4); tests/parallel/test_1067_* asserts both, with those values as fail-before.

MB-3 — the free surface is configuration B, and Crameri is safe. Reproduced the review's numbers exactly (mean(held.p) −1.55e+07 → 2.17e-02, |held.u| +23.8%, and |held.u| now smooth where it used to jump). But the benchmark itself is unchanged to ten significant figures: _case2(res=40), 8 steps, amplitude in metres — 360.01, 376.51, 392.31, 404.88, 419.65, 433.80, 450.58, 463.86, identical before and after, crossing the published ≈398 m checkpoint between steps 2 and 3 as before. Case-2's slope is ~1e-4, far below where the weighting bites. Agreement unchanged, not degraded; no re-baselining required. PR body corrected — free surface is touched, and improved.

H-2 — fault_contact fixed, and it was more exposed than a wall. A constant pressure gives F⁻ = −F⁺ exactly, so it cancels in the mean rows and doubles in the jump rows — and the jump-tangential row is the slip row. The bisector left √2·p·sin(δ/2)·(|f₁|−|f₂|)/6 there: lost gauge plus pressure-driven spurious slip scaling with lithostatic pressure. The existing leak diagnostic contracts with the same normal that defined the constraint, so it is true by construction and cannot see it. Fault pair nodes are rank-local by construction, so no cross-rank sum is needed there.

MB-2 — re-measured after MB-1 and written as the measurement, reconciled in place across all six locations rather than adding a seventh statement: geometric 6.6e-20 at every rank count; analytic 1.1e-14 (uniform) / 3.1e-07 (skewed), unchanged by the fix. The advice now says the two normals answer different questions — assembly-consistent versus geometry-consistent. H-1 — docstring now claims exact for simplices and all of 2-D, leading-order on non-affine 3-D quads. M-3/L-1/L-2/M-1 all addressed ("flat" corrected to "axis-aligned"; the deform test's now-redundant projection removed).

Verification: 27 rotated · 13 parallel at np=2 and np=4 · faults 6/8/6/10 · free surface 6+1skip/1/2 · scripts/test.sh 1845 passed, 0 failed.

One gap found and not fixed here: the nondimensional→Myr time scale is recorded nowhere (hscale is, tscale is not), so a published checkpoint like "398 m at 3 Ma" can only be checked as a trajectory, not at matched physical time. Filed separately.

Underworld development team with AI support from Claude Code

[lmoresi]

Adversarial re-review — PR #561, response commits (fixes #560)

Re-reviewed at f6e27044 (fc81ba32 / 91dc8fda / 7eb92dd5 / f6e27044 on top of the
original three), in the same independent worktree as round 1
(.claude/worktrees/r561-review, own pixi amr-dev env, rebuilt at the new head;
uw.__file__ verified in site-packages). Every number below is our own measurement.
Our pre-fix control is still our own re-implementation of the a25dd1e4 accumulation —
and this round we also ran the merge-base modules themselves (a25dd1e4's
rotated_bc.py + fault_contact.py copied over site-packages) as a negative control.

Verdict

The parallel mechanism is sound, and we verified it rather than accepting it. We
attacked the cross-rank sum on four geometries at np = 1, 2, 3, 4, 8 and could not break
it: no facet is double-counted, none is labelled away from its owner, ranks owning zero
boundary facets take part without hanging, every shared node's normal is bit-identical
across ranks, and no node loses its constraint row. The fault claim, the Crameri claim
and the deleted-shortcut claim all reproduce exactly. The three round-1 merge-blockers
are properly closed.

What we are blocking on is again what the PR does not say, plus one wiring gap. The
new outward test silently flips the sign of σ_nn on every concave boundary — the
annulus/spherical-shell inner arc, i.e. the CMB — and leaves the geometric and analytic
normals disagreeing in sign there while the docs still recommend X/|X|. And the one
parallel test that guards this entire fix lives in a directory the project's own test
script does not run.

Request changes. Merge once MB-1 and MB-2 are addressed — both are edits, not
re-derivations.

---

MERGE-BLOCKERS

MB-1 — The new outward test flips σ_nn on every concave boundary, silently

fc81ba32 replaces "away from the mean of the mesh coordinates" with "away from the
facet's own support cell". On a convex boundary those agree. On a concave one
they are opposite, and the commit message sells this only as a parallel-correctness
improvement ("also correct on a non-convex domain where the coordinate mean sits in the
hole"). It is more than that: it is a user-visible sign change in a shipped reading.

Annulus (radiusOuter=1.0, radiusInner=0.55, cellSize=0.2), rotated free-slip on
both arcs, mean radial component of the nodal normal:

arc	NEW (f6e27044)	OLD (a25dd1e4)	analytic X/|X|
Upper (convex)	+1.000000	+1.000000	+1.000000
Lower (concave)	−1.000000	+1.000000	+1.000000

boundary_normal_traction is R_i = n̂_i · r_c and its docstring says "n̂_i is THIS
boundary's outward normal" (rotated_bc.py:1970; :103 likewise). So the sign flips
with it — measured, same mesh, same solve, first-5-node mean of σ_nn:

boundary	before	after	analytic override, after
Upper	+5.810571e-02	+5.810571e-02	+5.810571e-02
Lower	−5.233110e-02	+5.233110e-02	−5.233110e-02

Three things follow, none of which is in the PR body, the docstrings or a test.

• Dynamic topography on a concave boundary flips sign against the merge base.
  dynamic_topography_field is h = −σ_nn/(Δρ g) on top of exactly this number. Any
  annulus/shell run reading the inner boundary on the geometric path changes sign.
• The two normal sources now disagree in sign on a concave boundary — magnitude
  identical to every digit, sign opposite. The same file (rotated-freeslip.md), the
  add_rotated_freeslip_bc docstring and the adapt-on-top-faults skill all still offer
  X/|X| as the analytic normal "on a spherical cap"; on a shell's inner boundary that
  is inward, and now contradicts the default.
• A non-zero prescribed wall-normal datum reverses there. build_rotation:423-426
  anchors the datum sign to M.sum(axis=0) — the accumulated normal — so
  add_rotated_freeslip_bc(+u_n, "Lower") means outflow after and inflow before.

We think the new behaviour is the right one: it is what the docstrings have always
claimed, and it is what makes the geometric path self-consistent. That is why this is a
"say it and reconcile it" blocker and not a "revert it" one. Required: state the sign
change in the PR body; say in rotated-freeslip.md that the geometric normal is the
domain's outward normal on concave boundaries while an analytic X/|X| is not, so the
two give opposite-signed σ_nn there; and add one assertion (an annulus inner arc, radial
component < 0) so the convention is pinned by a test rather than by a comment.

MB-2 — The parallel test that guards this fix is not run by the project's test script

tests/parallel/test_1067_rotated_nodal_normal_parallel.py is the response's evidence
that MB-1 (round 1) is closed. scripts/test.sh --p N runs exactly one parallel glob:

mpirun -n $PARALLEL_RANKS python -m pytest --with-mpi tests/parallel/test_075*py
# Add other parallel test categories as they're created:
# mpirun -n $PARALLEL_RANKS python -m pytest --with-mpi tests/parallel/test_10*py

The second line is commented out, and .github/workflows/build_uw3_and_test.yaml:68
runs ./scripts/test.sh --p 2. So the new test — and test_1062…test_1066 with it —
never executes in CI at any rank count. The un-wiring predates this PR; shipping the
fix's only regression guard into it does not.

We confirmed the test is a real fail-before guard: with the merge-base modules installed
over site-packages it fails at np=2 with "4 of 4 shared boundary nodes disagree between
ranks, max |dn| = 2.215e-01", and passes at np=2, 3, 4 at the PR head. It is worth
running. Uncomment the line (or add test_1067* specifically) so it does.

---

HIGH

H-1 — _local_boundary_candidates states an invariant that is measurably false

Its docstring justifies the whole helper with:

an interior partition facet keeps support 2 through the overlap, so support 1 still
means the domain boundary.

It does not. The facets on the outer fringe of the overlap layer have support 1 locally
and are interior globally. Measured — candidates mapped to coordinates and checked
against the true (serial) boundary-node set:

mesh	np	candidates/rank	of which not on the domain boundary
skewed annulus	2	78, 78	10, 10 (12.8%)
skewed annulus	4	48, 48, 48, 44	12 each (25–27%)
3-D tet box	4	170, 178, 190, 190	34, 32, 31, 43 (16–23%)

The result is still correct, but not for the stated reason: those points receive no
contribution from any labelled facet, so if w.any() drops them. That filter is the
load-bearing safety property and it is documented as an aside about constrained DOFs.
Fix the docstring to say what actually holds ("support 1 over-collects at the overlap
fringe; the zero-sum test discards the extras"), because the next person to touch the
node set will trust the sentence that is there now.

H-2 — The constrained-DOF fallback is justified for one consumer and used by two

_sum_facet_normals_across_ranks:303-307 falls back to the rank-local contribution when
a node's velocity DOFs are constrained out of the global vector, on the grounds that
"build_rotation skips such nodes anyway". That is true — build_rotation:402 drops
l2g < 0. But boundary_normal_traction:2009 consumes the same list, reads
n̂·r_c at every returned node including those, and hands them to _desmear with
partial_reaction=False (overwrite-by-coordinate). If such a node is shared, two ranks
overwrite the same coordinate with normals summed from different facet subsets.

This configuration is shipped, not hypothetical — essential BC on Bottom + rotated
free-slip on the other walls is tests/test_1070::_case2. We measured it (simplex box
1.0 × 0.25, cellSize=0.25/6): 2 nodes with l2g < 0 at np = 2, 3, 4 and 8, at
every rank count. We could not make them disagree — at the partitions we got, the two
constrained corners were never shared, so all 75 node normals stayed bit-identical across
ranks and against serial. So this is a documented-reason gap rather than a demonstrated
defect: either extend the comment to cover boundary_normal_traction, or drop
constrained nodes from the returned list outright (nothing downstream wants them).

---

MEDIUM

M-1 — The subsystem doc still carries the claim the same commit series corrected

f6e27044 correctly narrowed "flat walls are unchanged to the last bit" to
axis-aligned in rotated_bc.py:153-156 and in
tests/test_1018_rotated_nodal_normal.py:227. docs/developer/subsystems/rotated-freeslip.md:40-41
still says:

Flat walls are unchanged to the last bit: every facet there shares a normal, so the
weighting cancels.

Both halves are now wrong: "flat" should be "axis-aligned", and the reason describes the
deleted co-planarity shortcut ("the weighting cancels"), not the surviving mechanism
(0/±1 components normalising identically). One-line edit, but it is the user-facing
document and it contradicts the docstring 100 lines away.

M-2 — The boundary setup got 1.5–5.7× more expensive, and no longer scales down

3-D tet box, cellSize=0.12 (7505 local facets, 1689 labelled boundary facets), all six
rotated boundaries, best of 3:

np	BEFORE (a25dd1e4)	AFTER (PR)	ratio	of which _local_boundary_candidates
1	9.5 ms	14.3 ms	1.51×	n/a (skipped at size 1)
2	4.8 ms	23.2 ms	4.87×	6.4 ms (28%)
4	2.5 ms	14.3 ms	5.71×	3.7 ms (26%)

The before column halves with each doubling of rank count; the after column does not,
because _local_boundary_candidates walks the entire local facet stratum (with a
getTransitiveClosure per support-1 facet, in Python) and the scatter is over the whole
DM vector — both O(local mesh) per boundary, where the old code was O(boundary facets).
Mitigated in practice by the rotation cache (rotated_bc:1061-1068, and its verdict is
made collective by an allgather), so this only bites when the geometry changes — but
that is exactly the free-surface / adaptive loop. Worth a note; not worth blocking. The
candidate walk can be restricted to the closure of the labelled facets' vertices' stars
if it ever matters.

M-3 — The fault rank-locality claim is verified in 2-D only

We tested the configuration the claim has to survive: the spanning fault
([0.06,0.08] → [0.94,0.92], the geometry ptest_0845 uses precisely because it crosses
any contiguous partition), built through the user path mesh.add_fault, with
add_frictionless_fault_bc + solve_with_fault.

np	pair nodes	shared across ranks	normals differing from serial	local facet-count mismatches
1	63	—	—	—
2	63	0	0 (max |dn| 0.0)	0
3	63	0	0	0
4	63	0	0	0

Not one pair node is shared at any rank count, and every node sees exactly its serial
complement of fault facets locally — so the measure weighting is complete rank-locally
and 91dc8fda needs no cross-rank sum, as claimed. build_rotation:471-482 also keeps
an explicit RuntimeError if a pair ever straddles. Claim verified. The residual is
that we did not repeat this in 3-D (ptest_0848); the 3-D splitter uses the same
redistribution, but the PR asserts the property for both and only 2-D is measured here.

---

LOW / NITS

• L-1. Cross-rank identity is exact; serial-vs-parallel is exact except one 3-D ulp.
  On the 3-D tet box at np=4, 4 of 618 (wall, node) pairs differ from the serial normal by
  3.469e-18 — the scatter's ADD order is not the serial dict's. Every other geometry
  and rank count is 0.0 to the last bit, and cross-rank agreement is 0.0 everywhere,
  which is what the new test asserts. Nothing to fix; worth knowing before someone writes
  a serial-vs-parallel bit-identity assertion.
• L-2. Test number collision. tests/parallel/test_1067_rotated_nodal_normal_parallel.py
  and tests/test_1067_newton_cold_start.py both claim 1067.
• L-3. Deleting the shortcut moved the tilted-wall count. Round 1 measured 20 of 68
  nodes moving by 1 ulp on a 30°-tilted flat wall (shortcut present); with it gone we
  measure 27 of 68, max |dn| unchanged at 1.570e-16 — which is the number
  fc81ba32's message quotes. Magnitude unchanged; the claim is stated correctly.
• L-4. A node whose weighted normals cancel is now silently dropped
  (return ... if np.linalg.norm(v) > 0.0). Pre-fix such a node reached the SVD and got
  r=0, i.e. also no constraint — so no behaviour change, but the drop is now invisible
  where it used to be traceable.
• L-5. _boundary_stratum_is(...).getSize() on a null IS segfaults — we hit it in our
  own probe at np=4. The PR code guards it correctly everywhere; a helper returning an
  empty IS instead of NULL would stop this recurring (this is the Stokes_Constrained segfaults at np>1 in the interior-multiplier section reduction #291 class).
• L-6. test_1067's two tests both build and solve the annulus (2 Stokes solves for
  what is really one measurement); _ConstantPressureResidual duplicates the serial
  file's capture with a different normaliser, correctly flagged in its own docstring.

---

What we verified of the response's own claims

Independently reproduced in our worktree; verified, not accepted.

claim	response	ours	verdict
no boundary facet labelled twice	asserted	0 duplicates, np=2/3/4/8, on 4 geometries	✓
none labelled away from its owner	asserted	0, all runs (checked against the point SF)	✓
labelled facet count = serial	—	66/66, 40/40, 258/258, 36/36 at every np	✓
reduction is collective (rank with no facet)	asserted	at np=8, 6 of 8 ranks own no Left facet, 7 of 8 no Right, one owns no Top — no hang, complete normals	✓
shared-node normals bitwise identical	asserted	max |dn| = 0.0 at np=2/3/4/8, geometric and analytic	✓
no node loses its constraint row	asserted	rows = serial exactly (132/132, 84/84, 618/618, 73/73), 0 missing, 0 duplicated	✓
skewed annulus |A z|/|A|_F	6.62e-20 (1), 6.62e-20 (2), 6.69e-20 (4)	6.624e-20, 6.623e-20, 6.690e-20; and np=3 6.629e-20, np=8 6.713e-20	✓
uniform annulus was the worst in parallel	5.71e-15 / 1.49e-06 / 1.40e-05 → ~6.7e-20	5.711e-15 / 1.485e-06 / 1.396e-05 → 6.62/6.78/6.72e-20	✓
non-convex orientation now correct	asserted	inner-arc radial component −1.000000 on all 48 nodes, 0 wrong-signed	✓ (see MB-1)
shortcut deletion changes no axis-aligned node	68/68, 68/68, 390/390	68/68, 68/68, 390/390, max |dn| 0.0	✓
fault pair nodes rank-local	asserted	63/63, 0 shared, 0 differing, np=2/3/4 on a spanning fault	✓
Crameri unchanged to 10 s.f.	360.01 … 463.86	360.01, 376.51, 392.31, 404.88, 419.65, 433.80, 450.58, 463.86; rel. diff before/after 7.9e-11 … 3.1e-10	✓
test_1067 at np=2/4	13 parallel	2 passed at np=2, 3 and 4; fails on merge-base modules (negative control fires)	✓
rotated parallel suite	—	12 passed at np=2 and at np=4 (1062+1064+1066)	✓
test_101*+test_102*	27 rotated	142 passed (197 s)	✓
test_08* (faults)	6/8/6/10	605 passed, 11 skipped, 11 xfailed (348 s)	✓
free surface + spherical (1070/1071/1072/1011/1110)	6+1skip/1/2	15 passed, 1 skipped (117 s)	✓

One parallel failure, and it is not this PR.
tests/parallel/test_1063_constrained_freeslip_parallel.py::test_constrained_freeslip_partition_independent[iso]
and [ti] fail at np=2 and np=4 — velocity L2 0.6194547793955 (golden) vs
0.5982807168536982 (np=2), a 3.4% partition dependence in the multiplier free-slip
path. We attributed it rather than reporting it: with a25dd1e4's rotated_bc.py and
fault_contact.py copied over site-packages the failure reproduces to the last digit,
and the whole PR touches only those two modules plus one .pyx docstring.
Stokes_Constrained never calls _boundary_velocity_nodes. Pre-existing; needs its
own issue. (The same swap makes test_1067 fail, which is how we know it took effect.)

---

Carried forward from round 1 — not re-verified

Cleared last round and untouched by the response commits, so we did not re-measure:

• the derivation itself, and the serial |A z| numbers for configs A / B / T / T2 — the
  scope wording is now corrected exactly as H-1 asked (exact for simplicial facets and
  all of 2-D, leading-order on non-affine 3-D quads), and the deform test now passes at
  full strength inside the 142-test batch;
• the 3-D hex caveat (H-1) — the docstring now states it; not re-measured;
• the analytic path's serial consistency error (3.082e-07 skewed / 1.9e-13 uniform);
• L-3 (nsp[0] ordering), L-4 (test cost), L-6 (no 3-D coverage in the serial file).

Re-checked because a response commit touched them:

• straight-wall bit identity — the co-planarity shortcut is deleted, so we re-ran it
  (68/68, 68/68, 390/390, 27/68 tilted at 1 ulp);
• analytic-path byte identity — re-verified across ranks and against serial, since the
  path now shares the node-set code with the reduced one;
• the uniform-annulus geometric change (round-1 M-2) — re-measured, now 5.711e-15 →
  6.621e-20 serial, and partition-independent;
• round-1 M-3's two test holes — both closed (assert fixed and fixed.keys() == ... on
  the analytic block; _real=rbc._boundary_velocity_nodes bound at import);
• round-1 M-1 — the projection, the probe solve and the pre-fix prose are gone; the
  surviving mentions at :404-407 and :515-525 are correctly written as history.

---

Artifacts

• Re-review: /private/tmp/claude-501/-Users-lmoresi--Underworld-underworld3-pixi/64e34dd9-4af2-4774-9d55-4238a998eb6e/scratchpad/r561_rereview.md
• Round 1: .../scratchpad/r561_review.md
• Probes and logs: /private/tmp/claude-501/-Users-lmoresi--Underworld-underworld3-pixi/64e34dd9-4af2-4774-9d55-4238a998eb6e/scratchpad/r561/
  — pa_audit.py + pa_report.py (label / orientation / node / row audit, 4 geometries
  × np 1–8), pb_sign.py and pc_concave.py (the concave sign flip and σ_nn),
  pd_cand.py (spurious candidates), pf_crameri.py (+ crameri_before.json,
  crameri_after.json), pg_fault.py (spanning fault, np 1–4), ph_cost.py,
  probeE_annulus_par.py (metric sweep), probeG_bitident.py, common.py (the
  independent a25dd1e4 control), logs g_101x.log, g_08.log, g_fs.log,
  g_par2.log, g_rot{2,4}.log.
• Worktree: /Users/lmoresi/+Underworld/underworld3-pixi/.claude/worktrees/r561-review
  (branch bugfix/r561-review, reset to pr561-head2 = f6e27044). The site-packages
  module swap used for the negative control was restored and byte-compared against the
  worktree source before the final runs.

Timing caveat: another session was running pytest -m level_1 -n 8 on this machine during
part of this review, so read M-2's absolute milliseconds as best-of-3 under contention;
the ratios were measured back to back.

Nothing was posted, pushed or merged.

Underworld development team with AI support from Claude Code

[lmoresi]

Response commit a9fd1d16.

MB-1 — sign reconciled, and no golden moves. The geometric normal keeps the corrected sign ("outward" = out of the domain, which on an inner arc points toward the centre); boundary_normal_traction("Lower") on an annulus goes −5.233110e-02 → +5.233110e-02, magnitudes identical to every digit, with dynamic_topography_field and a prescribed datum following. Convex boundaries untouched.

The analytic override is deliberately not reoriented: normal= means "use exactly this direction", and auto-orienting would silently change the meaning of a user's prescribed datum. So X/|X| on an inner arc stays inward-of-domain and now disagrees in sign with the default, on purpose — the docs say so and tell users to pass -X/|X| for the domain convention.

Consumer audit: every curved-boundary golden already goes through an explicit normal=, including the CMB values (test_1018:699-704 = 0.77060, test_1064 GOLDEN_SPHERICAL3D_TOPO), and those already pass normal=-unit_r on Lower — the direction the new default now produces. The two geometric-path goldens are on a box Top and are sign-normalised in the assertion. Nothing to re-baseline. (Correcting round 1's audit: the geometric path is not straight-box-only — the new tests drive it on a concave annulus Lower; they passed regardless because |A z| is sign-blind.)

Pinned by test_geometric_normal_points_out_of_the_domain_on_a_concave_boundary, whose negative control asserts the frozen pre-#560 rule gave +1 on the inner arc — without that the assertion would pin nothing.

A defect the audit found in my own change, now fixed: the support-cell flip had no getSupportSize(f) == 1 guard. On an internal boundary both support cells exist and support[0] is arbitrary, so neighbouring facets of one surface could be oriented oppositely and cancel in the weighted sum. Latent (never exercised today) but introduced here; both sibling implementations guard the same way.

MB-2 — the parallel solver tests now run. scripts/test.sh --p N ran only test_075*; the test_10*py line was commented out, so test_1017 and test_1062–test_1068, including this fix's own guard, executed at no rank count. Enabled. test_1063 then fails for real — #495, 3.4% partition dependence in the multiplier path, which never calls _boundary_velocity_nodes — so it is xfail(strict=False) with the measured numbers in the reason, and the other five now run. My file renamed 1067 → 1068 to clear a collision with test_1067_newton_cold_start.

HIGH — _local_boundary_candidates claimed "support 1 still means the domain boundary"; it does not, the overlap fringe also has support 1 (13–27% of its output, measured). The docstring now says it over-collects by design and that the caller's w.any() zero-sum test is what discards the extras, marked load-bearing at both ends, with boundary_normal_traction named as the second consumer.

Verification: scripts/test.sh --p 2 end to end, exit 0 — 13 serial batches (1846 passed) plus the parallel phase (test_075* 38 passed 1 skipped; test_10* 25 passed, 2 xfailed), and the test_10* batch standalone at np=4: 25 passed, 2 xfailed. |A z|/|A|_F unchanged by the support guard: 6.624e-20 at np=1, 6.690e-20 at np=4.

Marked, not silently fixed: two stale copies of the old orientation rule carry TODO(BUG) — boundary_flux._node_normals (dead: its only caller guards it, but rank-local and now inconsistent with rotated_bc on concave boundaries) and FreeSurface._normal_direction (unconditionally +radial, would fight the now −radial σ_nn if an inner free surface were ever built; unreachable today).

Underworld development team with AI support from Claude Code

[lmoresi]

Response commit 564bcad8 — and the PR is exonerated by measurement, not by argument.

We built probe/parallel-batch-baseline: development at this PR's merge base with a single change, the scripts/test.sh line this PR uncomments, dispatched on the same CI. Result: 7 failures at the baseline, with numbers identical to every digit to what this branch shows. The five still failing here plus the two already xfailed are exactly those seven. This PR causes none of them; it makes them visible.

Filed as #564 — and it is wider than #495: the ROTATED path is partition-dependent too (annulus 1.7e-04, spherical 3-D 1.1e-03, spherical topography 5.9e-03, prescribed datum 9.4e-06), not just the constrained/multiplier path. Those rotated cases pass an explicit analytic normal= whose values are bit-identical before and after the nodal-normal fix, so the mechanism is excluded as well as the attribution.

All seven are now xfail(strict=False) against #564, each carrying its own measured serial-vs-np=2 numbers and the note that it reproduces at the merge base. The two #495 ones are re-pointed at #564 "of which #495 is one member".

The batch stays enabled, which is now the more valuable half of this PR: test_1017 and test_1062–test_1068 had executed at no rank count in CI since they were written, so the assertions that catch #564 were dormant in the tree. Twenty parallel solver tests now run on every CI run.

scripts/test.sh --p 2 end to end: exit 0, "All tests passed" — 13 serial batches (1846) plus parallel test_075* 38 passed 1 skipped and test_10* 20 passed, 2 xfailed, 5 xpassed; the same at np=4 independently. The 5 XPASS are exactly the cases #564 records as passing on macOS/arm64 — which is why strict=False is the correct marker, and why a local green must not be read as those being fixed.

Underworld development team with AI support from Claude Code