feat(Rosenbrock): unified RosenbrockAverage substepping framework + 2M+P3 exact-AD mode

[haakon-e]

What this adds

Introduces the RosenbrockAverage{Jacobian, GrowthTreatment, TendencyLimiter} options framework for averaged-tendency microphysics substepping, with its first use: the 2-moment warm-rain + P3 ice (Microphysics2Moment) exact-AD mode.

Each call divides Δt into nsub linearized-implicit (Rosenbrock-Euler) substeps; per substep it solves (I/h − J) Δ = f(x), clamps the updated state to nonnegative values, and advances temperature from the latent heat of the realized increment. f is the instantaneous tendency and J is the ForwardDiff Jacobian of f(x); logλ and q_tot are held fixed across the interval.

The mode is configured by:

• Jacobian, the substep matrix. Options: ExactJacobian, the exact ForwardDiff derivative. (DonorJacobian and CoupledDonorJacobian, the donor-based linearizations, are 1M-only and arrive in feat(Rosenbrock): 1M donor / coupled / exact modes on the unified framework #731.)
• GrowthTreatment, how the positive growth diagonal enters the operator. Options: ImplicitGrowth; ExplicitGrowthDiagonal, which zeros the positive growth diagonal so a growth mode is taken explicitly.
• TendencyLimiter, a limiter on the realized increment. Options: NoLimiter; EndStateSaturationAdjustment, which scales the increment so the latent-heated end state stays at or above ice saturation.

The preset rosenbrock_exact() = ExactJacobian + ExplicitGrowthDiagonal + EndStateSaturationAdjustment is the 2M+P3 mode. On 2M+P3 only ExactJacobian is supported.

Coarse-step ice-growth instability

At cold, ice-supersaturated states carrying supercooled cloud liquid, the ice-growth tendency has an autocatalytic mode: rime mass grows by collecting cloud droplets, and the exact Jacobian carries this as a positive diagonal in the rime-mass (q_rim) row. At a representative isolated cell — ρ = 1.0 kg m⁻³, T = 263 K, q_tot = 10⁻², q_lcl = 2 × 10⁻³, n_lcl = 10⁸ m⁻³, q_ice = 2 × 10⁻³, n_ice = 10⁴ m⁻³, unrimed, ice supersaturation S_ice ≈ 1.8 — this diagonal is +5 × 10⁻² s⁻¹ (time scale ≈ 20 s, identical in Float32 and Float64), and it grows past 10⁻¹ s⁻¹ at colder, more liquid-rich states. The mode requires supercooled liquid: with the same ice state but no cloud liquid the rime-mass diagonal falls to order 10⁻⁵ s⁻¹, and the pure-deposition diagonal is negative there (vapor depletion opposes further deposition).

With ImplicitGrowth the implicit operator loses definiteness once the substep is coarse relative to the growth time scale, and a single coarse step overshoots saturation. ExplicitGrowthDiagonal removes that mode from the implicit operator and EndStateSaturationAdjustment supplies the nonlinear saturation bound; together they make the exact scheme robust across the resolved range of time steps. The crash itself is a property of the single-column convective configuration, not of an isolated cell: the growth diagonal above is an isolated-cell measurement, but the non-physical blow-up reproduces only in the coupled single-column run. docs/src/RosenbrockNumerics.md documents the mechanism, the resolution, the API and how to extend it, and the experiments that did not resolve it.

Conditioning

Per substep, a diagonal equilibration brings the number and mass rows (which differ by about nine orders of magnitude) to a comparable scale before the solve, and max(x + Δ, 0) is applied afterward; a non-finite Jacobian falls back to a forward-Euler substep. For the ImplicitGrowth configuration, a mask can advance near-empty species by forward Euler rather than including them in the implicit solve (the tendency gradients are very large at near-zero mass, which degrades the conditioning and accuracy of the solve).

Float32 robustness

The exact ForwardDiff Jacobian requires the tendency's derivatives to stay finite in single precision, where two sites did not:

• sgs_weight_function (the rime-density blend): the sigmoid argument u can round to ±1 within the transition band, so atanh(±1) = ±Inf returns a finite value with a NaN derivative at rimed pure-ice states. This is chiefly a single-precision hazard, but the same band-edge rounding also occurs in Float64. The branch now uses the identity tanh(2·atanh(u)) = 2u/(1 + u²), which differentiates finitely across the band and matches the original value to within one to two ULP.
• gamma_inc_moment: given an analytic ForwardDiff derivative (∂/∂α = −gamma_inc_moment(·, p+1, α), ∂/∂D = D^p e^{−αD}), so the closed-form rain integrals differentiate through the exact integrand rather than through a difference of near-equal incomplete gammas, which loses precision on the narrow subintervals from the fall-speed-crossing split. Values are unchanged.

Testing

Tests of the 2M+P3 exact mode: finite tendencies across Δt × nsub in Float32 and Float64; a non-ExactJacobian RosenbrockAverage and warm-rain-only (with_ice = false) parameters both throw; the rosenbrock_exact() bulk tendencies broadcast over ClimaCore fields; and an assertion that the exact 2M+P3 substep is type-stable (@inferred, JET @test_opt) and allocation-free. The exact-AD 2M+P3 substep is type-stable and allocation-free on the CI Julia versions; the allocation check uses an interpolated BenchmarkTools benchmark (as in bench_press), which measures the call's own allocations rather than a surrounding closure, so the assertion runs unconditionally.

Foundation for the 1M modes (#731) and the verbose diagnostics (#732).

[haakon-e]

This change is part of the following stack:

• feat(P3): Gauss-Legendre quadrature, integrand breakpoints -> GL(6) default #741
  • feat(Rosenbrock): unified RosenbrockAverage substepping framework + 2M+P3 exact-AD mode #726 ◀
    • feat(Rosenbrock): 1M donor / coupled / exact modes on the unified framework #731
      • feat(Rosenbrock): post-solve per-process verbose diagnostics #732
        • fix(Rosenbrock): saturation adjustment on max(S_ice, S_liq) #737
          • perf(BMT): recover the primal tendency from the Rosenbrock Jacobian pass #742
            • feat(BMT): analytic ManualJacobian for the 2M+P3 Rosenbrock substep #743
              • feat(BMT): coupled Tier-3 donor linearization for the 2M+P3 ManualJacobian #744

_{Change managed by git-spice.}

[codecov]

Codecov Report

❌ Patch coverage is 88.49558% with 13 lines in your changes missing coverage. Please review.
✅ Project coverage is 92.65%. Comparing base (355b541) to head (8ae90c4).

Additional details and impacted files

@@                     Coverage Diff                     @@
##           he/ice-selfcol-triangle     #726      +/-   ##
===========================================================
- Coverage                    92.81%   92.65%   -0.16%     
===========================================================
  Files                           56       57       +1     
  Lines                         2952     3065     +113     
===========================================================
+ Hits                          2740     2840     +100     
- Misses                         212      225      +13     

Components	Coverage Δ
src	93.40% <88.49%> (-0.20%)	⬇️
ext	69.47% <ø> (ø)

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.