Differentiate raw 2M+P3 tendency in RosenbrockAverage

The Rosenbrock mode previously linearized a limited instantaneous
tendency: a supersaturation cap and 1/dt coupled mass-number sink
limits were applied before the ForwardDiff Jacobian was taken. Those
clips are explicit-Euler stabilization devices, not model physics, and
they do not belong inside an implicit update.

The published P3 condensation/deposition scheme (Morrison and Milbrandt
2015, J. Atmos. Sci., 72, 287-311, doi:10.1175/JAS-D-14-0065.1,
Appendix C) is an analytic time-averaged relaxation with no tendency
clip, so the raw relaxation ODE is the model physics. The one-stage
Rosenbrock update is L-stable and damps the stiff vapor-exchange
subsystem monotonically, so the saturation overshoot and ringing those
limiters suppress cannot occur. The clips are also harmful: 1/dt
tendency limits inject dt-independent error and break convergence under
refinement (Wan et al. 2020, J. Adv. Model. Earth Syst., 12,
e2019MS001982), and a hard saturation clip structurally forbids the
mixed-phase quasi-steady vapor pressure, which lies between liquid and
ice saturation (Korolev and Mazin 2003, J. Atmos. Sci., 60, 2957-2974).

The differentiated functor now evaluates the unmodified
Microphysics2Moment process rates. The discrete safeguards that remain
are h-free conditioning and projection of the linear solve, not model
terms: the channel mask, equilibration, and positivity clamp. The cap
and sink helpers, used only by the Rosenbrock path, are removed.

Convergence-test references are rebuilt against the raw tendency and
tolerances re-derived from measured errors. The rime-volume and
rime-mass error floors are raised to genesis scale: once ice has
melted, a single-precision residual in the collapsing coupled
(q_rim, b_rim) pair is physically negligible but unbounded against an
eps-scale floor.

Author: haakon-e

docs/bibliography.bib

@@ -819,4 +819,15 @@ @article{Horn2012
   pages={2635--2660},
   year={2011},
   publisher={G{\"o}ttingen, Germany}
+}
+
+@article{Wan2020,
+  author = {Wan, Hui and Zhang, Shixuan and Rasch, Philip J. and Larson, Vincent E. and Zeng, Xubin and Yan, Huiping},
+  title = {Quantifying and attributing time step sensitivities in present-day climate simulations conducted with EAMv1},
+  journal = {Journal of Advances in Modeling Earth Systems},
+  year = {2020},
+  volume = {12},
+  number = {12},
+  doi = {10.1029/2019MS001982},
+  pages = {e2019MS001982}
 }

docs/src/BulkTendencies.md

@@ -162,11 +162,14 @@ For the 2-moment + P3 configuration, `RosenbrockAverage` replaces the hand-built
 \left(\frac{I}{h} - P J P\right) \Delta x = f(x), \qquad x \leftarrow \max(x + \Delta x,\, 0)
 ```
 
-where ``f`` is the *limited* instantaneous tendency (the supersaturation cap and the coupled mass–number sink limits are inside the differentiated function, so the implicit update respects them), and ``J = \partial f / \partial x`` is its exact 8×8 `ForwardDiff` Jacobian.
+where ``f`` is the **raw instantaneous tendency** — the unmodified `Microphysics2Moment` process rates, with no timestep-dependent clipping — and ``J = \partial f / \partial x`` is its exact 8×8 `ForwardDiff` Jacobian.
 
-Two structural safeguards:
+The differentiated tendency is the model physics, not a stabilized variant of it. The P3 condensation/deposition scheme is an analytic time-averaged relaxation [MorrisonMilbrandt2015](@cite) with no tendency clip; a supersaturation cap and ``1/h`` sink limits are explicit-Euler stabilization devices, not physical terms. They are unnecessary here because the one-stage Rosenbrock update is L-stable: it damps the stiff vapor-exchange subsystem monotonically, so the saturation overshoot and ringing those limiters suppress cannot occur. They are also harmful: ``1/h`` tendency clips inject ``h``-independent error and break convergence under refinement [Wan2020](@cite), and a hard saturation clip structurally forbids the mixed-phase quasi-steady vapor pressure, which lies between liquid and ice saturation [KorolevMazin2003](@cite).
 
-- **Channel projection** ``P = \mathrm{Diag}(z)``: channels whose condensed mass is below ``10^{-10}`` are projected out of the Jacobian. Their rows of ``I/h - PJP`` reduce to the identity, so the solve returns exactly a forward-Euler update for those channels while healthy channels stay implicit — an IMEX-style splitting at channel granularity. Near-empty channels otherwise produce finite but very large Jacobian entries (through the sink limiter) whose linearized steady state fabricates phantom number concentrations that substep refinement cannot remove.
+The discrete safeguards are therefore conditioning and projection devices, all ``h``-free and applied to the linear solve rather than to the physics:
+
+- **Channel projection** ``P = \mathrm{Diag}(z)``: channels whose condensed mass is below ``10^{-10}`` are projected out of the Jacobian. Their rows of ``I/h - PJP`` reduce to the identity, so the solve returns exactly a forward-Euler update for those channels while healthy channels stay implicit — an IMEX-style splitting at channel granularity. Near-empty channels otherwise produce finite but very large Jacobian entries whose linearized steady state fabricates phantom number concentrations that substep refinement cannot remove.
 - **Equilibration** ``S = \mathrm{Diag}(|x| + h|f| + \epsilon)``: the linear system is solved as ``S^{-1} A S`` so the rows, which span roughly nine orders of magnitude across number and mass species, become O(1)-conditioned. This keeps single-precision roundoff relative to each species' own scale; an unscaled Float32 factorization deposits roundoff from the large rows into empty species as phantom mass.
+- **Positivity clamp** ``x \leftarrow \max(x + \Delta x, 0)``: a projection onto the physical nonnegative orthant after each substep.
 
-The local temperature is advanced each substep from the latent heating of the realized increments. `logλ` and `q_tot` are held fixed across the interval, matching the explicit-substepping semantics; non-finite states or Jacobians fall back to forward-Euler substeps of the limited tendency. The implicit update makes the stiff ice-process path insensitive to the substep length; at very large substeps (``h \gtrsim 100`` s) the limiter switching degrades the single linearization, so increase `nsub` rather than relying on the limiters.
+The local temperature is advanced each substep from the latent heating of the realized increments. `logλ` and `q_tot` are held fixed across the interval, matching the explicit-substepping semantics; non-finite states or Jacobians fall back to forward-Euler substeps of the raw tendency. The implicit update makes the stiff ice-process path insensitive to the substep length; at very large substeps (``h \gtrsim 100`` s) the single linearization carries the usual first-order error of a one-stage method, so increase `nsub` to refine.

src/BMT_rosenbrock.jl

@@ -23,140 +23,41 @@ SA.similar_type(::Type{<:MicroState2MP3}, ::Type{FT}, ::SA.Size{(8,)}) where {FT
     MicroState2MP3{FT}
 
 """
-    _sink_limit(q, dt, n = 3)
-
-Largest admissible sink rate for `q`: deplete at most `1/n` of it over `dt`.
-"""
-@inline _sink_limit(q, dt, n = 3) = max(zero(q), q) / dt / n
-
-"""
-    _coupled_sink_factor(S1, S2, q1, q2, dt)
-
-Common scaling factor in `(0, 1]` for a coupled (mass, number) sink pair
-`(S1, S2)`, such that neither `f S1` nor `f S2` exceeds its
-[`_sink_limit`](@ref) over `dt`. Scaling the pair by one common factor
-preserves the mean particle mass implied by the pair; limiting the members
-independently would not.
-"""
-@inline function _coupled_sink_factor(S1, S2, q1, q2, dt)
-    M1 = _sink_limit(q1, dt)
-    M2 = _sink_limit(q2, dt)
-    f1 = ifelse((S1 < zero(S1)) & (-S1 > M1), M1 / -S1, one(S1))
-    f2 = ifelse((S2 < zero(S2)) & (-S2 > M2), M2 / -S2, one(S2))
-    return min(f1, f2)
-end
-
-"""
-    _cap_supersaturation(tend, tps, T, ρ, q_tot, q_lcl, q_rai, q_ice, dt)
-
-Cap the net condensation/deposition tendencies in `tend` so that one substep of
-length `dt` cannot overshoot saturation: the liquid channel
-(`dq_lcl + dq_rai`) is capped against the analytic saturation-adjustment
-condensation increment, then the ice channel (`dq_ice + dq_rim`) against the
-deposition increment computed from the vapor remaining after the liquid
-step. Each channel is scaled by a common, sign-preserving ratio in `[0, 1]`
-(both members of a channel are scaled together, like
-[`_coupled_sink_factor`](@ref)). Branch-free: `ifelse` evaluates both arms,
-which is safe here and keeps the kernel GPU-friendly.
-"""
-@inline function _cap_supersaturation(tend, tps, T, ρ, q_tot, q_lcl, q_rai, q_ice, dt)
-    (; dq_lcl_dt, dn_lcl_dt, dq_rai_dt, dn_rai_dt,
-        dq_ice_dt, dn_ice_dt, dq_rim_dt, db_rim_dt) = tend
-    FT = typeof(q_tot)
-    T_safe = max(150, T)
-    q_vap = max(0, q_tot - q_lcl - q_rai - q_ice)
-
-    Rᵥ = TDI.Rᵥ(tps)
-    cp_d = TDI.TD.Parameters.cp_d(tps)
-    L_v = TDI.Lᵥ(tps, T_safe)
-    L_s = TDI.Lₛ(tps, T_safe)
-    qv_sat_liq = TDI.saturation_vapor_specific_content_over_liquid(tps, T_safe, ρ)
-    qv_sat_ice = TDI.saturation_vapor_specific_content_over_ice(tps, T_safe, ρ)
-
-    # liquid channel vs the saturation-adjustment condensation increment
-    qcon_satadj =
-        (q_vap - qv_sat_liq) / (1 + L_v^2 * qv_sat_liq / (cp_d * Rᵥ * T_safe^2)) / dt
-    net_liq = dq_lcl_dt + dq_rai_dt
-    target_liq = ifelse(
-        net_liq > 0, min(net_liq, max(zero(qcon_satadj), qcon_satadj)),
-        ifelse(net_liq < 0, max(net_liq, min(zero(qcon_satadj), qcon_satadj)), zero(net_liq)),
-    )
-    ratio_liq = clamp(
-        ifelse(abs(net_liq) > eps(FT), target_liq / net_liq, one(net_liq)),
-        zero(net_liq), one(net_liq),
-    )
-    dq_lcl_dt *= ratio_liq
-    dq_rai_dt *= ratio_liq
-    dn_lcl_dt *= ratio_liq
-    dn_rai_dt *= ratio_liq
-
-    # ice channel vs deposition of the vapor remaining after the liquid step
-    qv_after_liq = max(0, q_vap - target_liq * dt)
-    qdep_satadj =
-        (qv_after_liq - qv_sat_ice) /
-        (1 + L_s^2 * qv_sat_ice / (cp_d * Rᵥ * T_safe^2)) / dt
-    net_ice = dq_ice_dt + dq_rim_dt
-    target_ice = ifelse(
-        net_ice > 0, min(net_ice, max(zero(qdep_satadj), qdep_satadj)),
-        ifelse(net_ice < 0, max(net_ice, min(zero(qdep_satadj), qdep_satadj)), zero(net_ice)),
-    )
-    ratio_ice = clamp(
-        ifelse(abs(net_ice) > eps(FT), target_ice / net_ice, one(net_ice)),
-        zero(net_ice), one(net_ice),
-    )
-    dq_ice_dt *= ratio_ice
-    dq_rim_dt *= ratio_ice
-    dn_ice_dt *= ratio_ice
-    db_rim_dt *= ratio_ice
-
-    return (; dq_lcl_dt, dn_lcl_dt, dq_rai_dt, dn_rai_dt,
-        dq_ice_dt, dn_ice_dt, dq_rim_dt, db_rim_dt)
-end
-
-"""
-    _limited_2mp3_tendency(mp, tps, ρ, T, q_tot,
-        q_lcl, n_lcl, q_rai, n_rai, q_ice, n_ice, q_rim, b_rim, logλ, dt)
-
-The instantaneous 2M+P3 tendency projected onto the eight prognostic
-species, with [`_cap_supersaturation`](@ref) and the coupled
-[`_coupled_sink_factor`](@ref) mass-number limits applied for a substep of
-length `dt`.
-
-This is the function the Rosenbrock step differentiates: the limiting sits
-inside the Jacobian so the implicit update respects it. Interrogated
-empirically: at substep lengths up to tens of seconds the supersaturation
-cap improves accuracy substantially (the coupled-sink limits are neutral
-there and protect the explicit fallback paths); at very large substeps
-(h ≳ 100 s) the limiter switching degrades the single linearization —
-increase `nsub` rather than relying on the limiters.
-"""
-@inline function _limited_2mp3_tendency(mp, tps,
+    _instantaneous_2mp3_tendency(mp, tps, ρ, T, q_tot,
+        q_lcl, n_lcl, q_rai, n_rai, q_ice, n_ice, q_rim, b_rim, logλ)
+
+The raw instantaneous 2M+P3 tendency projected onto the eight prognostic
+species: the unlimited process rates of the `Microphysics2Moment` entry,
+without timestep-dependent clipping.
+
+This is the function the Rosenbrock step differentiates. It is the model
+physics — an analytic time-averaged relaxation (Morrison and Milbrandt
+2015, J. Atmos. Sci., 72, 287-311, Appendix C) carries no tendency clip;
+the L-stable one-stage Rosenbrock update damps the stiff vapor-exchange
+subsystem monotonically, so the overshoot an explicit-Euler supersaturation
+cap suppresses cannot arise, while a 1/dt clip would add dt-independent
+error and break convergence under refinement (Wan et al. 2020, J. Adv.
+Model. Earth Syst., 12, e2019MS001982).
+"""
+@inline function _instantaneous_2mp3_tendency(mp, tps,
     ρ, T, q_tot,
     q_lcl, n_lcl, q_rai, n_rai, q_ice, n_ice, q_rim, b_rim, logλ,
-    dt,
 )
     full = bulk_microphysics_tendencies(Microphysics2Moment(), mp, tps,
         ρ, T, q_tot,
         q_lcl, n_lcl, q_rai, n_rai, q_ice, n_ice, q_rim, b_rim, logλ,
     )
-    tend = (; full.dq_lcl_dt, full.dn_lcl_dt, full.dq_rai_dt, full.dn_rai_dt,
+    return (; full.dq_lcl_dt, full.dn_lcl_dt, full.dq_rai_dt, full.dn_rai_dt,
         full.dq_ice_dt, full.dn_ice_dt, full.dq_rim_dt, full.db_rim_dt)
-    tend = _cap_supersaturation(tend, tps, T, ρ, q_tot, q_lcl, q_rai, q_ice, dt)
-    f_liq = _coupled_sink_factor(tend.dq_lcl_dt, tend.dn_lcl_dt, q_lcl, n_lcl, dt)
-    f_rai = _coupled_sink_factor(tend.dq_rai_dt, tend.dn_rai_dt, q_rai, n_rai, dt)
-    return (;
-        dq_lcl_dt = tend.dq_lcl_dt * f_liq, dn_lcl_dt = tend.dn_lcl_dt * f_liq,
-        dq_rai_dt = tend.dq_rai_dt * f_rai, dn_rai_dt = tend.dn_rai_dt * f_rai,
-        tend.dq_ice_dt, tend.dn_ice_dt, tend.dq_rim_dt, tend.db_rim_dt,
-    )
 end
 
 """
-    Limited2MP3Tendency(mp, tps, ρ, T, q_tot, logλ, dt)
+    Instantaneous2MP3Tendency(mp, tps, ρ, T, q_tot, logλ)
 
 Callable bundling the frozen per-substep context; applying it to the species
-vector evaluates [`_limited_2mp3_tendency`](@ref). A top-level struct rather
-than a closure, so `ForwardDiff` differentiates a concretely-typed callable.
+vector evaluates [`_instantaneous_2mp3_tendency`](@ref). A top-level struct
+rather than a closure, so `ForwardDiff` differentiates a concretely-typed
+callable.
 
 The frozen `q_tot` is promoted to the state's element type at the call: a
 zero-partial Dual is exact for a constant, and it keeps
@@ -166,21 +67,19 @@ slots never carry derivatives (the forward `gamma_inc` has no shape rule);
 `T` and `ρ` stay plain because they feed working-type computations that must
 remain floats.
 """
-struct Limited2MP3Tendency{P, H, F}
+struct Instantaneous2MP3Tendency{P, H, F}
     mp::P
     tps::H
     ρ::F
     T::F
     q_tot::F
     logλ::F
-    dt::F
 end
-@inline function (g::Limited2MP3Tendency)(x::SA.StaticVector{8})
+@inline function (g::Instantaneous2MP3Tendency)(x::SA.StaticVector{8})
     (q_lcl, n_lcl, q_rai, n_rai, q_ice, n_ice, q_rim, b_rim) = x
-    tend = _limited_2mp3_tendency(g.mp, g.tps,
+    tend = _instantaneous_2mp3_tendency(g.mp, g.tps,
         g.ρ, g.T, eltype(x)(g.q_tot),
         q_lcl, n_lcl, q_rai, n_rai, q_ice, n_ice, q_rim, b_rim, g.logλ,
-        g.dt,
     )
     return MicroState2MP3(values(tend)...)
 end
@@ -196,10 +95,9 @@ This is a scheme partition, not a preconditioner: with `P J P` in place of
 `J`, a masked channel's row of `(I/h - P J P)` is the identity row, so the
 solve returns exactly the forward-Euler update for that channel while
 healthy channels stay implicit (an IMEX-style splitting at channel
-granularity). The gate is load-bearing: through the coupled-sink limiter,
-near-empty channels produce finite but enormous Jacobian entries whose
-linearized steady state fabricates phantom number concentrations that
-substep refinement cannot remove.
+granularity). The gate is load-bearing: near-empty channels produce finite
+but enormous Jacobian entries whose linearized steady state fabricates
+phantom number concentrations that substep refinement cannot remove.
 """
 @inline function _rosenbrock_channel_mask(x::MicroState2MP3{FT}) where {FT}
     ϵ_empty = FT(1e-10)
@@ -253,13 +151,14 @@ end
         Δt, nsub = 1)
 
 Compute average 2M+P3 microphysics tendencies over `Δt` using `nsub`
-linearized-implicit (Rosenbrock-Euler) substeps.
+linearized-implicit (Rosenbrock-Euler) substeps of the raw instantaneous
+tendency.
 
 # Algorithm
 
 For each substep of `h = Δt / nsub`:
 
-1. Evaluate the limited tendency `f` ([`_limited_2mp3_tendency`](@ref)) and
+1. Evaluate the raw tendency `f` ([`_instantaneous_2mp3_tendency`](@ref)) and
    its exact 8×8 Jacobian `J` via `ForwardDiff` at the current state.
 2. Advance with [`_rosenbrock_update`](@ref): solve `(I/h - P J P) Δx = f`
    in equilibrated variables, where the projection `P` routes near-empty
@@ -269,11 +168,11 @@ For each substep of `h = Δt / nsub`:
 
 A non-finite state cannot be differentiated and a non-finite Jacobian (an
 exotic state escaping the channel mask) cannot be linearized; both fall back
-to a forward-Euler substep of the limited tendency, which is safe to step
-explicitly. `logλ` and `q_tot` are held fixed across substeps, matching the
-explicit-substepping semantics; the limiting (supersaturation cap and
-coupled mass-number sinks) applies the explicit-timestepping semantics
-unconditionally and sits inside the differentiated function.
+to a forward-Euler substep of the raw tendency. `logλ` and `q_tot` are held
+fixed across substeps, matching the explicit-substepping semantics. The
+discrete safeguards are h-free conditioning and projection (channel mask,
+equilibration, positivity clamp), not model terms: the differentiated
+tendency is the unmodified relaxation physics.
 
 Returns the net change in the species over `Δt` divided by `Δt`, in the same
 fields as the `Instantaneous` entry (without the activation diagnostic).
@@ -293,7 +192,7 @@ fields as the `Instantaneous` entry (without the activation diagnostic).
     x₀ = x
     Tsub = T
     for _ in 1:nsub_eff
-        g = Limited2MP3Tendency(mp, tps, ρ, Tsub, q_tot, logλ, h)
+        g = Instantaneous2MP3Tendency(mp, tps, ρ, Tsub, q_tot, logλ)
         f = g(x)
         x_prev = x
         x = if all(isfinite, x)

src/BulkMicrophysicsTendencies.jl

@@ -122,7 +122,7 @@ struct LinearizedAverage <: TendencyMode end
 
 Return time-averaged tendencies computed via repeated linearized-implicit
 (Rosenbrock-Euler) substeps, linearizing with the exact `ForwardDiff` Jacobian
-of the limited pointwise tendency at each substep. Implemented for the
+of the raw instantaneous pointwise tendency at each substep. Implemented for the
 2-moment + P3 configuration.
 """
 struct RosenbrockAverage <: TendencyMode end