Add average bulk microphyiscs 1m tendency based on linearized microphysics: donor-based linearization (sinks as D*q, vapor sources constant); add operator construction and average tendency solve.

Author: sajjadazimi

Project.toml

@@ -1,6 +1,6 @@
 name = "CloudMicrophysics"
 uuid = "6a9e3e04-43cd-43ba-94b9-e8782df3c71b"
-version = "0.34"
+version = "0.35"
 authors = ["Climate Modeling Alliance"]
 
 [deps]

src/BulkMicrophysicsTendencies.jl

@@ -279,6 +279,279 @@ This is a pure function of local thermodynamic state, suitable for:
     return (; dq_lcl_dt, dq_icl_dt, dq_rai_dt, dq_sno_dt)
 end
 
+"""
+Construct a local linear approximation of 1-moment microphysics tendencies:
+
+    dq/dt ≈ M * q + e
+
+using a donor-based linearization:
+
+- Donor → receiver transfers are linearized as `D * q_donor`
+- Vapor → condensate sources are treated as constants (`e`)
+- Condensate → vapor sinks are treated as linear sinks (`-D * q`)
+
+All coefficients use `D = S / max(eps(T), q_donor)` for robustness.
+
+With this formulation, sink terms take the form
+
+    dq/dt = -D q
+
+which corresponds to exponential decay over a timestep in the implicit solve.
+
+The resulting operator is sparse with nonzeros only in:
+- Diagonal: `M11`, `M22`, `M33`, `M44`
+- Couplings: `M31`, `M34`, `M41`, `M42`, `M43`
+- Sources: `e1`, `e2`, `e4` (rain has no constant source)
+
+Returns a `NamedTuple` with these entries.
+"""
+@inline function bulk_microphysics_linearized_operator(
+    ::Microphysics1Moment, mp::CMP.Microphysics1MParams, tps,
+    ρ, T, q_tot, q_lcl, q_icl, q_rai, q_sno, N_lcl = zero(ρ),
+)
+    ρ = UT.clamp_to_nonneg(ρ)
+    q_tot = UT.clamp_to_nonneg(q_tot)
+    q_lcl = UT.clamp_to_nonneg(q_lcl)
+    q_icl = UT.clamp_to_nonneg(q_icl)
+    q_rai = UT.clamp_to_nonneg(q_rai)
+    q_sno = UT.clamp_to_nonneg(q_sno)
+
+    lcl = mp.cloud.liquid
+    icl = mp.cloud.ice
+    rai = mp.precip.rain
+    sno = mp.precip.snow
+    ce = mp.collision
+    aps = mp.air_properties
+    vel = mp.terminal_velocity
+    var = mp.autoconv_2M
+
+    FT = promote_type(
+        typeof(ρ), typeof(T), typeof(q_tot), typeof(q_lcl),
+        typeof(q_icl), typeof(q_rai), typeof(q_sno),
+    )
+
+    M11 = zero(FT)
+    M22 = zero(FT)
+    M31 = zero(FT)
+    M33 = zero(FT)
+    M34 = zero(FT)
+    M41 = zero(FT)
+    M42 = zero(FT)
+    M43 = zero(FT)
+    M44 = zero(FT)
+    e1 = zero(FT)
+    e2 = zero(FT)
+    e4 = zero(FT)
+
+    # ------------------------------------------------------------
+    # Vapor -> cloud condensate: constant sources
+    # ------------------------------------------------------------
+    S_lcl_cond = CMNonEq.conv_q_vap_to_q_lcl_icl_MM2015(
+        lcl, tps, q_tot, q_lcl, q_icl, q_rai, q_sno, ρ, T,
+    )
+    if S_lcl_cond >= zero(FT)
+        e1 += S_lcl_cond
+    else
+        D = (-S_lcl_cond) / max(eps(FT), q_lcl)
+        M11 -= D
+    end
+
+    S_icl_dep = CMNonEq.conv_q_vap_to_q_lcl_icl_MM2015(
+        icl, tps, q_tot, q_lcl, q_icl, q_rai, q_sno, ρ, T,
+    )
+    if S_icl_dep >= zero(FT)
+        e2 += S_icl_dep
+    else
+        D = (-S_icl_dep) / max(eps(FT), q_icl)
+        M22 -= D
+    end
+
+    # ------------------------------------------------------------
+    # Autoconversion: donor-based transfer
+    # ------------------------------------------------------------
+    S_acnv_1M = CM1.conv_q_lcl_to_q_rai(rai.acnv1M, q_lcl, true)
+    S_acnv_2M = isnothing(var) ? S_acnv_1M : CM2.conv_q_lcl_to_q_rai(var, q_lcl, ρ, N_lcl)
+    S_acnv_lcl = ifelse(N_lcl > zero(N_lcl), S_acnv_2M, S_acnv_1M)
+    D = S_acnv_lcl / max(eps(FT), q_lcl)
+    M11 -= D
+    M31 += D
+
+    S_acnv_icl = CM1.conv_q_icl_to_q_sno_no_supersat(sno.acnv1M, q_icl, true)
+    D = S_acnv_icl / max(eps(FT), q_icl)
+    M22 -= D
+    M42 += D
+
+    # ------------------------------------------------------------
+    # Accretion: donor-based transfer
+    # ------------------------------------------------------------
+    S_accr_lcl_rai = CM1.accretion(lcl, rai, vel.rain, ce, q_lcl, q_rai, ρ)
+    D = S_accr_lcl_rai / max(eps(FT), q_lcl)
+    M11 -= D
+    M31 += D
+
+    S_accr_lcl_sno = CM1.accretion(lcl, sno, vel.snow, ce, q_lcl, q_sno, ρ)
+    D = S_accr_lcl_sno / max(eps(FT), q_lcl)
+    M11 -= D
+    is_warm = T >= sno.T_freeze
+    if is_warm
+        M31 += D
+    else
+        M41 += D
+    end
+
+    α = warm_accretion_melt_factor(tps, sno, T)
+    S_accr_melt = α * S_accr_lcl_sno
+    # this is snow -> rain melt induced by warm collected liquid
+    # donor is snow
+    D = S_accr_melt / max(eps(FT), q_sno)
+    M44 -= D
+    M34 += D
+
+    S_accr_icl_rai = CM1.accretion(icl, rai, vel.rain, ce, q_icl, q_rai, ρ)
+    D = S_accr_icl_rai / max(eps(FT), q_icl)
+    M22 -= D
+    M42 += D
+
+    S_accr_icl_sno = CM1.accretion(icl, sno, vel.snow, ce, q_icl, q_sno, ρ)
+    D = S_accr_icl_sno / max(eps(FT), q_icl)
+    M22 -= D
+    M42 += D
+
+    S_accr_rai_icl = CM1.accretion_rain_sink(rai, icl, vel.rain, ce, q_icl, q_rai, ρ)
+    D = S_accr_rai_icl / max(eps(FT), q_rai)
+    M33 -= D
+    M43 += D
+
+    S_accr_rai_sno = CM1.accretion_snow_rain(sno, rai, vel.snow, vel.rain, ce, q_sno, q_rai, ρ)
+    S_accr_sno_rai = CM1.accretion_snow_rain(rai, sno, vel.rain, vel.snow, ce, q_rai, q_sno, ρ)
+
+    if is_warm
+        # snow -> rain
+        D = S_accr_sno_rai / max(eps(FT), q_sno)
+        M44 -= D
+        M34 += D
+
+        S_accr_melt2 = α * S_accr_rai_sno
+        D = S_accr_melt2 / max(eps(FT), q_sno)
+        M44 -= D
+        M34 += D
+    else
+        # rain -> snow
+        D = S_accr_rai_sno / max(eps(FT), q_rai)
+        M33 -= D
+        M43 += D
+    end
+
+    # ------------------------------------------------------------
+    # Rain evaporation: sink to vapor if negative, else constant source
+    # ------------------------------------------------------------
+    S_evap_rai = CM1.evaporation_sublimation(
+        rai, vel.rain, aps, tps, q_tot, q_lcl, q_icl, q_rai, q_sno, ρ, T,
+    )
+    D = (-S_evap_rai) / max(eps(FT), q_rai)
+    M33 -= D
+
+    # ------------------------------------------------------------
+    # Snow sublimation/deposition
+    # ------------------------------------------------------------
+    S_subl_sno = CM1.evaporation_sublimation(
+        sno, vel.snow, aps, tps, q_tot, q_lcl, q_icl, q_rai, q_sno, ρ, T,
+    )
+    if S_subl_sno < zero(FT)
+        D = (-S_subl_sno) / max(eps(FT), q_sno)
+        M44 -= D
+    else
+        e4 += S_subl_sno
+    end
+
+    # ------------------------------------------------------------
+    # Snow melt: snow -> rain
+    # ------------------------------------------------------------
+    S_melt_sno = CM1.snow_melt(sno, vel.snow, aps, tps, q_sno, ρ, T)
+    D = S_melt_sno / max(eps(FT), q_sno)
+    M44 -= D
+    M34 += D
+
+    return (
+        M11 = M11, M22 = M22,
+        M31 = M31, M33 = M33, M34 = M34,
+        M41 = M41, M42 = M42, M43 = M43, M44 = M44,
+        e1 = e1, e2 = e2, e4 = e4,
+    )
+end
+
+"""
+Compute time-averaged microphysics tendencies over a period of time Δt.
+
+Solves the linearized implicit system:
+
+    (q* - q⁰) / Δt = M q* + e
+
+and returns the average tendency
+
+    dq/dt = (q* - q⁰) / Δt
+
+The system is solved using a specialized sparse structure (specific to the microphysics 1M model):
+- `q_lcl`, `q_icl`: independent scalar solves
+- `(q_rai, q_sno)`: 2×2 coupled solve
+
+Because sinks are linearized as `-D q`, this method effectively integrates
+them as exponential decays over the timestep, improving stability for stiff processes.
+
+Returns a `NamedTuple` of averaged tendencies for all species.
+"""
+@inline function average_bulk_microphysics_tendencies(
+    ::Microphysics1Moment, mp::CMP.Microphysics1MParams, tps,
+    ρ, T, q_tot, q_lcl, q_icl, q_rai, q_sno, Δt;
+    N_lcl = zero(ρ),
+)
+
+    FT = typeof(q_tot)
+
+    lin = bulk_microphysics_linearized_operator(
+        Microphysics1Moment(), mp, tps,
+        ρ, T, q_tot, q_lcl, q_icl, q_rai, q_sno, N_lcl,
+    )
+
+    invΔt = one(FT) / Δt
+
+    # A = I/Δt - M
+    a11 = invΔt - lin.M11
+    a22 = invΔt - lin.M22
+    a31 = -lin.M31
+    a33 = invΔt - lin.M33
+    a34 = -lin.M34
+    a41 = -lin.M41
+    a42 = -lin.M42
+    a43 = -lin.M43
+    a44 = invΔt - lin.M44
+
+    # rhs = e + q_0/Δt
+    # e3 = 0 by the 1m model
+    b1 = lin.e1 + invΔt * q_lcl
+    b2 = lin.e2 + invΔt * q_icl
+    b3 = invΔt * q_rai
+    b4 = lin.e4 + invΔt * q_sno
+
+    q_lcl_new = b1 / a11
+    q_icl_new = b2 / a22
+
+    # Reduced 2x2 system for q_rai_new, q_sno_new
+    r3 = muladd(-a31, q_lcl_new, b3)
+    r4 = muladd(-a41, q_lcl_new, muladd(-a42, q_icl_new, b4))
+
+    det = muladd(-a34, a43, a33 * a44)
+    q_rai_new = (r3 * a44 - a34 * r4) / det
+    q_sno_new = (a33 * r4 - r3 * a43) / det
+
+    dq_lcl_dt = (q_lcl_new - q_lcl) * invΔt
+    dq_icl_dt = (q_icl_new - q_icl) * invΔt
+    dq_rai_dt = (q_rai_new - q_rai) * invΔt
+    dq_sno_dt = (q_sno_new - q_sno) * invΔt
+
+    return (; dq_lcl_dt, dq_icl_dt, dq_rai_dt, dq_sno_dt)
+end
+
 """
     bulk_microphysics_derivatives(
         ::Microphysics1Moment,