Rotated free-slip: prescribed-normal datum through the SNES path; one unified Newton loop (#438, #403)

[lmoresi]

What this is

Workstream B of the free-surface roadmap (joint session with @lmoresi): the prescribed wall-normal datum u·n̂ = ũ_n now runs through the full nonlinear rotated machinery, and the separate linear/nonlinear rotated solve paths are unified into one Newton/Picard loop. Closes #438's acceptance checklist and #403 items 2 & 4.

Design

• Datum through Newton = feasible iterates. The constraint is affine, so accepted iterates carry v_n = ũ_n exactly (an affine snap in the rotated frame) and every increment satisfies the homogeneous constraint n̂·δ = 0 — the datum never touches the tangent, the increment RHS, the FMG prolongation, or the null-space handling.
• Cold starts take the affine lift, not the snap. Snapping a zero state onto a datum manufactures a boundary-strip strain state whose shear-thinning tangent produces a first Newton step ~4 orders too large (measured; no line-searched step descends). The first increment instead carries the datum jump via zeroRowsColumns(x̂, b̂) at the rest-state tangent — which is exactly the old linear path's operation sequence, and is what makes unification natural:
• One path. A linear model converges after that first increment; the loop self-terminates, so the up-front _residual_is_nonlinear probe leaves the dispatch (two Jacobian assemblies saved per linear rotated solve; the probe survives as a lazy guard in the picard+pure-Newton corner). Direct LU becomes a per-increment diagnostic option. One result-dict shape (reaction = F(u); the linear-only A/b keys retire). Linear solutions change at solver-tolerance level only (bit-preservation deliberately relaxed — maintainer decision in-session).
• API: add_rotated_freeslip_bc(conds, boundary, normal=) accepts the non-zero scalar datum (value-first, per the Charter); the FreeSurface manager uses it and its penalty-fallback catch is gone — consistent_constraint="strong" works for nonlinear rheologies.

Three latent bugs found and fixed (each was masking the next)

1. Rank-local datum branching — both paths branched on the rank-local datum_map; a rank owning no datum boundary nodes desyncs the collective sequences (np>1 deadlock class; latent in the linear path since the datum landed). The activity flag now rides a collective PETSc Vec norm of the lift vector.
2. Essential-pinned rotation modes admitted to the null space — _mode_satisfies_constraints checked only the rotated constraint rows, so a rigid-rotation mode pinned by an essential BC on another boundary was projected out of every increment RHS: Newton converged superlinearly to rel ~2e-5 and floored there, irreducibly. Modes are now also verified as null vectors of the assembled operator (‖J·m‖ against the velocity diagonal scale); after the fix the same solve reaches rel 2e-13 in 11 iterations and correctly leaves the gauge in place.
3. Step-norm false convergence — a tiny Newton step was reported as converged even at rel |F̂| = 0.85 (stiff tangents make tiny steps far from the solution). The exit is now verified against the reference scale max(r0, ‖F̂(0)‖), which also fixes the warm-start rtol trap (rtol relative to a good warm start's own small r0 demands ever-more absolute accuracy).

Plus: FreeSurface's consistent solve warm-starts from the free solve's converged fields (a cold start at the power-law regularisation floor stalls at rel 0.13 while the free solve converges beside it).

Evidence

Acceptance runs in ~/+Simulations/FreeSurface/annulus_fs_convection/b_snes_datum_acceptance/ (README carries the full table and attempt history):

• Power-law annulus FS convection, consistent_constraint="strong", matched steps: material-boundary error 5.3e-3 (strong-linear reference 7.8e-3; penalty 3.7e-2), net surface flux 1.4e-4 (penalty 2.0e-3).
• TI fault-bearing smoke (power-law background, weak-plane band ratio 0.01, cos θ datum): converges, 24 its to rel 6e-8, σ_nn recovery clean.
• Nonlinear cos θ datum imposed to 3e-16 through genuine Newton iteration (new regression in test_1018); nonlinear datum at np2/np4 (new test in test_1066).

Tests: test_1018 18/18 serial · test_1066 np2 + np4 · test_1070 serial + np2 (7) · test_1071 serial.

Known limitations / follow-ups (all recorded)

• Deformed-ring datum compatibility floor (rel ~2e-3, honestly reported unconverged): confirmed by block-split to be 100% pressure-row residual — the arc-length (trapezoid) demean is not FE-quadrature flux-free on a deformed ring; the momentum equations are at machine precision. Fix = trace-mass-weighted demean (also the 2D instance of workstream C's datum gauge). Surface metrics unaffected.
• TransverseIsotropicFlowModel: consistent (Newton) tangent inconsistent with the residual when anisotropy is active #457 (found during acceptance): the TI consistent tangent is inconsistent with its residual when anisotropy is active — native-path defect, rotated machinery exonerated; TI "Newton" runs at Picard pace. Separate constitutive-model work.
• The naive LU pressure pin stays serial-only (TODO(BUG) retained); _eval_datum inherits the uw.function.evaluate returns wrong values for P1 fields at points exactly on cell edges (3D) #432 3D point-evaluation exposure; np4 is suggestive rather than a guaranteed empty-rank exercise (annulus partitions give every rank boundary).

New subsystem documentation: docs/developer/subsystems/rotated-freeslip.md; changelog entry included.

An adversarial review follows as the first comment, per house practice. Draft until @lmoresi signs off.

🤖 Generated with Claude Code

Underworld development team with AI support from Claude Code

[lmoresi]

Adversarial review (pre-ready, per house practice)

Ordered by how much each item should worry a reviewer.

1. The unconditional lift-step acceptance has no safety net. On a cold start with a datum, iteration 0 is accepted without any line search. For the tested rheologies the lift lands on a smooth state, but a rheology whose rest-state tangent is a bad global model (e.g. strongly pressure-dependent viscosity, or yield already active at rest) could get a first iterate far uphill, and the loop would then need the line search to walk back from a worse place than zero. Mitigations: bounded by max_it + stall exit + honest reporting; no observed failure. But the reviewer should know this is argued-safe, not proven-safe.

2. converged=True semantics changed twice in one PR. The reference scale is now max(r0, ‖F̂(0)‖) and the step-norm exit verifies against it. This is more honest, but it means: (a) a warm-started solve can report converged at an absolute residual that is rtol·‖F̂(0)‖ even if its own start was much better — by design, but it is a contract change; (b) any downstream code that compared rnorm/rnorm0 to rtol will see different verdicts. I audited the in-repo consumers (_capture_rotated_report, FreeSurface) — clean — but out-of-repo scripts reading _rotated_freeslip_info will see behaviour change.

3. Linear rotated solves are no longer bit-identical to development. Deliberate (maintainer ruling in-session: within-tolerance is the acceptance), and test_1066's hardcoded ∫v·v reference still passes at 1e-6 — but that pass is luck of magnitude for the nullspace fix (the spurious projection's imprint on the linear isoviscous case happened to be < 1e-6). A different fixture could have needed a re-baseline; anyone bisecting a downstream tolerance-level change will land on this PR.

4. The occasional polish increment is a real cost on linear solves. When the first residual check sits at the tolerance boundary the loop pays a second (cheap, but nonzero) KSP solve. Measured: FS plume wall-time unchanged, test_1018 faster overall (probe removal wins). But a user with many marginal linear solves could see it; the mitigation (inner KSP slightly tighter than outer) was considered and deferred.

5. The nullspace mode verification depends on an assembled J. The build was moved inside the loop to guarantee that, but _finalize_rotated_solution's gauge-removal check can run with an unassembled J on one exit path (trivially-converged warm start before any Jacobian assembly this call) — it then skips the operator test (treat-as-pass, matching old behaviour). In practice J persists from a previous solve; still, that path is weaker than the main one.

6. The FreeSurface warm-start copies fields via .array every step. Correct and cheap relative to the solves, but it assumes the free solve converged — if the free solve diverged, the consistent solve now inherits its garbage rather than starting clean. The manager already warns on free-solve divergence; the coupling is new though.

7. Parallel coverage is honest but not exhaustive. np2/np4 pass for both datum tests; np4 on an annulus still gives every rank a boundary segment, so the true empty-datum-rank path (a 3D cap datum at high rank counts) is exercised only by construction (collective flags), not by test. The _eval_datum empty-rank + #432 3D exposure is documented, unresolved.

8. Deleted code / behaviour surface. The linear one-shot's removal takes with it: the Uhat result key (nobody consumed it), the internal nsp/custom_Pl self-build in _solve_rotated_iterative (callers audited: one), and the A/b fallback in boundary_normal_traction (dead by construction). If an external script called solve_rotated_freeslip_nonlinear by name, it breaks loudly (name gone) — accepted, the pyx dispatch is the only in-repo caller.

9. What we could not test. The 3D spherical datum (no 3D FreeSurface yet — the guard stands); TI Newton-rate claims (blocked by #457, found during this work); behaviour under _rotated_use_lu at np>1 (still guarded by the naive-pin TODO(BUG), serial-only as before).

Net: the riskiest items are (1) and (2), both semantic rather than correctness regressions, both documented in the subsystem page. The three fixed bugs each have a regression guard (deadlock class → np2/np4 tests; nullspace floor → the 1018 nonlinear datum test converges to tolerance, which the floor made impossible; false convergence → the TI smoke asserts the honest verdict).

Underworld development team with AI support from Claude Code