adjoint_solve inherits the forward solve's INEXACT Newton tolerance: every sensitivity on an iteratively-solved problem is silently short by 15-20% #769

What

adjoint_solve reuses the SNES's own KSP:

ksp = self.snes.getKSP()
ksp.setOperators(J, P)
ksp.solve(b, x)

so it inherits the forward problem's linear tolerance. A Newton solve's linear solves
are deliberately inexact — the Stokes tolerance property enables Eisenstat–Walker
(snes_ksp_ew), which re-picks the outer ksp_rtol every step and starts from PETSc's
default rtol0 = 0.3. That is correct for a forward step, because the outer Newton
iteration corrects the error.

mu has no outer iteration. The adjoint solve's error goes straight into
sensitivity(), as a fixed fraction — independent of any finite-difference step used
to check it, which is exactly what makes it look like a modelling result rather than a bug.

The KSP still reports CONVERGED, because it did converge — to 0.3.

Measured

Shear-thinning Stokes (eta = eta_0/(1 + edot/edot_0)), 4 Newton iterations, J = 1/2 int edot_II^2, dJ/d(eta_0) against central finite differences. Forced onto an iterative KSP
(fgmres, pc_type none) so the tolerance binds:

On the Spiegelman notch's viscosity floor (fieldsplit-Schur multigrid, 1e3 viscosity
contrast) the same defect reads 0.8088 / 0.8062 / 0.8040 at FD steps of 10 / 3 / 1 per
cent — step-independent, as predicted.

Why no test caught it

A small problem's KSP is DIRECT — 1 iteration, residual 1e-15 — so no tolerance binds
and the adjoint is exact whatever the setting. Every existing adjoint test, and five
purpose-built reproducers (shear-thinning; smooth_max floor; pressure-dependent
viscosity; Drucker–Prager + floor; P0 discontinuous pressure) all returned ratios of
0.9986–1.0144 for this reason. The defect needs an iteratively solved problem, which is
every problem anyone actually runs.

And the user cannot work around it

Setting ksp_rtol does not stick, because Eisenstat–Walker overwrites it every step. Five
runs requesting rtol between 1e-12 and 3e-1 all reported rtol 3.00e-01. The tolerance
property documents this for the forward path ("to steer the linear solve via ksp_rtol you
must first switch snes_ksp_ew off"); nothing said it applies to the adjoint, where it is
not a tuning preference but a correctness requirement.

Fix

adjoint_solve(rhs, target=None, rtol=None) now pins the transpose solve:
snes.setUseEW(False), ksp.setTolerances(rtol=_ADJOINT_KSP_RTOL, atol=1e-50, max_it>=2000), and restores the KSP's saved tolerances and EW state in a finally — the
KSP belongs to the SNES and is reused by the next forward solve, so leaving it pinned would
silently make every subsequent Newton step solve far harder than the inexact-Newton design
intends. Default _ADJOINT_KSP_RTOL = 1e-10, overridable per call.

After the fix, the same iterative rig gives 0.999078 with the forward tolerance left at
0.3, and the direct-solve path is unchanged.

Note for anyone who has used this API

Any sensitivity or inversion gradient taken on an iteratively-solved problem before this
fix is short by a fixed fraction — order 15-20% on the cases measured. Gradient directions
are much less affected than magnitudes, so an inversion may well have converged anyway,
but Taylor tests and FD cross-checks on real problems would have failed at the few-per-cent
level and been attributed to discretisation.