Rotated free-slip: the custom-FMG velocity sub-KSP must converge, not just apply

CLAUDE.md

@@ -326,7 +326,7 @@ The PETSc-based solvers are carefully optimized and validated. **NO CHANGES with
 
 ## Boundary Conditions: Free-slip
 
-**Prefer rotated strong free-slip** (`solver.add_rotated_freeslip_bc(conds, boundary, normal=None)`, value-first: `conds=0` is the only implemented datum)
+**Prefer rotated strong free-slip** (`solver.add_rotated_freeslip_bc(conds, boundary, normal=None)`, value-first: `conds=0` is free-slip; a non-zero scalar/expression `conds` prescribes the wall-normal datum `u·n̂ = ũ_n` strongly)
 to impose `v·n̂ = 0`:
 
 - Enforces zero wall-normal flow to **machine precision** (Nitsche / penalty leak ~1e-3).

docs/developer/CHANGELOG.md

@@ -6,6 +6,55 @@ This log tracks significant development work at a conceptual level, suitable for
 
 ## 2026 Q3 (July – September)
 
+### One Rotated Free-Slip Path, Now With a Prescribed Wall-Normal Velocity (July 2026)
+
+**Rotated strong free-slip now takes a prescribed wall-normal velocity datum
+(`add_rotated_freeslip_bc(conds, boundary, ...)` with non-zero `conds`) through
+the full nonlinear Newton machinery, and the separate linear and nonlinear
+solve paths have been unified into one** (#438; #403 items 2 and 4).
+
+The rotated constraint `u·n̂ = ũ_n` is the primitive behind both the held free-slip
+lid and the free-surface material-boundary condition — they differ only in the
+constraint right-hand side. Previously the non-zero datum existed only on a
+linear one-shot path, so a power-law or anisotropic rheology silently fell back
+to a weak penalty (which leaks worst exactly where anisotropy makes it matter).
+Now every rotated solve runs a single Newton/Picard loop in which accepted
+iterates carry the datum exactly; a cold start imposes it through the first
+increment's affine lift at the rest-state tangent (snapping a zero state onto a
+datum creates a boundary strain state a shear-thinning tangent cannot recover
+from — measured, not assumed); a linear model simply converges after that first
+increment, so the up-front nonlinearity probe is gone from the dispatch and
+every linear rotated solve saves two Jacobian assemblies.
+
+Three latent solver defects were found and fixed by making the loop report
+honestly along the way:
+
+- Branching on rank-local datum bookkeeping desynchronised the ranks'
+  collective sequences (an np>1 deadlock class); the datum-activity decision is
+  now a collective PETSc reduction.
+- A rigid-rotation mode pinned by an essential condition on *another* boundary
+  could still enter the solver null space (only the rotated rows were checked),
+  silently projecting an irreducible component out of every increment — Newton
+  converged superlinearly and then floored, far above tolerance. Candidate
+  modes are now verified as null vectors of the assembled operator, which
+  catches any form of pinning.
+- A tiny Newton step was reported as convergence even when the residual was
+  still large (a stiff tangent also produces tiny steps); the step-norm exit is
+  now verified against the problem's rest-state residual scale, which is also
+  the convergence reference for warm starts (rtol relative to a good warm
+  start's own small initial residual demands ever-more absolute accuracy).
+
+The free-surface manager's `consistent_constraint="strong"` therefore works for
+nonlinear rheologies: the penalty fallback is removed, and the consistent solve
+warm-starts from the free solve's converged fields. Acceptance: power-law
+annulus free-surface convection holds the material boundary at the strong-datum
+level (5e-3, versus 4e-2 for the penalty) with net surface flux 1e-4; a
+transversely isotropic fault-bearing smoke test converges through the same
+path. Direct LU per Newton increment remains available as a serial,
+preconditioner-free diagnostic (`solver._rotated_use_lu`).
+
+New subsystem documentation: `subsystems/rotated-freeslip.md`.
+
 ### Local Interpolation That Reproduces Linear Fields (July 2026)
 
 **The local scattered-point interpolator now has a linear-reproduction

docs/developer/index.md

@@ -36,6 +36,7 @@ same topic are reference or historical material subordinate to the governing doc
 | Coding style & API conventions | [UW3 Style Charter](UW3_STYLE_CHARTER.md) (detailed reference: [Style Guide](UW3_Style_and_Patterns_Guide.md)) |
 | Data access | [subsystems/data-access.md](subsystems/data-access.md) (internals reference: [NDArray System](UW3_Developers_NDArrays.md)) |
 | Local scattered-point interpolation | [subsystems/interpolation.md](subsystems/interpolation.md) |
+| Rotated free-slip & wall-normal datum | [subsystems/rotated-freeslip.md](subsystems/rotated-freeslip.md) |
 | Units | [design/UNITS_SIMPLIFIED_DESIGN_2025-11.md](design/UNITS_SIMPLIFIED_DESIGN_2025-11.md) |
 | Testing tiers | [TESTING-RELIABILITY-SYSTEM.md](TESTING-RELIABILITY-SYSTEM.md) |
 | Branching & releases | [guides/branching-strategy.md](guides/branching-strategy.md) |
@@ -186,6 +187,7 @@ subsystems/mesh-shape-relaxation
 subsystems/discretisation
 subsystems/solvers
 subsystems/boundary-stress-and-projection-postprocessing
+subsystems/rotated-freeslip
 subsystems/petsc-jacobian-layout
 subsystems/constitutive-models
 subsystems/constitutive-models-theory

docs/developer/subsystems/rotated-freeslip.md

@@ -0,0 +1,178 @@
+# Rotated strong free-slip and the prescribed-normal datum
+
+The rotated free-slip boundary condition imposes
+
+$$\mathbf{u}\cdot\hat{\mathbf{n}} = \tilde u_n$$
+
+strongly — as an exact per-node constraint, not a weak (Nitsche/penalty) term —
+by rotating each boundary node's velocity components into a (normal,
+tangential) frame and constraining the rotated normal component. `conds = 0`
+(the default datum) is free-slip; a non-zero scalar `conds` prescribes the
+wall-normal velocity, which is what the free-surface manager uses to keep the
+surface a material boundary.
+
+```python
+stokes.add_rotated_freeslip_bc(0, "Upper", normal=nhat)          # free-slip
+stokes.add_rotated_freeslip_bc(h_dot.sym[0], "Upper", normal=nhat)  # u·n̂ = field
+```
+
+`normal=None` uses the geometric facet normal; a sympy `1×dim` matrix in
+`mesh.X` supplies an analytic normal (exact `X/|X|` on curved boundaries — the
+preferred choice there); a constant array is also accepted. The datum must be a
+*scalar* (a number, an expression of `mesh.X`, or a scalar field read); on an
+enclosed boundary it must be discretely flux-free for incompressibility. A
+corner or 3D-edge node shared between rotated boundaries has no single normal
+and stays at the free-slip pinning (the datum is ignored there).
+
+Why strong rather than Nitsche/penalty: the constraint holds to machine
+precision (a penalty leaks ~1e-3, and the leak grows exactly where anisotropy
+makes the boundary condition matter), it is correct on curved/tilted/deformed
+boundaries, and the constraint **reaction is the boundary normal traction**
+σ_nn (`solver.boundary_normal_traction`, `solver.dynamic_topography`) — the
+quantity the free-surface machinery consumes. Reserve Nitsche for conditions
+that must morph in time (Dirichlet→Neumann ramps).
+
+**Implementation**: `src/underworld3/utilities/rotated_bc.py`; registration and
+dispatch in `petsc_generic_snes_solvers.pyx` (`add_rotated_freeslip_bc`).
+
+## One solve path
+
+There is a single driver, `solve_rotated_freeslip`: a manual outer
+Newton/Picard loop that rotates the residual and tangent every iteration
+(`F̂ = Q F`, `Ĵ = Q J Qᵀ`), imposes the constraint on the rotated normal rows,
+and solves each increment with a self-contained fieldsplit-Schur KSP (geometric
+FMG on the custom prolongation when a hierarchy is registered, else GAMG; the
+native 1/μ pressure mass as the Schur preconditioner). There is no separate
+linear path and no up-front nonlinearity probe: a linear model converges after
+its first increment and the loop self-terminates.
+
+### The velocity sub-KSP must converge, not just apply
+
+The velocity block is preconditioned by multigrid, but the KSP *around* that
+multigrid is FGMRES to `0.1 × tolerance` — never `preonly`. This is not a
+tuning preference. `PCFieldSplit` forms the Schur complement
+`S = A₁₁ − A₁₀ A₀₀⁻¹ A₀₁` and applies `A₀₀⁻¹` through this same velocity KSP,
+so `preonly` replaces `A₀₀⁻¹` with a single multigrid cycle and hands the
+pressure Krylov a *different* system `S̃ ≠ S`, preconditioned by a 1/μ mass
+matrix built for `S`.
+
+Measured on an annulus with a weak plane reaching the constrained boundary
+(`η₁/η₀ = 1e-3`, transversely isotropic): under `preonly` the pressure residual
+falls by 4.4e4 in about 16 iterations and then **stagnates at a floor ≈ 3.1e-7**,
+spending the remaining 184 iterations of its cap for nothing — on every outer
+iteration, 9 outer iterations in total. Under FGMRES it converges by 1.1e8
+monotonically in 17 pressure iterations, and the outer solve takes 1. The
+isotropic control moves the same way (5 → 1), so this is a property of the Schur
+application rather than of the anisotropy.
+
+Note the Schur application is bitwise reproducible under both settings, so the
+floor is *not* "the operator changes between applications" — `S̃` is a fixed
+linear operator, just the wrong one. The precise origin of the floor (most
+plausibly a range/null-space inconsistency between `S̃` and the constant-pressure
+null space attached to `S`) has not been isolated.
+
+The native (non-rotated) Stokes path and the GAMG fallback have always wrapped
+their multigrid this way; the rotated custom-FMG branch was the sole exception.
+`max_it` matches both (200); `rtol` matches the GAMG fallback (0.1 × tol), where
+the native path asks for 0.033 × tol.
+
+Two things make this easy to miss, and both are now instrumented. The outer
+iteration count does not reveal it — a full Schur factorisation with a good
+pressure mass still reports ~1 outer iteration while the inner solve grinds
+underneath. And a velocity FGMRES that exhausts its cap returns
+`KSP_DIVERGED_ITS`, which `KSPCheckSolve` deliberately does not escalate, so it
+would degrade silently where `preonly` simply could not fail; the rotated path
+now warns on it.
+
+`solver._rotated_freeslip_info` therefore reports `velocity_pc`, `schur_pre` and
+`velocity_pc_type` (PETSc's own view of the sub-PC), plus `vel_its_last` and
+`pres_its_last`. The last two are **last-application samples, not work**:
+`KSPGetIterationNumber` reports only the most recent solve, and the velocity KSP
+is applied once per Schur `MatMult`. They are named `_last` so they are not
+confused with the summed counts `solver_health` uses as its work axis; wiring
+the rotated path into `SolverInstrumentation.sub_reports()` is the proper fix and
+has not been done. On the opt-in `solver._rotated_use_lu` path there are no
+sub-KSPs at all and both fields record `None`.
+
+The manual loop exists because the rotated operator `Q A Qᵀ` carries no DM
+field information, so PETSc's DM-coupled fieldsplit cannot precondition it;
+the increment is solved by an IS-built fieldsplit instead, driven from the
+loop. The loop honours `zero_init_guess`, `picard`, and the solver's
+`consistent_jacobian` tangent policy (frozen / Newton / continuation).
+
+### How the datum passes through Newton
+
+The constraint is affine, and the design keeps every **accepted iterate
+feasible** (`u·n̂ = ũ_n` exactly, via an affine snap in the rotated frame), so
+each Newton increment satisfies the *homogeneous* constraint `n̂·δ = 0` — the
+datum never touches the tangent, the increment right-hand side, the FMG
+prolongation, or the null-space handling.
+
+The one deliberate exception is a **cold start with a non-zero datum**: the
+first increment carries the datum jump through the affine lift
+(`zeroRowsColumns(rows, diag, x̂, b̂)`) at the rest-state tangent, and is
+accepted without a line search. Snapping the zero state onto the datum instead
+manufactures an extreme boundary-strip strain state whose shear-thinning
+tangent produces a first step orders of magnitude too large — no line-searched
+step descends (measured). A warm start is assumed smooth and takes the exact
+snap directly.
+
+### Convergence semantics
+
+The loop converges on `‖F̂‖ ≤ rtol·ref + atol` with
+`ref = max(‖F̂(u₀)‖, ‖F̂(0)‖)`: the **rest-state residual is the intrinsic
+forcing scale**, so a good warm start (whose own initial residual is small) is
+not punished with an ever-stricter absolute target. A step-norm exit (tiny
+Newton step) is *verified* against the same reference before it may report
+convergence — a stiff tangent also produces tiny steps far from the solution,
+and an unverified tiny step goes back through the line search (a slow crawl
+continues; a genuine stall ends the loop with an explicit warning). An
+unconverged exit always warns; the fields hold the last iterate.
+
+### Null-space handling
+
+Rigid-rotation candidates (one mode on a closed circle, three on a spherical
+shell) are admitted to the increment null space and the post-solve gauge
+removal only if they pass **two** tests: tangential to every rotated
+constraint row (`Q·m ≈ 0` there), **and** a null vector of the assembled
+operator (`‖J·m‖ ≈ 0` against the velocity diagonal scale). The second test is
+what catches pinning the first cannot see — an essential condition on another
+boundary leaves the mode tangential to the rotated wall while the eliminated
+operator is *not* null on it, and admitting it projects an irreducible
+component out of every increment (the residual then floors far above
+tolerance instead of converging). The null space is therefore built after the
+first Jacobian assembly.
+
+### Parallel notes
+
+- Every branch the loop takes is decided **collectively**. The datum-activity
+  flag rides a PETSc `Vec` norm of the lift vector (the datum bookkeeping
+  itself is rank-local — only owners of datum boundary nodes hold entries —
+  and branching on it per-rank desynchronises the collective sequences: the
+  np>1 deadlock class).
+- All row surgery on vectors uses ownership-relative indices
+  (`_zero_rows_local` / `_set_rows_local`); indexing a local slice with global
+  rows is the np>1 crash class.
+- Direct LU per increment (`solver._rotated_use_lu = True`) is a **serial**
+  preconditioner-free diagnostic (the pressure-gauge pin is a naive per-rank
+  scan, marked `TODO(BUG)`); use it to separate "the operator/constraint is
+  wrong" from "the preconditioner is struggling".
+
+## Result and reaction
+
+The solve fills the solver's fields and stores a result dict
+(`solver._rotated_freeslip_info`) with the rotation, the constrained rows, the
+per-increment KSP iteration counts, the convergence verdict, and the
+**reaction** — the converged Cartesian residual `F(u)`, which for a linear
+residual equals `A·u − b` exactly. `boundary_normal_traction` projects the
+nodal reaction onto the boundary normal (corner-correct) and de-smears it with
+the shared boundary-mass machinery in `utilities/boundary_flux.py`;
+`dynamic_topography_field` writes `h = −(σ_nn − σ̄_nn)/(Δρ g)` onto a surface
+field for the free-surface integrator.
+
+## Tests
+
+`tests/test_1018_rotated_freeslip.py` (serial: essential-equivalence, FMG,
+tangent policies, datum linear + nonlinear),
+`tests/parallel/test_1066_rotated_datum_parallel.py` (np≥2: partition
+independence of the linear datum and the nonlinear Newton datum path).