The discrete adjoint: built into the solvers, driven from the run transcript

CLAUDE.md

@@ -344,9 +344,12 @@ to impose `v·n̂ = 0`:
   a consistency error against the faceted assembly. See
   `docs/developer/subsystems/rotated-freeslip.md` ("Which normal to use").
 - Works **inside the nonlinear SNES** and with **geometric FMG**. It honours
-  `solver.consistent_jacobian`: use `True` (consistent Newton) for smooth nonlinear
-  rheologies; `"continuation"` (staged Picard→Newton) for robustness far from the
-  solution. The rotated constraint is transparent to the tangent.
+  `solver.consistent_jacobian`: `True` (consistent Newton, the default — the
+  residual is symbolic, so the tangent is exact); `"continuation"` (staged
+  Picard→Newton) for robustness far from the solution; `False` (Picard) only where
+  a hard-yield viscoplastic solve needs it as an entry requirement. The rotated
+  constraint is transparent to the tangent, and the rotated path takes no warm-up
+  sweep before Newton.
 - The constraint **reaction** is the boundary normal traction σ_nn
   (`solver.boundary_normal_traction(boundary)` / `solver.dynamic_topography(...)`) — no
   augmented-Lagrangian splitting.

docs/developer/design/run-plan-and-transcript.md

@@ -146,6 +146,61 @@ insists that it is.** The signature requires a `dt`, so there is no container
 for "the next task". If the event clock is the general thing, the timestep is
 the common case rather than the definition.
 
+## Where the adjoint lives, and where it stops
+
+The transcript now supplies two of the three things a discrete adjoint needs:
+the ordered operator list, and the state each operator was linearised about
+(a snapshot before the operator, bit-exact on restore). The third — the
+linearisation itself — is a contract on each operator, not a pass over the
+record: an operator provides it or declines with a reason.
+
+The declining is recorded first. Every `solve`, `history_shift` and
+`swarm_advect` event carries `adjoint: {supported, reason}`, written when it
+ran. The verdicts are structural: an implicit step is a residual (Jacobian
+transpose for the state, symbolic derivative for a parameter); a rotated
+constraint solves inside its own Krylov loop with no transpose path; an
+unconverged solve is linearised about a state it never reached; a
+semi-Lagrangian trace is differentiable in the velocity but its interpolation
+at the departure points is not materialised; a particle step is adjointable
+exactly when the particle set is fixed across it, which `swarm.advection`
+checks by counting.
+
+Read back, the verdicts partition the run (`transcript_adjoint_segments`).
+That partition is what data assimilation needs rather than perfect
+invertibility: strong-constraint adjoint within a segment where every operator
+is smooth, and across a refusal a control variable with an error covariance —
+weak-constraint 4D-Var, with the joins chosen by the run. The optimiser needs
+a descent direction that is the same inexact direction each iteration, not an
+exact gradient; the exact discrete adjoint is the verification anchor where
+the operators admit it, and the segments say where that anchor holds.
+
+Two things follow from the residual being symbolic. First, every first
+derivative is always available: ∂R/∂u and ∂R/∂m are differentiated, not
+approximated, so the gradient is never in question — and the tangent the
+forward *iteration* used is irrelevant to it. Picard iterations spoil
+nothing; the converged state is the same, and the adjoint assembles ∂R/∂u at
+that state itself. Second, the same is not automatically true at second
+order. A Hessian — for posterior covariance, or a Newton step on the outer
+optimisation — needs ∂²R/∂u², ∂²R/∂u∂m, and a yield law written with `Min`
+or a softmin has a second derivative that is a distribution at the yield
+surface. Those terms exist symbolically, but they have to be handled with
+care rather than differentiated and trusted.
+
+What follows from it, in order: `adjoint_solve`, `dual_of` and
+`sensitivity` on the solvers — landed, checked against finite differences
+on Poisson, on a non-symmetric SUPG step, and on Stokes with a linear and a
+strain-rate-dependent viscosity, with the consistent tangent assembled for
+the adjoint whichever tangent the forward iteration used; the reverse driver
+(`uw.adjoint.TranscriptAdjoint`) — landed: it walks the transcript backwards,
+restores each step's snapshot, replays each solve to its own input state, and
+reads what each solve depends on from its residual, checked to 1e-7 against
+finite differences on a two-solver run, including a field read through its
+gradient (the Crank–Nicolson old flux), assembled as a FEM load rather than
+by parts; the two transport operators materialised — interpolation at
+departure points and ∂X_dep/∂v, which lift the semi-Lagrangian refusal; and
+a Taylor test in the library (`test_0020`, and the sinker example through
+the library at 1.00000).
+
 ## Inferred plan, then declared plan
 
 The plan is **inferred** today — the figure takes the most common step as the

docs/developer/guides/HOW-TO-WRITE-UW3-SCRIPTS.md

@@ -642,6 +642,129 @@ shown, and your own filters still apply. "The velocity block fell back to
 gamg" changes what the numbers mean, and a record that kept the residual norms
 but not that line would be an account of the run with the explanation removed.
 
+**Every operator says whether it admits a discrete adjoint.** Each `solve`,
+`history_shift` and `swarm_advect` event carries a verdict, written when the
+operator ran:
+
+```json
+{"kind": "solve", "name": "SNES_Stokes(v)",
+ "adjoint": {"supported": true,
+             "reason": "implicit residual: Jacobian transpose for the state, symbolic derivative of the residual for a parameter"}}
+```
+
+The verdict is structural — about the operator as configured, not about
+whether a driver exists yet — so a run says where its adjoint breaks *while it
+runs*. What refuses, and why:
+
+- a rotated constraint (free-slip or fault contact): the solve runs on a
+  rotated operator inside its own Krylov loop, with no transpose path;
+- a solve that did not converge: a linearisation about a state the solve
+  never reached is not the adjoint of anything — the outcome overrides the
+  structural verdict after the fact;
+- a semi-Lagrangian history: the departure-point trace is differentiable in
+  the velocity, but the interpolation at the departure points is not
+  materialised as an operator;
+- a swarm step whose particle set changed — `swarm.advection` records the
+  count before and after, and a particle removed on leaving the domain
+  changes the dimension of the state. The rule is one line: a particle step
+  is adjointable exactly when the particle set is fixed across it.
+
+An Eulerian or SUPG history is supported — an implicit step is a residual,
+and the SUPG adjoint that passed its Taylor test at 1.00000 is exactly that
+case. The text transcript notes where the adjoint breaks, once per change
+rather than on every step.
+
+`uw.transcript_adjoint_segments(source)` reads the verdicts back as the
+partition they imply — maximal runs of steps whose every operator admits an
+adjoint, separated by the steps where one refused. That partition is the
+assimilation window's structure: strong-constraint adjoint within a segment;
+across a refusal, a control variable and an error covariance, which is
+weak-constraint 4D-Var with the joins chosen by the run rather than by hand.
+Nothing is approximated silently — the refusal says what the model was
+allowed to be wrong about.
+
+**The adjoint of one solve is built in.** For a solver whose verdict is
+"supported", the discrete adjoint is two calls, with no hand algebra:
+
+```python
+b = -solver.dual_of(T.sym[0] - T_target.sym[0])   # -dJ/dT for J = 1/2 int (T - T*)^2
+mu, reason = solver.adjoint_solve(b, target=mu_var)  # K^T mu = b, K the SNES Jacobian
+dJ_dkappa = solver.sensitivity(mu_var, kappa)        # int (dF/dkappa) . mu, symbolic dF/dkappa
+```
+
+`dual_of` assembles the right-hand side on the solver's own space, so the
+Dirichlet nodes are excluded and the multiplier comes back zero there — the
+homogenised adjoint conditions, without stating them. `sensitivity` follows
+the parameter through the constitutive model's own symbol (the residual holds
+`\upkappa`, whose value is your `kappa`), so the chain rule reaches it.
+
+One thing to get right, because `solve()` moves it: a time step's residual is
+`F(u_new; u_old, v, dt)`, and the history manager shifts `u_old` out of its
+slot in the post-solve hook. Put the step's input back before linearising —
+`solver.DuDt.psi_star[0].array[...] = u_old` — or the sensitivity is a few
+per cent wrong on a SUPG step (measured).
+
+Stokes takes the same transpose on its composite (u, p) system, with
+`target=(u_adj, p_adj)` and `dual_of` taking a velocity-space expression.
+With a linear viscosity the operator is symmetric, and this reproduces the
+second-solver construction in `docs/examples/adjoint`. With a strain-rate- or
+pressure-dependent viscosity the adjoint is the transpose of the **consistent
+tangent** ∂R/∂u, which that construction cannot build. Picard iterations in
+the forward solve spoil nothing — the converged state is the same, and ∂R/∂u
+is a function of that state alone — but they leave the SNES holding the
+frozen-viscosity Jacobian *kernel*. So when the forward ran Picard on a
+nonlinear residual, `adjoint_solve` switches the kernel to the consistent
+tangent for its assembly (a JIT rebuild; the DM and KSP are kept), transposes
+that, and puts the Picard kernel back for the next forward solve. The
+verdict says so.
+
+**The whole run, backwards.** With `model.record_every = 1` the transcript is
+a forward tape — the operators per step, and the state each step started
+from — and `uw.adjoint.TranscriptAdjoint` walks it in reverse with no
+problem-specific wiring:
+
+```python
+final = model.save_state()                       # the N+1th level
+back = uw.adjoint.TranscriptAdjoint(model, final)
+result = back.gradient(misfit_integrand, parameters=[eta0], fields=[beta])
+result["parameters"][eta0]                       # dJ/d eta0
+result["fields"][beta]                           # dJ/d beta_0, as a dual field
+```
+
+For each solve, in reverse order of the record, it restores the step's
+snapshot, replays the solves before it, replays it, and puts each history's
+input back where the post-solve hook shifted it — so the residual is
+linearised at the solve's own input state without anyone touching
+`psi_star`. The residual then says what the solve read: every field in
+`F0`/`F1` other than the unknown gets the dual `(dR/df)^T mu`, a history
+slot's dual goes to the field it tracks at the previous level, and every
+parameter gets `mu^T dR/dm`. A field read through its *gradient* — a
+Crank–Nicolson step (θ = 0.5, the `AdvDiffusion` default) reads the old
+level as `κ∇T_old` — gets the gradient part of the load too: the dual is
+assembled as the FEM load `∫ g₀ φⱼ + g₁·∇φⱼ` by a generic solver's residual
+at zero, so there is no integration by parts and no boundary term to drop.
+A dual is held as a field (one coefficient per node); pair it with a
+direction using `uw.adjoint.inner(field, dual, direction)`, which sums over
+the owned degrees of freedom (a NumPy dot on `.array` counts a partition's
+ghost nodes twice), so a control `c` with `f_0 = f_0(c)` finishes with
+`inner(f, dual, d f_0 / d c)`.
+
+Two things the tape has to contain. Every solve must be inside a step — a
+Stokes solve taken before the loop to make `v_0` is invisible to the walk,
+and its dependence on the parameters with it. And a driver that runs the
+forward model more than once must reset the Eulerian history each time it
+sets the initial condition (`adv.DuDt.initialise_history()`), or the second
+run reads the first run's history. Checked in `tests/test_0020` on a
+two-solver, two-step sinking blob: the viscosity gradient and the dual on
+the initial level set both match central finite differences to 1e-4
+(measured 1e-7), at θ = 1 and at the default θ = 0.5, serially and on two
+ranks.
+
+All of it is checked against central finite differences in `tests/test_0019`:
+Poisson; one SUPG step, where the Jacobian is not symmetric and a transpose
+taken the wrong way round would show; Stokes with a constant viscosity; and
+Stokes with η(ε̇) under the consistent tangent.
+
 The figure marks the same three states per solve — converged, converged with a
 fieldsplit block that hit its iteration cap, and diverged. The middle one is
 worth the separate mark: a capped block did not solve, so the Schur operator

docs/developer/guides/adversarial-review.md

@@ -82,6 +82,12 @@ silently understates the run. See
 [HOW-TO-WRITE-UW3-SCRIPTS](HOW-TO-WRITE-UW3-SCRIPTS.md) and
 `docs/developer/design/run-plan-and-transcript.md`.
 
+**Every operator gives an adjoint verdict.** A `solve`, `history_shift` or
+`swarm_advect` event carries `adjoint: {supported, reason}`, written when it
+ran. A new history scheme without `_adjoint_support()` fails
+`tests/test_0018_adjoint_support_record.py`; a new operation on model state
+that records no verdict lets a run claim invertibility it does not have.
+
 **Named quantities keep their names.** A coefficient written as
 `uw.expression(r"\rho_0 \alpha g", ...)` appears in the description under that
 name. An anonymous float collapses into the assembled product and the