solver.sensitivity() returns dJ/dm in SI, not in the non-dimensional system: wrong by the reference scales #763

What

solver.sensitivity(mu, wrt) assembles mu^T dR/dm as

total = float(uw.maths.Integral(self.mesh, self.adjoint_integrand(mu, wrt)).evaluate())

Integral.evaluate() re-dimensionalises: PETSc integrates in the model's
non-dimensional space and evaluate then attaches
integrand_units * coordinate_units**dim and converts, returning a UWQuantity.
float() on that takes the magnitude in those units, while mu, the residual
templates and the parameter value the chain rule differentiates against are all
non-dimensional. The returned number therefore carries whatever reference scales the
surviving units name.

Two separate leaks:

1. A dimensionless knob multiplying a dimensional scale. _peel_except deliberately
   keeps dimensional values as leaves (substituting them raises SympifyError). Where
   the parameter is dimensional this is harmless, because it is differentiated away; but
   a dimensionless parameter leaves the dimensional factor standing in dR/dm, so the
   integral comes back in e.g. km**2 * Pa * s.
2. The coordinate scale, always. The integral carries coordinate_units**dim
   regardless of the parameter, so any model whose length reference is not 1 is off by
   L_ref**dim in every sensitivity. A 30 km reference is a silent factor 900.

Measured

Linear Stokes is an exact oracle: with the viscosity constant and the body force
independent of it, u ~ 1/eta, so for any J linear in u, dJ/d(eta) = -J/eta.

Rig: unit box, shear_viscosity_0 = alpha * eta with alpha dimensionless = 2 and
eta = 1e22 Pa*s; viscosity reference 1e21 Pa*s, length reference 1 km.

1.000000e+21 is exactly the viscosity reference. The manifest is not at fault —
_pack_constants packs \\eta -> 1.000000e+01, correctly non-dimensional. The integral
itself returns 207999810040165.88 [kilometer ** 2 * pascal * second], and
uw.non_dimensionalise of that gives 2.07999810e-07, the right answer.

Where it showed up

The Spiegelman notch viscosity floor, shear_viscosity_min = m * floor * eta_bg with
floor dimensionless. dJ/d(floor) read 3.25185e+23 against a central finite
difference of -4.05e-01. eta_bg is the viscosity reference in that model
(1e24 Pa*s), which is why the corruption looked like a plausible physical scale rather
than a bug.

Sites

• petsc_generic_snes_solvers.pyx sensitivity() — volume term and the natural-BC facet term
• adjoint.py integral(), and TranscriptAdjoint.gradient()'s J and explicit dJ/dm

Fix

Integral.evaluate() re-dimensionalising is documented behaviour that other callers
depend on, so the conversion belongs at the adjoint boundary: uw.adjoint.nd_float(),
used at every site above. Regression test in
tests/test_0025_adjoint_sensitivity_is_non_dimensional.py, asserting against the closed
form with the length reference deliberately 30 km and the viscosity reference
deliberately different from the viscosity, so neither scale can cancel by accident.