CBF sigma_nn recovery on deformed/faceted boundaries: ~2.4% amplitude deficit (chord-vs-arc boundary mass suspected)

[lmoresi]

CURRENT SCOPE (2026-07-27) — the report below is the ORIGINAL filing and its
headline mechanism is retracted (probe artifact: radius-mask surface selection on a
deformed ring; see the correction comment below). The recovery is exonerated at ~0.976
of exact. What this issue now tracks is only:

The CBF σ_nn recovery on a deformed, faceted boundary under-reports amplitude by
~2.4% (load-only 0.976, driving 0.972 on a 5%-deformed ring at cellSize 0.06).
Chord-vs-arc boundary mass on the faceted ring is the suspected bulk. Likely related to
the #414 curved-boundary error family. Next step before any fix: measure how the
deficit scales with deformation amplitude and resolution. The free-surface pathology
originally attributed here was a sign error in the FreeSurface manager, fixed in PR #435.

---

Original report (headline mechanism retracted — kept for the record)

Summary

On a deformed boundary, the CBF reaction recovery behind dynamic_topography /
boundary_normal_traction returns only ~0.16 of the current-shape load component
rho_0gh of the constraint reaction, while recovering the dynamic (driving) component
exactly. Any sigma_nn measurement on a deformed boundary under a full-density body force
is biased by this; in the FreeSurface manager it made the full-Boussinesq scheme park the
surface at 24% of isostasy with steady O(vrms) flow through a stationary surface.

Decisive probes (annulus r_in=0.5, r_out=1.0, cellSize=0.06, rotated free-slip lid)

A — pure statics, the recovery has nothing to hide behind: T=0, body force
-rho_0grhat (rho_0g=1e6), mode-4 boundary deformation amplitude 0.05. The discrete
solve is exact hydrostatics (u=0, p'=const); the constraint reaction is rho_0g*h(theta)
to solver precision. Recovered (mass="lumped"): m=4 amplitude 0.158 of exact, phase
correct (-179.9 deg for the -(sigma-mean)/rho_g convention).

B — the dynamic leg is clean: driving-only vs driving+full-background on the
UNDEFORMED mesh (self-load = 0): recovered m=4 identical to 4 decimals (ratio 1.0000).
The huge hydrostatic RHS does not corrupt the recovery; the loss is specific to the
deformed-boundary load.

Suspected mechanism (unconfirmed)

The de-smear credits the constrained (rotated-normal) rows; on an irregular faceted
deformed ring part of the hydrostatic boundary term assembles into rows the recovery
does not read (tangential-rotation split and/or P2 edge-node rows). Needs a source-level
pass through boundary_flux / the reaction extraction.

Workaround in place (997a6b0)

FreeSurface(background_buoyancy=...) runs a twin held solve with the background force
alone on the same geometry and converts the DIFFERENCE of the recoveries — the defective
load leg cancels identically. Validated: relaxation decays 264 km -> exactly 0; the
parked equilibrium restarts to 97% of the isostatic referee in 25 steps with the
through-flow collapsing 30x in one step.

Why fixing it properly matters

The load is analytic (rho_0gh from the geometry), so with a faithful recovery
h_inf = h + (sigma_rec - rho_0gh)/(rho_0*g) needs NO extra solve per step — the twin
solve is a stopgap. Fixing kappa -> 1 also un-biases sigma_nn / dynamic-topography
measurements on deformed boundaries generally (they are a primary model OUTPUT).

Probes: ~/+Simulations/FreeSurface/annulus_fs_convection/teaching/ (session scripts,
reproducible in ~2 min serial each).

Underworld development team with AI support from Claude Code

[lmoresi]

CORRECTION — the recovery is exonerated; this issue's mechanism was wrong

The 0.158 measurement in the report above was an artifact of the probe itself: its
surface selection used a radius mask (r > r_max - 0.03), which on a 5%-deformed ring
keeps only the CREST arcs — an m=4 Fourier fit over crest-only nodes is meaningless.

Redone with label-based ring selection (full ring, boundary velocity nodes):

• load-only on the 5%-deformed ring: recovered/exact = 0.976
• driving, deformed vs undeformed: 0.972
• full density: exact superposition of the two (0.0580 @ -32.6 vs predicted 0.0579 @ -32.5)
• raw nodal reaction n.(A u - b) vs sigma*mass: 0.986, uniform across vertices and midpoints

The actual defect was in the FreeSurface manager, not the recovery: under full
density the recovered field contains the self-load component +h_current, and the
manager's reduced-formulation negation flips it (h_inf = -h + driving/rho_g: fixed
point at HALF the equilibrium; the -1 eigenvalue is the period-2 ringing seen in every
full-density run). Fixed with NO extra solve in 9896114: background_buoyancy="analytic"
subtracts the geometric current height before the standard conversion. Validated
bit-identical (at reported precision) to the exact two-reaction reference on both the
relaxation benchmark and the equilibrium restart.

Remaining scope for this issue (small): the recovery's residual inaccuracy on a
deformed ring is ~2.4% (0.976) — worth understanding (chord-vs-arc mass on the faceted
ring is the likely bulk of it) but no longer implicated in any instability.

Lesson recorded: never fit boundary harmonics through a radius mask on a deformed
mesh — select by boundary label.

Underworld development team with AI support from Claude Code

[lmoresi]

Final scope after PR #435: the free-surface pathology originally attributed here is fully resolved in that PR (one-sign self-load fix; this issue's original mechanism retracted above). What remains here is only the ~2.4% residual inaccuracy of the recovery on a deformed ring (chord-vs-arc boundary mass is the likely bulk). Low priority; measure its scaling with deformation/resolution before deciding whether it warrants a fix.

Underworld development team with AI support from Claude Code