Free surface in 3D: spherical shell support, dimension-general surface machinery, and the 2D compatibility-floor fix

[lmoresi]

What this is

Workstream C first pass (Louis + Claude; @gthyagi for review): uw.systems.FreeSurface runs end to end on a 3D spherical shell, built by making the surface machinery dimension-general rather than porting it piecewise. On the way, the 2D strong-datum compatibility floor was fully resolved (three stacked causes, each measured before fixed).

The 2D fix (rides first, benefits everything)

Strong-datum consistent solves on a deforming ring stalled honestly at rel ~2e-3, residual 100% in the pressure rows (block-split evidence). Three causes, closed in order, each verified by rerunning the power-law acceptance:

1. Trace-mass gauge — the demean weights become the P1 lumped boundary mass accumulated per owned DMPlex facet (necessary; measured not dominant).
2. Directed flux strip (maintainer ruling) — the datum's mean volume flux through the deformed facet normals is stripped with the same FE surface integral the residual uses (Φ = ∮ũ_n(n̂·Γ)ds, S = ∮(n̂·Γ)ds, mesh.Gamma), so the constrained field carries exactly zero discrete net flux.
3. Pressure-gauge quotient convergence — the constant-pressure mode is projected out of every increment by the inner KSP (it is the gauge of an enclosed domain) and is not exactly a left-nullvector on a deformed faceted boundary; the outer loop now measures the quotient-space residual, mirroring PETSc's projected residual for singular systems.

Result: every step of the power-law annulus acceptance converges (3, 3, then 1 warm-started increment per step; previously honest stalls).

The 3D machinery

• Surface mean/demean: reduction over owned-facet partial trace-mass weights — each facet counted exactly once globally, two scalar allreduces, no gather, no ordering. Same code is the 2D line gauge and the 3D area gauge.
• σ_nn on the 3D P2 trace: new boundary_fluxmass="p1" — P1-projected recovery (edge-midpoint loads folded exactly onto vertices via the P1⊂P2 interpolation transpose; sound P1 lumped triangle mass; midpoints read back as the P1 interpolant). Chosen over mass="consistent" because the P2 vertex-integral checkerboard (Fix 3D boundary reaction flux recovery #404 hold) sits exactly at the vertices the P1 h∞ field consumes. The manager translates its lumped default to p1 on 3D meshes.
• 2D-only features refused per piece (ring Taubin filter, tangential transport) instead of the blanket dim != 2 guard; the ring's gather generalises via lexicographic coordinate keys for its remaining diagnostic uses.
• Interior carrier reads the deform displacement by direct nodal row map (no on-node point evaluation — the uw.function.evaluate returns wrong values for P1 fields at points exactly on cell edges (3D) #432 3D exposure class — and no per-step global location).

Evidence (~/+Simulations/FreeSurface/spherical_fs/)

Spherical Y20 topographic relaxation (constant-density shell, background_buoyancy="analytic", resolved timestep): clean exponential decay at rate 6.1–8.2 vs the half-space Cathles reference 12.5 — an O(1) shell/fixed-base correction in the physically correct direction — with the equilibrium modal bias falling 16.3% → 2.2% of A0 over one resolution step (matching the 2D production-resolution figure of ~2.4%; the #431-class discrete recovery defect, resolution-convergent, not structural).

Tests: new test_1072 (3D capability, 17s); full 2D battery green throughout (test_1018 18/18, 1070, 1071 serial; 1070 + 1066 at np2 and np4 for the quadrature/datum tests).

Deliberately delegated to the review/benchmarking pass (maintainer decision)

This PR ships the capability with first-evidence physics; the detailed benchmarking is left open, and this is the note recording that decision:

• exact analytic shell-rate comparison (degree-l relaxation, radius ratio 0.547, no-slip base) and a rate/floor convergence study;
• low-Ra spherical convection with bounded T;
• 3D parallel validation — a spherical-cap datum finally exercises true empty-boundary ranks (np2/np4 annulus partitions never do);
• the CBF sigma_nn recovery on deformed/faceted boundaries: ~2.4% amplitude deficit (chord-vs-arc boundary mass suspected) #431 chord-vs-arc boundary-mass question, which this work's floor evidence feeds directly.

Related: #457 (TI consistent tangent, found during workstream B) still limits Newton-rate claims for anisotropic rheology in any dimension.

An adversarial review follows as the first comment.

🤖 Generated with Claude Code

Underworld development team with AI support from Claude Code

[Copilot]

Copilot was unable to review this pull request because the user who requested the review has reached their quota limit.

[lmoresi]

Adversarial review (pre-review, per house practice)

Ordered by how much each should worry the reviewer.

1. The quotient-space convergence measure is doing load-bearing work. The 2D floor fix's third leg declares convergence while a physically-real (small) incompatibility remains in the pressure-gauge component. This is standard for projected singular systems, and the component is minimized by the directed flux strip — but nothing REPORTS its size per solve. A user whose deformed geometry develops a large gauge-mode residual would see "converged" with no hint. Worth a follow-up: log the projected component's magnitude when it exceeds, say, rtol·reference.

2. The P1-projected σ_nn recovery is validated by physics, not by an exact oracle. The fold-down is exact as a load transfer (interpolation transpose), but its recovery ACCURACY on curved 3D traces is only evidenced by the Y20 benchmark's behaviour (right direction, resolution-convergent bias matching the 2D figure). No manufactured-traction test exists for the 3D path the way the SolCx topography test pins 2D. The delegated benchmarking should include one.

3. 3D parallel is untested. All 3D evidence is serial. The surface-mean reduction and the lex-key gather are collective-by-construction, but the empty-boundary-rank path (a rank owning no piece of the surface) has never actually executed — annulus partitions at np2/np4 always give every rank boundary. A spherical shell at np4+ is the first genuine test, and it is explicitly in the delegated list.

4. The benchmark's fitted rate is not resolution-converged (8.2 at res 0.2, 6.1 at res 0.15, half-space reference 12.5): the O(1) shell correction claim is qualitative. The floor's clean 16%→2.2% convergence is the strong quantitative result; the rate needs the delegated study before any analytic comparison is asserted.

5. Smaller items. The 3D lex-key gather assumes the surface deforms only along the normal (ordering built once — same assumption as the 2D ring, inherited, documented). The mass="p1" mode returns midpoint values as P1 interpolants — correct for the P1 h∞ consumer, but a caller wanting genuine P2 pointwise recovery must use "consistent" and inherit the checkerboard risk; the error message says so. estimate_dt is CFL-based and can hop the entire relaxation in one L-stable step (correct but surprising — the benchmark README documents it; a rate-based dt suggestion is a possible ergonomic follow-up). The refused-in-3D features raise at construction, not import, so a 2D script switched to a 3D mesh fails at the right moment with a named alternative.

What is well-covered: the 2D regression surface (every prior test green throughout, serial and np2, including the re-derived quadrature invariant); the three-cause floor fix each verified by measurement rather than assertion; capability pinned by test_1072 at 17s.

Underworld development team with AI support from Claude Code

[lmoresi]

Adversarial review — head fc4e3c3 (post-rebase), isolated worktree, own build

The PR body deferred benchmarking to the review pass, so this review verified the five headline claims empirically as well as attacking the code. All five hold; two MAJOR defects found, both fixed in response commit 26eec2d.

Findings

1. MAJOR (fixed in 26eec2d) — free_surface.py_surface_node_coords called facet_is.getIndices() unguarded: a rank owning NO surface facets gets a NULL-handle IS and any method call on it segfaults the whole job. Reproduced 3/3 at np=2 and np=3 (1×2 box, rank without Top facets dies in __init__). Every other stratum-IS call site carries the guard — including this PR's own _surface_weights_local. The line predates the PR verbatim, but the PR's "parallel-safe by construction" gauge is unreachable past it under partition asymmetry. With the guard: np=2/np=3 solve/advance/deform runs complete, gauge invariants identical across partitionings (step-0 shape_max 1.851923e-4 serial == np2 == np3; gauge 2e-16). Verified post-fix: np=2 --with-mpi plume suite 7 passed.
2. MAJOR (fixed in 26eec2d) — the pressure-gauge quotient projection silently passed an INCOMPATIBLE datum: enclosed annulus, prescribed constant outflow 0.01 on Upper (net flux 0.0628, unabsorbable) → converged=True at rel 7.5e-9 while the realized net flux was 0.0627 — a silent mass-conservation violation where the unprojected loop failed loudly (rel 3.0e-1). FreeSurface itself is safe (it strips its own datum flux), but generic add_rotated_freeslip_bc users lost a loud failure. Fix: the projected component is stashed as solver._rotated_pressure_gauge_residual and a converged solve warns when it dominates the converged residual. Compatible-case control: fields bit-identical with/without projection (v 1.4e-17, mean-free p 2.1e-16) — the projection changes measurement only, as designed.
3. MINOR (fixed: assertion added) — _surface_weights_local's wmap.get(key, 0.0) fallback could drop trace-mass weights silently on any ring-vs-plex coordinate divergence; now asserts exact conservation of the accumulated mass.
4. MINOR (fixed: docs) — mass="p1" was undocumented at both recovery entry points; two docstrings still said "arc length" for what is now the FE trace-mass gauge.
5. MINOR (open, tracked) — test_1072 is capability-only by design: the mean-free assertion is true by construction and the shape bound survives a sign flip, so CI does not pin the Y20 rate. Reference number from this review for the follow-up test: fitted decay rate 0.0842 vs half-space Cathles 0.204 (ratio 0.41) at R_i/R_o=0.547, cell 0.35.

Attacks that failed (independently measured)

• Gauge claim: 3D shell weights sum == BdIntegral(1) to 9e-15 serial / 2e-15 np3; weighted mean vs FE integral O(h²)-consistent.
• P1 σ_nn fold: no midpoint-pairing bug; linear-field oracle + quadratic negative control behave as claimed (consistent P2 machine-exact on linear, O(h²) on quadratic; p1 ~O(h) at patch-edge vertices — the accepted trade for killing the vertex checkerboard).
• Rotated free-slip: the custom-FMG velocity sub-KSP must converge, not just apply #465/One owner for the geometric-MG option bundle; rotated free-slip picks up a mesh-owned hierarchy (#468, #467) #471 composition after the rebase: geometric-FMG rotated tests pass with the projection active, strict post-Rotated free-slip: the custom-FMG velocity sub-KSP must converge, not just apply #465 asserts hold (velocity pc_type mg, ksp_its ≤ 3, counters populated); pres_is defined before every rotated_residual path including LU and rest-state.
• Parallel: test_1070 MPI np2 passed; test_1072 under mpirun -n 2 passed; _surface_mean is two allreduces, no gather, no rank-conditional collectives; 3D refusals fire collectively at construction.
• Y20 rate measured at 0.41× half-space Cathles — the claimed O(1) shell correction, right direction.
• 2D regression: the only test_1070 edit STRENGTHENS the ring-weights assertion; no tolerances loosened anywhere in the diff.
• Full level_1 and tier_a gate: 554 passed, 0 failed.

The rebase itself (4 commits replayed over the #458/#465/#471/#456 wave) was also under review: only the CHANGELOG needed manual resolution, and the semantic composition inside solve_rotated_freeslip verified live.