Fix 3D boundary reaction flux recovery

[gthyagi]

Summary

• recover pointwise 3D boundary traction by dividing rotated free-slip nodal reactions by an assembled lumped surface mass
• support triangular P1 traces and P2 traces through four-triangle facet subdivision
• remove the traction gauge with a surface-mass-weighted mean
• add focused serial and MPI regressions using the Zhong et al. (2008) degree-2 internal-load shell

Root cause

The previous 3D branch in _desmear() returned the raw nodal reaction after subtracting its arithmetic mean. A finite-element nodal reaction is an integrated force, not a pointwise traction. The recovery must solve

M_boundary * sigma_nn = R

or, for lumped mass,

sigma_nn[i] = R[i] / m_boundary[i].

Without the surface-mass division, recovered 3D topography was mesh-dependent and nearly zero. The arithmetic mean also did not represent the surface integral gauge.

Implementation

For 3D triangular boundaries, the recovery now:

1. gathers and de-duplicates boundary topology and nodal reactions by coordinate;
2. uses each P1 triangle directly, or subdivides a complete P2 triangle into four linear triangles;
3. adds one third of each subtriangle area to its nodes;
4. divides the assembled reaction by the lumped nodal mass; and
5. removes the mass-weighted surface mean when requested.

Unsupported 3D mass choices and element traces fail explicitly. The existing 2D implementation is unchanged.

Validation

Zhong et al. degree-2 isoviscous internal-load response, ri=0.55, rint=0.775, ro=1.0, P2/P1:

Run	Surface topography	CMB topography
Zhong et al. reference	0.41920	0.77060
UW3 serial	0.4205055	0.7704049
UW3, 2 MPI ranks	0.4205056	0.7704048

Focused checks:

• serial 3D Zhong reaction/topography regression: passed
• 2-rank 3D partition-independence regression: passed
• existing 2D SolCx normal-traction regression: passed
• existing 2D dynamic-topography field regression: passed
• existing scalar boundary-flux regression: passed

This PR intentionally excludes the higher-level spherical-shell geoid and response postprocessing work. It contains only the boundary-flux correction and focused tests.

Context: #248

[lmoresi]

Adversarial review (public record) — strong core, but the P2 recovery has a structural defect

Thanks @gthyagi — the approach is right and most of it verifies beautifully. Full skeptic pass with live probes at np1/np2/np4 (PR branch merged onto current development, dedicated build):

What HELD (verified live)

• The partial_reaction contract and seam assembly: boundary_flux passes raw per-rank residuals (SUM correct), boundary_normal_traction passes assembled reactions (overwrite correct). Exact discrete check: P1 tet-box conduction gives flux ≡ −1.000000 at every node at np1/np2/np4 — a double-count would print −2.
• np4 partition independence (the PR claimed np2): Zhong DIAG values match serial to 1.5e-8/6.3e-8. Better than the golden's tolerance.
• The lumped mass geometry is exact: sum(lumped_mass) matches the PETSc faceted sphere area to 1e-11, P1 and P2.
• Coordinate-key dedup: −0.0/+0.0 keys collide correctly (Python equality), cross-rank coords are bit-identical (DMPlex-communicated; midpoints are 2-term IEEE-commutative means). Sound.
• Empty-rank safety: ranks owning no Lower facets run clean at np4; empty xs returns empty.
• The P1 path is genuinely correct — exact uniform-flux recovery.

The blocker: P2 pointwise recovery is structurally wrong, and the Zhong test cannot see it

On a straight P2 triangle, ∫φ_vertex dA = 0 exactly — so the consistent nodal reaction at every VERTEX is ≈0 for smooth flux, and no positive lumped mass can recover a signal from a zero numerator. Live, exact evidence:

• Tet box, true uniform flux 1, P2: recovered flux is 0.0000 at every vertex and −4/3 at every midpoint (exactly, serial and np2/np4).
• Your own Zhong config, harmonic fit split by node class: vertices-only → 0.0108, midpoints-only → 0.5580, all-nodes → 0.4205 (the golden). The all-node fit matches literature only because sphere triangulations have E≈3V, so (4/3)·E/(V+E) ≈ 1 — an exact structural cancellation in the fitted coefficient, not a correct pointwise field.
• dynamic_topography on a P2 field: vertices sit on a flat gauge shelf, midpoints carry 4/3-inflated signal — a surface-wide checkerboard with ~0.4 jumps. Anything pointwise (free-surface driving, plots, local σ_nn) gets garbage even though the l=2 coefficient looks perfect.

Asks before merge: (a) a consistent surface-mass solve on the P2 trace (M_boundary σ = R), or an explicit P1-only restriction with a clear raise; (b) a pointwise exact-flux box regression — P1 passes it today, P2 fails it today, and it's the test the harmonic fit can't substitute for; (c) docstring updates — mass="consistent" is still advertised but raises in 3D, 3D hex meshes now hard-raise (release-note worthy: 3D box Nusselt-via-CBF was previously wrong-but-nonzero, now an error — arguably an improvement, but say so), and "partition-independent" should be qualified analytic normals only (with the default geometric facet normal the same diagnostic shifts 1.6e-4/9.1e-4 between np1 and np2 — pre-existing seam behaviour in the normal accumulation, not this PR's doing).

Filed separately (pre-existing, exposed by the probes — not blockers here)

• _assemble_volume_reaction produces garbage on inhomogeneous Dirichlet boundaries (g≠0), in 2D as well as 3D — sumR = −4.41 on a boundary whose true flux is ±1; existing tests pass only because they measure g=0 boundaries. → new issue below.
• The allgather materializes the full global boundary on every rank (~2 GB/rank for a 1M-node P2 surface, ×3 per component for full traction). Fine for coarse post-processing, needs a ticket before production 3D. → new issue below.

The core of this PR — P1 de-smear, global assembly, gauge — is exactly the missing 3D piece and worth landing promptly once the P2 story is honest.

Underworld development team with AI support from Claude Code

[gthyagi]

@lmoresi Thank you for the adversarial review. I reproduced the P2 vertex/midpoint defect and pushed bc62b870, which replaces the four-subtriangle P2 lumping with the exact consistent six-node triangular surface-mass solve.

Changes

• P1 triangular traces support lumped and consistent recovery.
• P2 triangular traces use M_boundary * sigma_nn = R with the exact 6 x 6 element mass.
• mass="auto" is now the public default: P1/2D retain lumped recovery, while 3D P2 selects consistent recovery.
• Explicit mass="lumped" on a 3D P2 trace raises because the vertex row sums are exactly zero.
• Unsupported non-triangular and higher-order 3D traces raise explicitly.
• Mean removal uses the assembled surface-mass weights.
• Docstrings and the developer changelog now describe the supported trace/mass combinations, non-triangular limitation, and analytic-normal qualification.

Pointwise box regression

Unit-cube conduction, exact outward Bottom flux -1, measured on the homogeneous Dirichlet wall:

Trace	Recovery	Serial maximum nodal error	np=2	np=4
P1	lumped	1.23e-9	pass	pass
P2	auto -> consistent	4.17e-9	pass	pass

The P2 assertion covers every vertex and edge midpoint. A separate regression verifies that explicit P2 lumping raises.

Zhong degree-2 validation

ri=0.55, rint=0.775, ro=1.0, P2/P1:

Sample class	Surface	CMB
Zhong et al. reference	0.41920	0.77060
UW3 all P2 nodes	0.41497	0.79322
UW3 vertices only	0.39523	0.84260
UW3 edge midpoints only	0.42159	0.77624

The regression now checks all-node, vertex-only, and midpoint-only coefficients independently, so the previous node-count cancellation cannot pass. Serial, np=2, and np=4 results agree within 1e-6 relative tolerance.

Focused local validation:

• 8 serial tests passed.
• 4 tests per rank passed at np=2.
• 4 tests per rank passed at np=4.
• Existing 2D SolCx traction, lumped/consistent behavior, dynamic-topography field, and scalar boundary-flux tests remain green.

Issues #407 and #408 remain separate as requested: the Zhong rotated-constraint reaction is not affected by #407, and this correctness patch retains the existing global topology reconciliation tracked by #408.

[gthyagi]

@lmoresi Follow-up resolution sweep for the exact P2 surface-mass recovery in bc62b870. Commit 1533574a parameterizes the diagnostic cell size while retaining cellSize=0.25 as the inexpensive CI regression default.

This is the no-self-gravity Zhong l=2 internal-load case (ri=0.55, rint=0.775, ro=1.0, P2/P1). Here 1/4, 1/8, and 1/16 denote the target cellSize, not a structured element count.

Target cellSize	Ranks	Surface, all P2 nodes	Error vs 0.41920	Surface vertex / midpoint	CMB, all P2 nodes	Error vs 0.77060	CMB vertex / midpoint	Wall time (s)
Zhong et al.	-	0.41920	-	-	0.77060	-	-	-
1/4	1	0.41497	-1.01%	0.39523 / 0.42159	0.79322	+2.94%	0.84260 / 0.77624	10.33
1/4	4	0.41497	-1.01%	0.39523 / 0.42159	0.79322	+2.94%	0.84260 / 0.77624	9.16
1/8	1	0.40717	-2.87%	0.37102 / 0.41924	0.80059	+3.89%	0.87703 / 0.77495	36.21
1/8	4	0.40718	-2.87%	0.37104 / 0.41924	0.80055	+3.89%	0.87701 / 0.77491	16.85
1/16	1	0.41404	-1.23%	0.39851 / 0.41921	0.78579	+1.97%	0.83037 / 0.77090	183.24
1/16	4	0.41404	-1.23%	0.39851 / 0.41921	0.78578	+1.97%	0.83036 / 0.77090	59.96

Observations:

• Serial and four-rank coefficients agree at every resolution. The largest relative serial/MPI difference among the six reported node-class coefficients is 0.0055% (1/8); at 1/16 it is 0.00056%.
• The P2 midpoint fits converge closely to the paper: at 1/16, surface is 0.41921 and CMB is 0.77090.
• The all-node fit is not monotonic over these three unstructured meshes. Its 1/16 errors are -1.23% at the surface and +1.97% at the CMB. The vertex/midpoint spread is smaller at 1/16 than at 1/8, but has not vanished.
• The exact uniform-flux P2 regression remains the pointwise correctness check for the mass solve; this Zhong sweep separately exposes resolution and node-class behavior rather than allowing the old node-count cancellation to pass unnoticed.

Validation after parameterization:

• targeted serial regression: passed;
• targeted four-rank regression with --with-mpi: passed;
• 1/16 serial peak RSS: 3.10 GB, no swaps;
• 1/16 four-rank wall time: 59.96 s versus 183.24 s serial (3.06x).

[lmoresi]

Adversarial re-review — round 2 (public record): the P2 fix is verified. Approving.

@gthyagi — the consistent trace-mass rework is exactly right. Full round-2 findings:

The round-1 defect is genuinely fixed

• The hard-coded 6x6 P2 (and 3x3 P1) triangle surface mass matrices were verified entry-for-entry against exact symbolic integration on the reference triangle, including the −4 opposite-vertex/midside pairings for your node ordering.
• The decisive probe from round 1 — unit flux through a tet-box wall, on an IRREGULAR mesh — now recovers ±1 pointwise at vertices AND midpoints to ~1e−11, serial, np2 and np4 (round 1: exactly 0 at vertices, −4/3 at midpoints).
• mass="auto" semantics verified: consistent for 3D-P2, lumped elsewhere; explicit lumped-on-P2 raises loudly; P3/quad traces refuse legibly.
• Partition independence holds: all six Zhong coefficients agree np1/np2/np4 to ~6e−7 and match the new goldens.
• Your own test batteries all pass in the round-2 environment (test_1018 15 passed incl. the new 6-way vertex/midpoint assertions; test_1064 np2 8 passed; test_1015/1019/1065 regression sweep green).

One genuine subtlety found — filed as #414, not a blocker

On CURVED boundaries the P2 vertex values converge only slowly: the vertex basis has zero surface mean, so vertex reactions carry only the O(h) facet-geometry error, and the (correct) consistent solve faithfully reconstructs that error. Measured on a spherical conduction probe: vertex error 93%→55% under one refinement while midpoints go 2.8%→0.7% (superconvergent). Your Zhong vertex-only fit drifts away from the reference with refinement (0.395→0.371) while the midpoint fit converges onto 0.419 — the wide test tolerances absorb this, which is fine, but consumers of pointwise σ_nn on curved boundaries should read midpoints or fitted quantities. A docstring caveat is the ask (#414).

Housekeeping from the round

• The review's independently-rediscovered 'flux garbage on g≠0 Dirichlet boundaries' is CBF boundary_flux: _assemble_volume_reaction is wrong on inhomogeneous Dirichlet boundaries (2D and 3D) #407, already fixed on development by PR CBF flux on driven boundaries, kdtree module identity, global velocity_error #412 (the reviewer's base predated that merge) — your branch will pick the fix up on merge; nothing to do here.
• The consistent solve's redundant global assembly extends the scaling debt tracked in boundary_flux 3D assembly gathers the full global boundary onto every rank (memory scaling) #408 (noted there).
• 2D P1 traces still KeyError in the old 2D lumped path (pre-existing) — filed as boundary_flux: 2D P1 traces crash with a bare KeyError in the lumped line-mass assembly #413 since the new docstring enumerates 2D behaviour.
• A stale pr404-head tag from round 1 was shadowing branch names in local reviews — deleted.

Merging. Thanks for the fast, correct turnaround — the four-subtriangle → consistent-mass pivot was exactly the right call.

Underworld development team with AI support from Claude Code