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

[lmoresi]

What this is

The discrete adjoint, built into the solvers and driven from the run transcript. Stacked on #716, which supplies the step boundary, the snapshots and the record it walks; the diff here is exactly the adjoint additions.

An implicit step is a residual R(u; m) = 0 and the residual is SymPy, so every first derivative is available symbolically: ∂R/∂u is the Jacobian the SNES already assembles, ∂R/∂m and the coupling to every field the residual reads are differentiated, not approximated. The adjoint of one solve is a transpose against that Jacobian; the adjoint of a run is the chain of them in reverse order of the record.

In every solver

b  = -solver.dual_of(T.sym[0] - T_target.sym[0])   # -dJ/dT for J = ½∫(T - T*)²
mu, reason = solver.adjoint_solve(b, target=mu_var) # Kᵀ mu = b
dJ_dkappa  = solver.sensitivity(mu_var, kappa)      # ∫ (∂F/∂κ)·mu, symbolic ∂F/∂κ

• adjoint_solve transposes against the assembled Jacobian via KSP.solveTranspose. Dirichlet conditions come out homogenised for free: the global vector holds only unconstrained dofs. Stokes does it on the composite (u, p) system.
• dual_of assembles ∫ e φⱼ on the solver's own constrained space as a generic solver's residual at zero — no linear solve.
• sensitivity differentiates the residual through the constitutive model's own symbol: the residual holds \upkappa whose value is the user's kappa, and uw.function.derivative substitutes constants' values before differentiating, so both naive routes give zero. _peel_except expands every named expression except the one being differentiated against.
• If the forward solve used the Picard tangent, the adjoint assembles the consistent tangent itself (a kernel rebuild; the DM and KSP are kept) and puts the Picard kernel back for the next forward. Picard iterations do not spoil the adjoint — the converged state is the same and ∂R/∂u is a function of the state alone; only the compiled kernel was the wrong matrix. Under "continuation" the adjoint assembles at α = 1.

Over a whole run

back   = uw.adjoint.TranscriptAdjoint(model, model.save_state())
result = back.gradient(misfit, parameters=[eta0], fields=[beta])
result["parameters"][eta0]                     # dJ/dη₀
uw.adjoint.inner(beta, result["fields"][beta], dbeta0_dc)   # dJ/dc through β₀

No problem-specific wiring. For each solve, in reverse order of the record: restore the step's snapshot, replay the solves before it, replay it, put each history's input back in the slot the post-solve hook shifted; transpose-solve with the accumulated dual on its unknown; then read the residual to see what the solve depended on. Every other field in F0/F1 gets (∂R/∂f)ᵀμ as a dual on its own space — through its value and, for a field read via its gradient (a Crank–Nicolson step reads κ∇T_old), through the gradient part of the load, assembled as ∫ g₀φⱼ + g₁·∇φⱼ with no integration by parts. A history slot's dual goes to the field it tracks at the previous level; every parameter gets μᵀ∂R/∂m, plus the explicit ∂J/∂m.

Every operator in the transcript carries a verdict, written when it ran: adjoint: {supported, reason}. A rotated constraint (no transpose path through its own Krylov loop), an unconverged solve, a semi-Lagrangian history (interpolation at departure points not materialised), a swarm step whose particle set changed — each refuses with the reason, and transcript_adjoint_segments reads the verdicts back as the strong/weak-constraint partition of the assimilation window.

The default tangent is Newton

consistent_jacobian = True is now the default; Picard is the opt-in, by name, where a hard-yield solve needs it as an entry requirement. The residual is symbolic so the tangent is exact and cheap; Picard converged linearly and left the SNES holding a Jacobian that was not ∂R/∂u for everything downstream. The design note had already named this the intended default and marked it "benchmark before flipping"; the full suite was that benchmark. What it turned up: the sharp-yield cold solve that Picard could not do converges under Newton (test_1057 now sets Picard explicitly on its cold half); and the cold-start "warm-up" on the saddle-point path degraded a linear ksponly solve under Eisenstat–Walker (12% vs 2% on the spherical-shell Nitsche response) — now skipped for ksponly, and honestly described: it is an nrichardson residual sweep, not a Picard step, and picard=-1 switches it off.

Verified

Every case against central finite differences:

case	agreement
Poisson, dJ/dκ	1e-4
one SUPG step (non-symmetric K, where a wrong transpose shows)	1e-3
Stokes, constant viscosity	1e-4
Stokes, η(ε̇), consistent tangent / continuation / Picard forward	1e-3
two-solver, two-step run: dJ/dη₀ and the dual on β₀ (θ = 1 and θ = 0.5)	1e-7 measured, serial and np = 2
the sinking-blob example, blob centre as control, through the library	Taylor ratio 1.00000 / 0.99997

The example (docs/examples/adjoint/sinker_transcript) keeps its hand-rolled adjoint beside the one-call version, and both Taylor tests give the same gradient to six digits.

Adversarial review

Three reviewers ran probes against the branch (posted on #716, fixed in fba44975 and later): two MPI deadlocks (a rank-local gate before a collective solve; rank-local size checks), a ghost double-count in the dual pairing, the missing explicit ∂J/∂m, sixteen scratch variables leaked per call, a second adjoint_solve refusing after the Picard rebuild, and the warm-up's misdescription. All fixed here except dropping the warm-up, which is a benchmark item in the planning file.

Not in this PR

Semi-Lagrangian schemes still refuse (their two transport operators are not materialised); second-order terms are noted in the design note and not built.

🤖 Generated with Claude Code

https://claude.ai/code/session_01Na7qBenCp67rDTZhFGTh5V

[lmoresi]

Adversarial review — adjoint machinery and the Newton default (2026-09-14)

Two independent reviewers running probes. Findings that survived, and what was done. The fixes were made on the build branch before the split; that history is on archive/timestepping-pattern-history (fixes in fba44975, further in dd928797), and this PR's single commit carries them.

Adjoint machinery (adjoint_solve / dual_of / sensitivity, uw.adjoint.TranscriptAdjoint)

Sev	Finding	Evidence	Status
HIGH	_nonzero was a rank-local max gating a collective solve; a misfit supported on one rank's cells deadlocked gradient()	np=2, misfit ½|v|²·[y>0.85]: rank 0 solves, rank 1 skips; alive at 100 s vs 10 s serial	fixed: reduced over ranks, empty-array safe
HIGH	the field-gradient test's dual @ direction was rank-local and counted ghosts	np=2: −0.00602 / −0.00614 vs fd −0.010876; owned-dof dot = fd to 7.7e-10	fixed: uw.adjoint.inner() over owned dofs; test passes at np=2
MED	gradient() omitted the explicit ∂J/∂m	misfit ½η₀|v|²: −0.004824 vs fd −0.002404	fixed, with a test
MED	16 MeshVariables leaked per gradient() call; call time 1.3 → 14.5 s by the sixth	vars 5→101 over 6 calls, RSS 335→656 MB	fixed: scratch pool; test asserts no growth
MED	a second adjoint_solve after the Picard-path rebuild raised "no forward solve"	is_setup=False after _restore_tangent	fixed: the consistent kernel stays installed until the next forward build
LOW	dual-space scratch name collided when the BC set changed without changing its count	_dual_SNES_Poisson_29_1 twice	fixed: hashed signature
LOW	size-mismatch raises were rank-local, before collectives	by inspection	fixed: decided with an allreduce
LOW	AdvDiffusion at θ=0.5 refused ("reads a derivative of psi_star_…") — names the slot, and test_0020 uses θ=1 without saying why	_reads on the default solver	refusal now names the tracked field and says θ=1 avoids it; the gradient-read term (integration by parts with its boundary part) is an open item

Attacks that failed: token prefix collisions (symbols print braced); a Projection with no BCs (rel 2.3e-9); Vector_Projection through mesh.vector.jacobian (2.4e-9); a parameter two levels deep (2.6e-8); serial gradient() bit-repeatable ×4; state after gradient() is level 0 and the next forward reproduces J to 1e-15; Picard→consistent rebuild cycles identical over 3 rounds.

Newton default (fa83b19d)

Sev	Finding	Evidence	Status
MED	the ksponly guard read snes.getType() — the previous solve's type — not the declared option	set ksponly after a solve, cold solve: picard = 1 ran	fixed: reads petsc_options
MED	"a cold start takes one Picard step" is true only on the saddle-point standard path — not scalar/vector solvers, not the rotated/fault-contact path	Poisson: no warm-up; rotated: log begins at Newton	wording fixed (docstring, CLAUDE.md)
MED	the warm-up is an nrichardson residual sweep, not a Picard step: no linear solve, no frozen tangent; 1–12% residual reduction on linear Stokes; the sharp-yield case converges under pure Newton without it (7 its)	snes_monitor; 1057 fixture on the rotated path	comment corrected; picard=-1 switches it off; dropping it is a planning benchmark item
LOW	Stokes_Constrained 2-D annulus with two add_constraint_bc(degree=2) crashes in _scatter_global_to_fields (pre-existing)	broadcast (12484,) into (6396,)	planning item; own issue

Attacks that failed: linear Stokes newtonls+EW new vs old default identical to 1e-9 on box and annulus; rotated + sharp yield cold: Newton 7 its where Picard did not converge in 50; Poisson κ(T) cold: Newton 8–10 its vs Picard 24–28, no window where the old default converged and the new one does not; continuation unchanged; the homotopy march forces its own tangent.

Also from self-review: the Stokes adjoint verdict called the two-assembly nonlinearity probe on every recorded solve; now the symbolic test, cached until the residual can change.

[Copilot]

Copilot was unable to review this pull request because the user who requested the review has reached their quota limit.

[lmoresi]

Adversarial review: the merge with development (4649653)

Thirteen hunks across six files. Ten were placement collisions — this branch
adds adjoint members where development added unrelated ones — and keep both
sides. Three were decisions, and taking both sides would have been wrong in
each.

Three places where keeping both sides was the wrong answer

_jacobian_unwrap. Both sides define it. development's signature is
(expr, route=None) for the #836 shared-graph route; this branch's is the
one-argument version from before #823 tier 2. This branch's copy is dropped.
Checked before dropping it: the only caller in the solver passes the route
explicitly (_jacobian_unwrap(expr, route=self._jit_route_in_use())), and the
only one-argument callers are in development's own
tests/test_1066_stokes_jacobian_layout.py, which route=None serves. No
adjoint code calls it. This matters for the adjoint specifically: the adjoint
must transpose the kernels the forward assembly actually produced, so it has
to use whichever route the forward solve used, not a second unwrap of its own.

The cold-start warm-up. This branch carries a hand-rolled nrichardson
sweep: declared = ...; if (picard == 0 and consistent_jacobian is True and declared != "ksponly" and zero_init_guess): picard = 1. development has
_plan_picard_warmup / _snes_solve_with_retries (#794 — Newton default with
one automatic frozen-tangent step). These are two implementations of one
policy. Keeping both would have applied it twice. development's stands and
this branch's block is dropped.

The consistent_jacobian docstring. development's text states the
current default and this branch's does not (it still says "flipped 2026-09").
development's is kept, plus one sentence it does not say and the adjoint
needs: the consistent tangent is the matrix the adjoint transposes, so an
adjoint solve inherits the choice made here. On the first attempt that
sentence landed between the True and "continuation" entries of the
definition list, breaking it; it is now a paragraph after the list.

Two failures the merge brought to light, both tests working

test_every_history_scheme_declares_a_verdict caught a real gap.
ForwardSemiLagrangian is a scheme development added after this branch
forked, and it declared no adjoint verdict — exactly what that test exists to
catch, since the base class refuses by naming the class so an undeclared
scheme cannot pass as either. It now refuses, with a reason read off its own
mechanism rather than copied from the backward scheme: the forward trace is
differentiable in the velocity, but the per-cell least-squares fit of the
arrivals is not materialised as an operator, and which cell an arrival lands
in — and so which fit reads it — is piecewise constant in the velocity.

test_a_swarm_step_that_loses_particles_refuses_with_the_count had stopped
testing anything. It set return_coords_to_bounds = None to switch the
restore off. On development None selects the general facet restore instead
(the deliberate fix for file meshes, which used to restore nothing at all), so
every particle was returned: measured 972 -> 972, with every particle clamped
to x = 0.99994. The assertion n_after < n_before failed, and would have gone
on failing for a reason that has nothing to do with the adjoint.

Switched off is now an identity closure. The step is also made partial
(972 -> 482 at dt = 0.05, measured) rather than total: with the original
dt = 0.5 an identity closure loses every particle (972 -> 0), and a refusal
read off an empty swarm is not the dimension change the test is about.

This is the fixture-no-longer-shows-the-effect failure mode, and it is worth
naming: the test was green when written, went silently void when a default
changed elsewhere, and only a merge surfaced it.

CI reach

test_1063_constrained_traction::test_constrained_topography_carries_the_augmentation_share
was the one real failure on this PR's last CI run (2026-09-24, SNES reason
-3). It passes on the merge.

What was run, and what was not

Run locally, at HEAD, after rm -rf build/ and ./uw build (a .pyx change):

• the eight adjoint test files — 36 passed
• test_1066_stokes_jacobian_layout, test_1067_newton_cold_start,
  test_1068_picard_warmup_is_a_frozen_tangent_step — 28 passed. These
  are the files that target the two dropped blocks directly: 1066 calls
  _jacobian_unwrap itself, 1067/1068 pin the warm-up policy.
• test_0030_capability_guides, test_0641_wave_c_api_shims,
  test_1057_yield_homotopy_solve (the conflicted test files) — 67 passed
• test_1063_constrained_traction — 1 passed

Not run locally: the full level_1 and tier_a suite. Three attempts
produced no output and exit 59, which is not a pytest exit code, and the cause
is not diagnosed. It is left to CI rather than reported as a pass. The
targeted runs above cover the three files this merge actually decided
something in; the broad regression check on the rest of the .pyx is CI's.

Issue state

No issue is marked on this merge. The adjoint issues this branch gates (#751
is stacked on it; #763, #766, #769, #773, #787, #804) are not probed here
because they need the branch merged to development before a probe against
development means anything. Per CLAUDE.md they stay as they are until then —
marking them on the word of a merge rather than a probe is what the policy
forbids.

Underworld development team with AI support from Claude Code

[lmoresi]

Correction and a second finding: the full tier_a suite did run

Two things in the review above are now wrong, and the second is a defect in
this branch.

The suite ran: 1449 passed, 1 failed

The review said the level_1 and tier_a suite produced no output and exit 59
with the cause undiagnosed. The cause is diagnosed and it was mine, not the
code's: setup.cfg's [tool:pytest] sets no testpaths, so
pytest -m "level_1 and tier_a" with no path argument collects from the
rootdir — which in a worktree includes .pixi/envs/ and its thousands of
vendored site-packages tests. Hence no output, a very long run and an exit
code that is not one of pytest's.

With the path given, pytest tests -m "level_1 and tier_a":

1 failed, 1449 passed, 56 skipped, 2116 deselected, 1 xfailed in 978.98s

CI never meets this because scripts/test.sh always passes explicit globs,
never a bare marker selection. A one-line testpaths = tests would retire it;
it is not this branch's territory, so it is left alone and noted here.

The one failure was real: Charter S7

FAILED tests/test_0023_field_realness_and_derivatives.py::
       test_no_solver_differentiates_with_plain_sympy_diff

development has zero sympy.diff( in
petsc_generic_snes_solvers.pyx. All five sites were this branch's adjoint
code, and the census test that forbids them landed on development after
this branch forked — so nothing on either side alone would have caught it.

It is not a style rule. sympy.diff differentiates with respect to anything
that is not a Symbol by swapping in a Dummy carrying no assumptions,
so the field is complex while it is differentiated: Abs(u - c) yields
sign(u - c) * du/du with the last factor left unevaluated, which cannot be
compiled. diff_wrt_field gives the stand-in the field's realness and is
otherwise the same swap.

Two of the five sites were exposed to exactly that. adjoint_integrand(mu, wrt) and the Dirichlet-datum sensitivity both differentiate with respect to
wrt, and both docstrings say wrt may be a field — "for a field, project
it to get the dual on that field". A yield law putting an Abs or a sign of
a field in its flux would have reached them.

Fixed in 33a4b198. Four sites are drop-in replacements. The fifth, in the
nonlinearity probe, differentiates twice by a Symbol the method constructs
itself, where diff_wrt_field delegates straight to sympy.diff; it is
converted anyway so that every derivative in the solvers reads the same and
the census stays a blunt textual check with nothing to audit.

Rebuilt with rm -rf build/ first, as a .pyx change needs. test_0023
plus the 36 adjoint tests pass — 45 in all.

Where that leaves the branch

• full level_1 and tier_a: 1449 passed, 1 failed, and that failure is now
  fixed
• the 36 adjoint tests, test_0023, test_1066/1067/1068,
  test_0030_capability_guides, test_0641, test_1057,
  test_1063_constrained_traction: all pass
• The discrete adjoint of a rotated free-slip solve #751, stacked on this branch, came through the push MERGEABLE/CLEAN with
  CI passing

The tier_a suite has not been re-run end to end since 33a4b198; the
targeted runs above cover what it changed, and CI has the rest.

Underworld development team with AI support from Claude Code