Split-level fault multigrid, the thin-sheet solver, and sheet extrusion

CLAUDE.md

@@ -40,6 +40,24 @@ fenced block for display, LaTeX not unicode). Private cross-project context,
 when a session has it, still arrives through `UW_AI_TOOLS_PATH` (colon-separated
 directories of `.md` files; set via `./uw setup`).
 
+**Pending work, before anything else.** Repeated work has come from investigating
+without looking at what is already in flight. Before an investigation, a new branch
+or a design decision, and again when a related PR merges mid-session:
+
+```bash
+git fetch origin development
+git log --oneline --merges -20 origin/development            # what merged recently
+```
+
+plus the open PRs (GitHub MCP `list_pull_requests`, or `gh pr list --state open`). Read
+every open or recently merged PR that touches the subsystem before writing a line, name
+them in the opening report, and base the work on current development (branch from it,
+or merge it into the base branch first). A related PR open for more than about a week
+is itself a finding: say so at the start and ask whether to merge it, build on it, or
+wait. Do not rely on the maintainer noticing the overlap. (Case: the fault multigrid
+session of 2026-10-13 re-diagnosed parallel fault defects that PR #672, open for
+weeks, had already fixed, then paid a seven-file merge.)
+
 **Inline TODOs** mark the *place*; the roadmap item or issue tracks the *work*:
 
 ```python

docs/developer/design/FAULT_MULTIGRID_HIERARCHIES.md

@@ -0,0 +1,185 @@
+# Multigrid hierarchies for split-node faults
+
+Technical note, 2026-10-13. The exploration that chose the hierarchy; the library
+implementation (C, below) is governed by
+[subsystems/fault-split-multigrid.md](../subsystems/fault-split-multigrid.md).
+The exploration drivers are kept with the campaign, outside this repository
+(`stage*.py` monkey-patch the solver; `library_c*.py` and `thick_parallel.py` go
+through the library).
+
+## The question
+
+A split-node fault lives on the finest mesh. The geometric-multigrid levels below it
+come from the add-on hierarchy: the static background mesh, which is never cut, and
+the add-on levels refined from it toward the fault. How the fault is represented on
+those coarser levels decides whether the velocity solve stays mesh-independent as the
+fault is refined.
+
+## Alternatives
+
+| | finest | add-on levels below it | background levels |
+|---|---|---|---|
+| **A** | split | welded, Galerkin | welded, Galerkin |
+| **B2** | split | welded, TI weak band, rediscretised | Galerkin from the TI levels |
+| **C** | split | split, Galerkin | welded, Galerkin |
+| **D** | split | split, Galerkin | broad TI weak zone, rediscretised |
+
+A is the current add-on behaviour (`mesh.adapt` then `add_fault`).
+
+TI band: a Gaussian in the distance to the trace, half-width equal to the level's
+near-fault cell size, director along the fault normal, $\eta_1 = \eta^{1-f}\eta_w^{f}$
+with $\eta_w = 0.1\,\eta$ for a frictionless fault and $\eta_w = \eta_f w$ (the
+compliance of the interface) for a viscous one.
+
+Split levels (C, D) need two things the current hierarchy does not provide.
+
+**A side-aware transfer.** For each fine node, one fine cell containing it is taken, and
+that cell's centroid is located in the coarse mesh. A centroid is never on the fault, so
+the coarse cell is on the node's own side. The coarse cell's Lagrange interpolant is
+then evaluated at the node (dual-basis identity $W = BM^{-1}$, as in
+`custom_mg.nested_refine_pair_prolongation`), so a node on the slit interpolates from
+its own side. The node-cloud builders cannot do this: they triangulate the coarse
+nodes, and coincident slit nodes make the transfer degenerate. The levels are not
+nested near the trace, because each cut pulls vertices onto it, so a fine node can lie
+slightly outside its located coarse cell. The evaluation point is therefore clamped
+into the cell (its barycentric coordinates clipped and renormalised). Without the
+clamp, extrapolated weights reached 5.4 near a tip, two coarse unknowns were reached
+only through them, and the Galerkin operator had a near-zero diagonal that made SOR
+return NaN. With it, every weight is at most 1.
+
+**A per-level pair frame.** Each split level carries its own pair rotation (mean and
+jump), and its opening (normal-jump) unknowns are removed:
+$\hat P_l = Q_{l+1} P_l Q_l^{T}$, restricted to admissible unknowns. Coarse unknowns
+that no fine unknown reaches at all (a few near the tips) are dropped too; they are null
+modes of the Galerkin coarse operator.
+
+## Test problem
+
+- Box $[0,1]^2$ with simple shear on the walls; incompressible Stokes, $\eta = 1$.
+- A straight fault, $(0.30, 0.45)$ to $(0.70, 0.55)$, frictionless or viscous ($\eta_f = 1$).
+- Background: `UnstructuredSimplexBox(cellSize=0.1, refinement=1)` (levels $h = 0.1, 0.05$).
+- Add-on level $j$: `base.adapt` to near-fault $h = 0.05/2^j$, built independently from the
+  background for each $j$, so the levels are conforming but not nested; split by
+  `add_fault` where the alternative asks.
+- $k$ add-on levels in total, the finest at $h = 0.05/2^k$.
+- Rotated fieldsplit solve: velocity block custom FMG, pressure Schur 1/μ mass.
+- The measure is velocity-block iterations per solve.
+
+## Results
+
+Frictionless fault:
+
+| $k$ | near-fault $h$ | A | B2 | C | D |
+|---|---|---|---|---|---|
+| 1 | 0.025 | 9 | 10 | 10 | 14 |
+| 2 | 0.0125 | 13 | 14 | 7 | 9 |
+| 3 | 0.00625 | 18 | 22 | 7 | 6 |
+| 4 | 0.0031 | 23 | 30 | 6 | 6 |
+| 5 | 0.0016 | 32 | 46 | 5 | 5 |
+
+Viscous fault:
+
+| $k$ | A | B2 | C | D |
+|---|---|---|---|---|
+| 1 | 9 | 9 | 9 | 17 |
+| 2 | 12 | 19 | 7 | 9 |
+| 3 | 17 | 53 | 6 | 8 |
+| 4 | 22 | 200 (cap) | 6 | 7 |
+| 5 | 30 | 200 (cap; stopped at 1 h) | 5 | 6 |
+
+The peak slip agrees to five figures across all alternatives at each $k$, so only the
+solver path changes.
+
+Coulomb fault ($\mu = 0.6$, $\sigma_n = 0.8$; four Newton steps in every case),
+velocity iterations per Newton step:
+
+| $k$ | A | C |
+|---|---|---|
+| 1 | 5, 15, 17, 17 | 5, 15, 17, 17 |
+| 3 | 9, 21, 24, 24 | 4, 14, 15, 16 |
+| 5 | 15, 38, 45, 50 | 3, 13, 15, 15 |
+
+Supporting measurements:
+
+- **A single TI level.** One TI level directly below the split finest level, with
+  Galerkin below it, is worse than A at every band strength. At $k = 4$ and a band
+  half-width of half a cell, A takes 27; the band at 0.1, 0.03 and 0.01 takes 29, 46
+  and 73. With a half-width of a full cell the band stalls.
+- **Transfer alone.** C's levels left welded but with the side-aware transfer behave
+  like A (35 iterations at $k = 5$). The gain comes from splitting the levels.
+- **3-D extruded layer (preliminary).** Only the through-thickness levels split, the
+  lateral levels welded: 11 to 16 iterations as $h$ halves, against 14 to 25 for A.
+
+## A moving fault
+
+The same box and fault, translated in 8 steps from $x = 0.3$ to $x = 0.7$. The
+background is built once; each step re-adapts the add-on part from it with one
+`mesh.adapt` call, splits the finest level and, for C and D, splits each intermediate
+level of that same call. The table gives the mean seconds per step. The runs shared a
+busy machine, so absolute times are good to a factor of 1.5 to 2; the ratios between
+strategies and the iteration counts are not affected.
+
+$k = 5$ (seven levels, near-fault $h = 0.0016$):
+
+| | adapt | split finest | split levels | level frames | background TI | solve | total | velocity its |
+|---|---|---|---|---|---|---|---|---|
+| A | 3.1 | 1.9 | – | – | – | 137 | 143 | 30–34 |
+| C | 2.7 | 1.7 | 2.8 | 1.0 | – | 38 | 46 | 5–6 |
+| D | 4.7 | 2.9 | 5.1 | 1.8 | 2.4 | 61 | 79 | 4–5 |
+| GAMG | 3.8 | 2.4 | – | – | – | 127 | 134 | 200 (cap) |
+
+At $k = 3$ the totals are 27 (A), 15 (C), 17 (D) and 28 (GAMG) seconds per step. The
+solve includes building the transfers: 2 s for A, 3.4 s for C. With the clamped
+transfer, C takes 5 to 6 velocity iterations at every position at $k = 3$ and 4 to 8 at
+$k = 5$.
+
+The peak slip is identical across A, C and D at every position. Level frames come from
+throwaway solvers in the drivers (`build_rotation`); they could be computed from the
+pair geometry directly.
+
+## Parallel
+
+A, the existing path, had three defects at np > 1:
+
+1. **Fixed** (`test_1086`). The rotated solve zeroed the opening rows of its fine
+   transfer by local position instead of global row.
+2. **Open, now harmless.** `mesh.adapt` keeps a different intermediate level for the
+   multigrid tail in parallel than in serial (1358 against 1520 cells; the finest mesh
+   is the same). The levels need not match: the transfer locates nodes.
+3. **Fixed** (`test_1089`). `add_fault` gathers the fault's cell star onto one rank, so
+   the split finest level is not co-partitioned with its tail. The node-cloud transfer
+   left 618 coarse unknowns with no fine image, their nearest-node injection
+   duplicated 512 columns, the Galerkin coarse operators became singular, and PETSc's
+   inode SOR returned NaN. Every pair touching a split level now takes the
+   cell-located transfer, which asks the other ranks for the nodes it cannot locate
+   locally, and unreached unknowns are dropped rather than repaired: A converges at
+   np = 2 (15 iterations against 13 serial), C in 6 at np = 2 and 7 at np = 4 (k = 2).
+
+The pairing rides every redistribution as a pair of DM labels. The migration star
+forest PETSc returns for an initial distribution does not index the stratified point
+numbering (three quarters of the points mapped to the wrong depth), so it is not used.
+
+## What works best
+
+**C: split every add-on level and leave the background welded.** The iteration count
+falls slightly with refinement, for frictionless, viscous and Coulomb faults alike.
+It needs nothing rediscretised and no TI parameters. With a moving fault its extra
+surgery and transfers cost about 5 s per step against a solve that is three times
+cheaper than A's at $k = 5$.
+
+D matches C at depth but is worse on shallow hierarchies, and it adds a rediscretised
+background operator to re-assemble at every step and band parameters to choose for
+every problem, for no gain. B2 is worse than A, and badly so for the
+viscous fault. A smeared band on a coarse level is a continuous shear that the split
+fine operator rejects, so the coarse space has to contain the discontinuity itself.
+
+## State (library, 2026-10-13)
+
+In the library: the add-on levels of one `adapt` call, each split by `add_fault`; the
+cell-located transfer (serial and parallel) and the per-level frames in
+`utilities/fault_mg.py`, installed through `custom_mg` and the rotated solve; the 3-D
+layer with every level split, extruded in parallel; 3-D sheet faults on add-on meshes.
+Not yet: fault networks through the levels (a network is replayed on every level, but
+a one-cell junction ligament closes at the coarser spacing and the replay then stops
+there, leaving that level and those beneath welded), faults near a boundary,
+co-locating the levels with the finest as an optimisation.

docs/developer/design/THIN_SHEET_VISCOELASTIC_DESIGN.md

@@ -0,0 +1,140 @@
+# Thin viscoelastic sheet with faults — flat and on the sphere
+
+**Status**: In progress (started 2026-10-08). Branch `feature/thin-sheet`.
+
+## Goal
+
+A vertically averaged (thin-sheet) model of the lithosphere carrying faults, on a
+flat 2-D mesh and on the surface of a sphere, with a Maxwell viscoelastic
+rheology and a viscous substrate. Headline result: the loading, slip and
+postseismic relaxation of a strike-slip fault in a viscoelastic sheet, flat first
+and then on the sphere.
+
+## The sheet
+
+With thickness $L$, tangential velocity $\mathbf u$ and the plane-stress closure
+$\sigma_{nn} = 0$ on the dynamic stress, the pressure of an incompressible
+material is eliminated: $p = \tau_{nn} = -\operatorname{tr}_s\boldsymbol\tau$, so
+
+$$
+\boldsymbol\sigma = \boldsymbol\tau + \operatorname{tr}(\boldsymbol\tau)\,\mathbf P,
+\qquad \mathbf P = \mathbf I - \hat{\mathbf n}\hat{\mathbf n}^T,
+\qquad
+-\nabla_s\cdot(L\boldsymbol\sigma) + c\,\mathbf u = \mathbf f .
+$$
+
+$\boldsymbol\tau$ is any deviatoric constitutive law of the in-plane strain rate
+$\dot{\boldsymbol\varepsilon}_s = \mathbf P\,\mathrm{sym}(\nabla_s\mathbf u)\,\mathbf P$
+— Newtonian, Maxwell (the stress history is the in-plane deviatoric stress;
+$\tau_{nn}$ follows from the trace), or yielding. $c$ is a basal drag
+($\eta_a/h_a$ for a viscous substrate); $\mathbf f$ a force per unit area.
+
+Consequences used as hard baselines:
+
+- uniaxial stretch with free lateral edges: $\sigma_{xx} = 3\eta\dot\varepsilon$ (Trouton);
+- rigid motion: $\boldsymbol\sigma = 0$ exactly;
+- Maxwell uniaxial loading: $\sigma = 3\eta\dot\varepsilon\,(1 - e^{-t\mu/\eta})$;
+- a fault in a thin sheet: Barr & Houseman (1996) Appendix (A10)–(A13) — the
+  plane-strain (A1)–(A9) solution is already in `underworld3.analytic`; the
+  thin-sheet one is to be added beside it;
+- stress diffusion with basal drag after imposed slip (Elsasser; Lehner, Li &
+  Rice 1981).
+
+Not in the first round (recorded so they are not forgotten): thickness
+evolution $\dot L = -L\,\nabla_s\cdot\mathbf u$ with the gravitational potential
+energy gradient (England & McKenzie's Argand number); yield in the sheet, whose
+invariant must include $\tau_{nn} = -\operatorname{tr}\boldsymbol\tau$ (the
+plane-strain code omits it).
+
+## Implementation
+
+`uw.systems.ThinSheet` (`systems/thin_sheet.py`) is an `SNES_Vector`: there is
+no pressure unknown, so no saddle point. Its inner `_Unknowns` hands the
+constitutive model the projected strain rate through `Unknowns.E` (the VEP model
+already reads it; `ViscousFlowModel.grad_u` now does too — identical on a volume
+mesh). Templates: `F1 = L (τ + tr τ P)`, `F0 = -f + c u`.
+
+### On the sphere
+
+The velocity is carried in the embedding space (3 components on a
+`SphericalManifold`, `dim=2`, `cdim=3`). PETSc already supplies the surface
+gradient along the flat facets in the derivative index; the outer projection in
+$\dot{\boldsymbol\varepsilon}_s$ supplies the curvature term
+$-\mathbf u\cdot\nabla_s\hat{\mathbf n}$. Two normals are involved and they
+are deliberately different:
+
+- the **projector** uses the facet's own normal (a cell-wise constant field
+  built once per mesh). The gradient is taken along the flat facet, so a rigid
+  rotation — a Killing field, strain-free on the sphere — is strain-free in the
+  discrete form only with the facet normal; with the analytic normal it keeps a
+  spurious strain of order $h$ (measured 0.7% of the velocity at $h = 0.25$);
+- the **tangency constraint** $\mathbf u\cdot\hat{\mathbf n}=0$ uses the
+  analytic normal of the surface, imposed strongly at every node through the
+  per-node rotation of `utilities/rotated_bc.py` (`SHEET_TANGENCY`, the same
+  machinery as rotated free-slip and the split-node fault constraint, so faults
+  compose with it). The nodes lie on the sphere, so a rigid rotation is exactly
+  admissible and is reproduced to solver tolerance.
+
+A penalty was used only to probe feasibility and is not an API. The toroidal
+manufactured solution converges at second order (the flat-facet geometry limit
+for P2) with tangency at machine precision.
+
+Vector solves on a manifold needed four `dim → cdim` corrections (the manifold
+vector-FE gap): the `SNES_Vector` FE component count, its PC block size and
+default unknown, `maths/vector_calculus` (component count), and the
+`sympy.vector` truncation in `_jitextension`.
+
+### Faults on the sphere (Phase 5)
+
+The split-node pipeline is topological and survives; what has to generalise is
+geometry: great-circle crossings in `line_cut`, cell orientation by the outward
+normal in `reconnect`, the contact frame $\{\hat\nu = \hat r\times\hat t,\ \hat t,\
+\hat r\}$ in `fault_contact` (PETSc's FVM edge normal is wrong for `cdim=3`), and
+the `mesh.dim` strides in `rotated_bc`, whose rigid-rotation gauge must offer the
+three Killing fields. `mesh.adapt` (NVB) runs on the manifold already; refined
+vertices need snapping back to the sphere.
+
+## Phases
+
+0. Feasibility probes (sheet operator on S² with penalty and with the strong
+   rotation; NVB adapt on S²; splitting an arc) — scratch, not committed.
+1. Flat viscous `ThinSheet` + tests (`tests/test_1075_thin_sheet.py`).
+2. Viscoelastic sheet (Maxwell; small-strain non-advected history for fault
+   problems; Elsasser check).
+3. Faults in the flat sheet (rotated/fault BCs hoisted off the saddle class; a
+   vector-only branch of the rotated solve; BH96 thin-sheet solution).
+4. Sheet on the sphere with strong tangency and the Killing null space.
+5. Faults on the sphere.
+6. Results and examples.
+
+## One model, thin or thick
+
+Principle (maintainer, 2026-10-10): the same physics is the same model script. A model
+must move between the thin sheet (`ThinSheet` on a surface mesh) and a thick sheet
+(`Stokes` on a thin 3-D mesh) by changing the mesh and the solver class, nothing else.
+The solver owns the kinematics, the closure and the geometry; the model owns the
+rheology. What is already shared: the constitutive models, the fault interface API, the
+boundary-condition methods, the stress history and `solve(timestep=)`.
+
+State of the alignment (2026-10-12):
+
+| | thin sheet | thick sheet | aligned form | state |
+|---|---|---|---|---|
+| construction | `ThinSheet(mesh, v, p, degree, p_continuous)` | `Stokes(mesh, v, p, ...)` | same signature; the sheet's `p` is the diagnostic $-\operatorname{tr}\tau$, refreshed after each solve | done |
+| mesh | the sheet mesh | `mesh.extrude(L, layers)` | the thick twin keeps the sheet's boundaries as walls (same names) and adds `Base` / `Surface` | done (serial); MG tail coarsens through the thickness first |
+| forcing | body force per unit volume, times $L$ | per unit volume | `bodyforce` means the same on both | done |
+| basal coupling | `add_basal_drag(c, "Base", reference_velocity=)`, `add_basal_traction(t, "Base")` | the same, as natural conditions on `Base` | one line serves both; the sheet ignores the boundary name | done |
+| fault geometry | `mesh.add_fault` | extruded with the mesh: a vertical fault through the layer, split at surface and base, tips shared | `add_fault_bc` on either | done flat (test_1083: velocity 0.5 %, slip 1.6 %); sphere: first order in L, 1.35 % at L = 0.02; the multigrid tail uses fault-free coarse levels |
+| stress | `solver.stress` (in-plane, tension positive) | `solver.stress` (3-D) | `uw.maths.tensor.resolved_traction` / `coulomb_failure_stress` on `solver.stress` | done |
+| nonlinear invariants | `Unknowns.second_invariant` adds $M_{nn} = -\operatorname{tr}M$ | plain contraction | every model invariant goes through the solver | done |
+| Newton tangent | `_flux_closure` | `_flux_closure` | the model's `flux_jacobian` is composed with the solver's closure (it used to replace the whole flux, and was ignored on vector solvers) | done |
+| topography | implied isostatic rate | free surface | with thickness transport (round 2) | open |
+
+Mesh-pairing notes. A thin 3-D fault cannot come from the 3-D splitter: it keeps
+the patch rim welded (the tip rule one dimension up), so a fault through a layer
+would be locked at the surface and the base, and it needs about nine cells of
+depth. Extruding the already-split sheet avoids both: the slit is duplicated at
+every level, so the fault is open where it meets the faces, and the two sides are
+cut into tetrahedra identically (the minus copy of a vertex takes its plus copy's
+place in the vertex order), so the pairing is rebuilt by joining mapped vertex
+tuples.