Free surface with nonlinear / anisotropic rheology: strong material-boundary constraint is linear-only (penalty fallback leaks)

[lmoresi]

Why this matters

The rotated-BC free surface (uw.systems.FreeSurface, PR #435) is the only working
free-surface strategy, and nearly every target application is nonlinear — with
anisotropic nonlinearity (TI faults) required for the flagship cases, potentially WITH
the free surface (faults reaching the surface). This must not be a blocker, and it is
not: everything runs today. But one leg runs on a fallback whose deficiency we have
now measured precisely, so it needs to be first-class work, recorded here.

Current compatibility matrix

combination	status
rotated free-slip (u.n=0) + nonlinear	WORKS (SNES-integrated rotated path, honours consistent_jacobian)
rotated free-slip + anisotropic (TI)	works at constraint level (rotation is rheology-blind); solver-config caveats only (Schur/mass pmat choice for TI)
FreeSurface free + held solves + nonlinear	WORKS (held auto-dispatches nonlinear rotated; sigma_nn recovery is constitutive-agnostic)
FreeSurface consistent solve (u.n = u~_n != 0) + nonlinear	GAP: linear-only. Nonlinear raises NotImplementedError with a clear message; fallback = consistent_constraint="penalty"
FMG hierarchy + nonlinear	works by construction (Galerkin RAP of whatever tangent is assembled; re-formed each Newton iteration; performance-only caveats for localised viscosity)

Why the penalty fallback is not acceptable long-term

The consistent solve exists to keep the surface a MATERIAL boundary. Measured on the
isoviscous annulus (PR #435): the penalty leaks ~4% net volume flux through the surface
and misses the prescribed rate by 6.4x more than the strong datum — the exact mechanism
that put semi-Lagrangian departure feet outside the domain and injected hot material
into the cold boundary layer. Nonlinear runs on the fallback inherit that class of
error at whatever amplitude their regime produces.

Required work

1. Extend the rotated prescribed-normal datum to the SNES path (Finish rotated free-slip prescribed-normal datum: wire FreeSurface consistent solve, pyx API, parallel, nonlinear #403 item 4): the
   linear path imposes the datum via zeroRowsColumns(x,b) with a solution-row restore;
   the Newton loop needs the equivalent per iteration — set the rotated-normal DOF of
   the iterate to the datum, zero the corresponding residual row, keep the tangent
   transparent (the machinery the u.n=0 nonlinear path already has, generalised to a
   non-zero, field-valued datum). Belongs with the pyx add_rotated_freeslip_bc datum
   argument (Finish rotated free-slip prescribed-normal datum: wire FreeSurface consistent solve, pyx API, parallel, nonlinear #403 item 2) and, ideally, with the planned linear/nonlinear rotated-path
   unification (the split itself is a design smell — flagged for the joint pyx session).
2. Anisotropic extras (small, config-level): a sensible scalar viscosity proxy for
   the Schur mass pmat under TI in the rotated solves; never-ksponly guidance carried
   over to the FreeSurface-derived solves.
3. Fault-at-the-surface interaction (open question, not solver work): the surface
   ring filter and tangential transport act along a ring that a TI weak zone may
   intersect; whether the filter should be fault-aware is a physics/design question to
   revisit once 1 exists and a fault-to-surface model runs.

Acceptance

Power-law annulus free-surface convection with consistent_constraint="strong" running
through the SNES path: material-boundary error and net surface flux at the strong-datum
level (not the penalty level) at matched steps; a TI fault-bearing smoke test compiling
and converging with the FreeSurface manager.

Related: #403 (items 2, 4), #435 (the 2D linear-complete toolkit), #423, #404 (3D).

Underworld development team with AI support from Claude Code

[lmoresi]

Rotated DOFs x anisotropy — decomposition of the concern (Louis, 2026-07-27)

Not SNES-specific, and worth separating into three different kinds of question:

1. Constraint algebra: no issue, by construction. The rotation is a per-node
orthogonal similarity on the ASSEMBLED operator (A_hat = Q A Q^T); it never touches the
constitutive law and imposes u.n = datum exactly for any tangent — symmetric TI viscous
or non-symmetric plasticity Jacobian alike. The reaction read-out is the assembled
residual: constitutive-agnostic. TI viscous anisotropy is linear and runs through the
existing ksponly rotated path today.

2. Formulation (the real worry): yield consistency at strongly-constrained nodes.
A strong constraint applies whatever reaction holds u.n — unbounded by material
strength. For a plastic material at the boundary the deviatoric part of that reaction
enters the yield criterion, so the constraint can demand stress states the material
cannot sustain; the nonlinear iteration then negotiates a hard kinematic condition
against a bounded-strength rheology at the same nodes (softening, reaction
redistribution, possible boundary-parallel localisation). Generic to strong Dirichlet
BCs on plastic materials — but sharper here because the constraint reaction IS the
physics we consume (sigma_nn -> h_inf, dynamic topography). Question to settle before
trusting fault-at-surface topography: what does the recovered sigma_nn of a YIELDING
constrained boundary mean, and does the free-surface update remain well-posed there?
(Candidate mitigations if it bites: consistent_constraint="penalty" locally near
yielding boundary segments; or a bounded-reaction variant of the constraint.)

3. Preconditioning (performance only): scalar viscosity proxy for the Schur mass
under a TI tensor; MG smoother robustness on anisotropic stencils. Already items here.

Cheap discriminating pre-test (no SNES work needed): linear TI viscosity with a
fixed weak-plane orientation + rotated free-slip on the annulus vs a Nitsche/penalty
reference — agreement on the flow AND the recovered sigma_nn validates the
rotation-x-anisotropy mechanics in isolation, leaving the joint session only the real
(formulation + datum-path) questions.

Underworld development team with AI support from Claude Code

[lmoresi]

Pre-test results: rotated x anisotropy mechanics EXONERATED; weak enforcement degrades

Linear TI (fixed 30-deg director), annulus, rotated strong free-slip (direct LU — no
preconditioner in the experiment) vs Nitsche weak reference (as implemented, gamma=10,
scalar-mu scaling). Leak measured at the machinery's own constrained DOFs (by label —
radius masks gave two successive false alarms and are banned from boundary metrics).

eta1/eta0	rotated constraint leak	sigma_nn amplitude	sigma_nn roughness	Nitsche leak (rel)
1.0 (iso)	2.0e-17	6.14	0.098	1.6e-2 (1.0e-1)
0.1	1.5e-16	5.51	0.142	3.4e-2
0.01	8.8e-16	5.12 -> 5.20 (res 0.06 -> 0.04)	0.451 -> 0.320 with resolution	0.37 (1.3e-1)

Findings:

1. The rotated constraint is exact under anisotropy — machine-precision u.n at 100x
   viscosity contrast, as the algebra (orthogonal similarity on the assembled operator)
   requires. The sigma_nn reaction is CONVERGENT under strong anisotropy: amplitude
   stable to 1.5% across resolutions, roughness falling with h — under-resolved
   anisotropic boundary stress, not a checkerboard. The mechanics half of this issue's
   concern is retired with numbers.

2. Weak (Nitsche) enforcement, as currently implemented, degrades badly with
   anisotropy: normal leak grows to 13% relative at eta1/eta0 = 0.01 (the penalty
   scales with a scalar viscosity that under-penalises when the normal-direction
   stiffness collapses). The large rotated-vs-Nitsche velocity gaps in the scan are the
   Nitsche solution being wrong. Caveat: a direction-aware penalty scaling might
   recover it — untested; "as implemented" stands.

3. Consequence for this issue's priority: the penalty fallback for the FreeSurface
   consistent solve will leak WORSE than its isoviscous 4% in anisotropic runs —
   exactly the target applications. The SNES-path strong datum (item 1 here) is
   therefore more urgent, not less. Remaining open question for the joint session is
   unchanged: yield consistency at strongly-constrained nodes (formulation).

Louis: "this is why I believed that we needed rotated dofs in the first place."

Script: ~/+Simulations/FreeSurface/annulus_fs_convection/teaching/ti_rotated_pretest.py

Underworld development team with AI support from Claude Code