perf(BMT): recover the primal tendency from the Rosenbrock Jacobian pass

ForwardDiff's Jacobian pass already evaluates the primal tendency f as the dual
.value and discards it, so a separate f = g(x) repeats a full (quadrature-
dominated) tendency evaluation. Recover f from the same pass: seed all N
partials of the state in a single call to g with static Duals, and read the
values and partials directly into the same FieldVector species type as x
(MicroState1M or MicroState2MP3). Applies to the 2M+P3 and the 1M
ExactJacobian paths; the analytic Donor/CoupledDonor 1M Jacobians (no f
by-product) keep their separate evaluation.

The DiffResults dependency is dropped: DiffResults.JacobianResult's mutable
result buffer heap-allocates on GPU, so the seeding is done directly with
StaticArrays instead. Both paths go through the shared
_tendency_and_jacobian(jacobian, g, x) method, removing the duplicated code
in the 2M+P3 substep loop.

Measured on MicroState2MP3 (N=8, quadrature-dominated) and MicroState1M
(N=4): the single full-width call is ~31% faster than the previous two-call
f = g(x); J = FD.jacobian(g, x) for 2M+P3, and ~39% faster for 1M, with zero
allocations in both cases. f and J are bitwise identical to the two-call
reference across representative states. The tendency and Jacobian are
unchanged, the path stays type-stable, and it composes with an outer AD pass.

Author: haakon-e

src/BMT_rosenbrock.jl

@@ -438,10 +438,10 @@ fields as the `Instantaneous` entry (without the activation diagnostic).
     Tsub = T
     for _ in 1:nsub_eff
         g = Instantaneous2MP3Tendency(mp, tps, ρ, Tsub, q_tot, logλ)
-        f = g(x)
         x_prev = x
         if all(isfinite, x)
-            J = _apply_growth(mode.growth, FD.jacobian(g, x))
+            f, J_raw = _tendency_and_jacobian(mode.jacobian, g, x)
+            J = _apply_growth(mode.growth, J_raw)
             z = _species_mask(mode.jacobian, mode.growth)(x)
             d = if all(isfinite, J)
                 _rosenbrock_update(x, f, J, z, h) - x
@@ -451,6 +451,7 @@ fields as the `Instantaneous` entry (without the activation diagnostic).
             d = _apply_limiter(mode.limiter, x, d, ρ, Tsub, q_tot, Lv_over_cp, Ls_over_cp, tps)
             x = max.(x .+ d, 0)
         else
+            f = g(x)
             x = _euler_update(x, f, h)
         end
         Δ = x - x_prev
@@ -756,6 +757,33 @@ end
 "Exact ForwardDiff Jacobian provider for [`_rosenbrock_average_1m`](@ref)."
 @inline _ad_jacobian_1m(g, x) = FD.jacobian(g, x)
 
+"""
+    _tendency_and_jacobian(jacobian, g, x)
+
+The raw substep tendency `f = g(x)` and the substep Jacobian (before the growth
+treatment) for a [`Jacobian`](@ref) option, returned as `(f, J)`.
+
+For [`ExactJacobian`](@ref) the primal `f` and the `N×N` Jacobian are both
+obtained from `ForwardDiff`. The donor-based matrices ([`DonorJacobian`](@ref),
+[`CoupledDonorJacobian`](@ref)) produce no tendency by-product, so `f = g(x)` is
+evaluated separately.
+"""
+@inline function _tendency_and_jacobian(::ExactJacobian, g, x::SA.FieldVector{N, FT}) where {N, FT}
+    Tag = typeof(FD.Tag(g, FT))
+    dx = SA.SVector(
+        ntuple(i -> FD.Dual{Tag}(x[i], ntuple(s -> ifelse(s == i, one(FT), zero(FT)), Val(N))...), Val(N)),
+    )
+    y = g(dx)
+    f = typeof(x)(ntuple(i -> @inbounds(FD.value(y[i])), Val(N))...)
+    J = SA.SMatrix{N, N, FT}(
+        ntuple(k -> @inbounds(FD.partials(y[(k - 1) % N + 1], (k - 1) ÷ N + 1)), Val(N * N)),
+    )
+    return f, J
+end
+@inline _tendency_and_jacobian(::DonorJacobian, g, x) = (g(x), _jacobian_1m_linearized(g, x))
+@inline _tendency_and_jacobian(::CoupledDonorJacobian, g, x) =
+    (g(x), _jacobian_1m_coupled(g, x))
+
 """
     _full_species_mask(x)
 
@@ -784,18 +812,17 @@ and the increment limiter through [`_jacobian_provider`](@ref),
     h = Δt / FT(nsub_eff)
     Lv_over_cp = TDI.TD.Parameters.LH_v0(tps) / TDI.TD.Parameters.cp_d(tps)
     Ls_over_cp = TDI.TD.Parameters.LH_s0(tps) / TDI.TD.Parameters.cp_d(tps)
-    jacobian = _jacobian_provider(mode.jacobian)
     mask = _species_mask(mode.jacobian, mode.growth)
 
     x = MicroState1M{FT}(q_lcl, q_icl, q_rai, q_sno)
     x₀ = x
     Tsub = T
     for _ in 1:nsub_eff
         g = Raw1MTendency(mp, tps, ρ, Tsub, q_tot)
-        f = g(x)
         x_prev = x
         if all(isfinite, x)
-            J = _apply_growth(mode.growth, jacobian(g, x))
+            f, J_raw = _tendency_and_jacobian(mode.jacobian, g, x)
+            J = _apply_growth(mode.growth, J_raw)
             z = mask(x)
             d = if all(isfinite, J)
                 _rosenbrock_update(x, f, J, z, h) - x
@@ -805,6 +832,7 @@ and the increment limiter through [`_jacobian_provider`](@ref),
             d = _apply_limiter(mode.limiter, x, d, ρ, Tsub, q_tot, Lv_over_cp, Ls_over_cp, tps)
             x = max.(x .+ d, 0)
         else
+            f = g(x)
             x = _euler_update(x, f, h)
         end
         Δ = x - x_prev