Added sink derivatives (#691)

* Added sink derivatives

* Replace ternary by GPU performant safe division in derivatives

Author: tapios

Project.toml

@@ -1,7 +1,7 @@
 name = "CloudMicrophysics"
 uuid = "6a9e3e04-43cd-43ba-94b9-e8782df3c71b"
 authors = ["Climate Modeling Alliance"]
-version = "0.31.10"
+version = "0.31.11"
 
 [deps]
 ClimaParams = "5c42b081-d73a-476f-9059-fd94b934656c"

src/BulkMicrophysicsTendencies.jl

@@ -234,6 +234,7 @@ This is a pure function of local thermodynamic state, suitable for:
     dq_sno_dt += ifelse(is_warm, zero(T), S_accr_lcl_sno)
     dq_rai_dt += ifelse(is_warm, S_accr_lcl_sno, zero(T))
     # Thermal melting (α=0 when cold, so these terms vanish)
+    # Physical consistency: melt is mass of snow melted, which is α * mass of warm liquid
     α = warm_accretion_melt_factor(tps, sno, T)
     S_accr_melt = α * S_accr_lcl_sno
     dq_sno_dt -= S_accr_melt
@@ -269,7 +270,8 @@ This is a pure function of local thermodynamic state, suitable for:
     dq_sno_dt -= ifelse(is_warm, S_accr_sno_rai, zero(T))
     dq_rai_dt += ifelse(is_warm, S_accr_sno_rai, zero(T))
     # Thermal melting from warm rain (α=0 when cold)
-    S_accr_melt = α * S_accr_sno_rai
+    # Physical consistency: melt is mass of snow melted, which is α * mass of warm rain
+    S_accr_melt = α * S_accr_rai_sno
     dq_sno_dt -= ifelse(is_warm, S_accr_melt, zero(T))
     dq_rai_dt += ifelse(is_warm, S_accr_melt, zero(T))
 
@@ -350,21 +352,61 @@ as `bulk_microphysics_tendencies`.
     aps = mp.air_properties
     vel = mp.terminal_velocity
 
-    # Condensation/evaporation of cloud liquid
+    # --- Cloud liquid: condensation + sink self-derivatives ---
     ∂tendency_∂q_lcl = CMNonEq.∂conv_q_vap_to_q_lcl_icl_MM2015_∂q_cld(lcl, tps, q_tot, q_lcl, q_icl, q_rai, q_sno, ρ, T)
+    # Autoconversion q_lcl → q_rai (sink, derivative ∂/∂q_lcl is conservative and safe for linear-above-threshold)
+    S_acnv_lcl = CM1.conv_q_lcl_to_q_rai(rai.acnv1M, q_lcl, true)
+    # Accretion q_lcl + q_rai → q_rai (rate is exactly linear in q_lcl)
+    S_accr_lcl_rai = CM1.accretion(lcl, rai, vel.rain, ce, q_lcl, q_rai, ρ)
+    # Accretion q_lcl + q_sno → q_sno (rate is exactly linear in q_lcl)
+    S_accr_lcl_sno = CM1.accretion(lcl, sno, vel.snow, ce, q_lcl, q_sno, ρ)
+    ∂tendency_∂q_lcl -=
+        (S_acnv_lcl + S_accr_lcl_rai + S_accr_lcl_sno) /
+        max(q_lcl, UT.ϵ_numerics(typeof(q_lcl)))
 
-    # Deposition/sublimation of cloud ice (INP limiter applied inside ∂conv_q_vap_to_q_lcl_icl_MM2015_∂q_cld)
+    # --- Cloud ice: deposition + sink self-derivatives ---
     ∂tendency_∂q_icl = CMNonEq.∂conv_q_vap_to_q_lcl_icl_MM2015_∂q_cld(icl, tps, q_tot, q_lcl, q_icl, q_rai, q_sno, ρ, T)
+    # Autoconversion q_icl → q_sno
+    S_acnv_icl = CM1.conv_q_icl_to_q_sno_no_supersat(sno.acnv1M, q_icl, true)
+    # Accretion q_icl + q_rai → q_sno (rate is exactly linear in q_icl)
+    S_accr_icl_rai = CM1.accretion(icl, rai, vel.rain, ce, q_icl, q_rai, ρ)
+    # Accretion q_icl + q_sno → q_sno (rate is exactly linear in q_icl)
+    S_accr_icl_sno = CM1.accretion(icl, sno, vel.snow, ce, q_icl, q_sno, ρ)
+    ∂tendency_∂q_icl -=
+        (S_acnv_icl + S_accr_icl_rai + S_accr_icl_sno) /
+        max(q_icl, UT.ϵ_numerics(typeof(q_icl)))
 
-    # Rain evaporation
+    # --- Rain: evaporation + sink self-derivatives ---
     ∂tendency_∂q_rai =
         CM1.∂evaporation_sublimation_∂q_precip(rai, vel.rain, aps, tps, q_tot, q_lcl, q_icl, q_rai, q_sno, ρ, T)
+    # Accretion rain_sink: ice + rain → snow (rate/q_rai approximation)
+    S_accr_rai_icl = CM1.accretion_rain_sink(rai, icl, vel.rain, ce, q_icl, q_rai, ρ)
+    # Snow-rain accretion (cold): rain → snow
+    is_warm = T >= sno.T_freeze
+    S_accr_rai_sno = ifelse(is_warm, zero(T),
+        CM1.accretion_snow_rain(sno, rai, vel.snow, vel.rain, ce, q_sno, q_rai, ρ))
+    ∂tendency_∂q_rai -=
+        (S_accr_rai_icl + S_accr_rai_sno) /
+        max(q_rai, UT.ϵ_numerics(typeof(q_rai)))
 
-    # Snow sublimation/deposition and melting
+    # --- Snow: sublimation/deposition + melting + sink self-derivatives ---
     dS_subl_sno =
         CM1.∂evaporation_sublimation_∂q_precip(sno, vel.snow, aps, tps, q_tot, q_lcl, q_icl, q_rai, q_sno, ρ, T)
     dS_melt_sno = CM1.∂snow_melt_∂q_sno(sno, vel.snow, aps, tps, q_sno, ρ, T)
     ∂tendency_∂q_sno = dS_subl_sno - dS_melt_sno
+    # Snow-rain accretion (warm): snow → rain
+    S_accr_sno_rai = ifelse(is_warm,
+        CM1.accretion_snow_rain(rai, sno, vel.rain, vel.snow, ce, q_rai, q_sno, ρ), zero(T))
+    # Thermal melting from warm rain collision -> snow sink proportional to S_accr_rai_sno
+    # (Since this is a diagonal derivative, we approximate ∂/∂q_sno by rate/q_sno, using the rain sink form S_accr_rai_sno)
+    S_accr_sno_rai_melt = ifelse(is_warm,
+        CM1.accretion_snow_rain(sno, rai, vel.snow, vel.rain, ce, q_sno, q_rai, ρ), zero(T))
+    # Accretion melt from warm liquid-snow collision
+    α = warm_accretion_melt_factor(tps, sno, T)
+    S_accr_lcl_sno_melt = α * CM1.accretion(lcl, sno, vel.snow, ce, q_lcl, q_sno, ρ)
+    ∂tendency_∂q_sno -=
+        (S_accr_sno_rai + α * S_accr_sno_rai_melt + S_accr_lcl_sno_melt) /
+        max(q_sno, UT.ϵ_numerics(typeof(q_sno)))
 
     return (; ∂tendency_∂q_lcl, ∂tendency_∂q_icl, ∂tendency_∂q_rai, ∂tendency_∂q_sno)
 end

src/Microphysics1M.jl

@@ -722,7 +722,7 @@ The approximation used is `∂(rate)/∂q_precip ≈ rate / q_precip`.
     T::FT,
 ) where {FT}
     rate = evaporation_sublimation(rain_params, vel, aps, tps, q_tot, q_lcl, q_icl, q_rai, q_sno, ρ, T)
-    return q_rai > UT.ϵ_numerics(FT) ? rate / q_rai : zero(rate)
+    return rate / max(q_rai, UT.ϵ_numerics(FT))
 end
 
 @inline function ∂evaporation_sublimation_∂q_precip(
@@ -739,7 +739,7 @@ end
     T::FT,
 ) where {FT}
     rate = evaporation_sublimation(snow_params, vel, aps, tps, q_tot, q_lcl, q_icl, q_rai, q_sno, ρ, T)
-    return q_sno > UT.ϵ_numerics(FT) ? rate / q_sno : zero(rate)
+    return rate / max(q_sno, UT.ϵ_numerics(FT))
 end
 
 """
@@ -815,7 +815,7 @@ The approximation used is `∂(rate)/∂q_sno ≈ rate / q_sno`.
     T::FT,
 ) where {FT}
     rate = snow_melt(snow_params, vel, aps, tps, q_sno, ρ, T)
-    return q_sno > UT.ϵ_numerics(FT) ? rate / q_sno : zero(rate)
+    return rate / max(q_sno, UT.ϵ_numerics(FT))
 end
 
 end #module Microphysics1M.jl

test/bulk_tendencies_tests.jl

@@ -710,6 +710,135 @@ function test_bulk_microphysics_1m_derivatives(FT)
         )
         @test all(isfinite, derivs)
     end
+
+    @testset "BulkMicrophysicsDerivatives 1M - Cloud liquid sink derivatives" begin
+        # With cloud liquid above autoconversion threshold and rain/snow present,
+        # the derivative should be negative (sinks from autoconversion + accretion)
+        ρ = FT(1.2)
+        T = T_freeze + FT(5)
+        q_lcl = FT(2e-3)  # above autoconversion threshold
+        q_icl = FT(0)
+        q_rai = FT(5e-4)  # rain present for accretion
+        q_sno = FT(5e-4)  # snow present for accretion
+        q_tot = get_saturated_q_tot(tps, T, ρ, q_lcl, q_icl, q_rai, q_sno)
+
+        derivs = BMT.bulk_microphysics_derivatives(
+            BMT.Microphysics1Moment(),
+            mp, tps, ρ, T, q_tot, q_lcl, q_icl, q_rai, q_sno,
+        )
+        # Derivative should be negative (sink terms dominate at saturation)
+        @test derivs.∂tendency_∂q_lcl < FT(0)
+
+        # Compare with no-sink case (no rain, no snow, below autoconv threshold)
+        q_lcl_small = FT(1e-5)  # below threshold
+        q_tot_small = get_saturated_q_tot(tps, T, ρ, q_lcl_small, q_icl, FT(0), FT(0))
+        derivs_nosink = BMT.bulk_microphysics_derivatives(
+            BMT.Microphysics1Moment(),
+            mp, tps, ρ, T, q_tot_small, q_lcl_small, q_icl, FT(0), FT(0),
+        )
+        # With sinks, derivative should be more negative
+        @test derivs.∂tendency_∂q_lcl < derivs_nosink.∂tendency_∂q_lcl
+    end
+
+    @testset "BulkMicrophysicsDerivatives 1M - Cloud ice sink derivatives" begin
+        # With cloud ice above autoconversion threshold and snow present,
+        # the derivative should include sink contributions
+        ρ = FT(1.2)
+        T = T_freeze - FT(20)  # cold
+        q_lcl = FT(0)
+        q_icl = FT(1e-3)  # above autoconversion threshold
+        q_rai = FT(0)
+        q_sno = FT(5e-4)  # snow present for accretion
+        q_tot = FT(0.012)
+
+        derivs = BMT.bulk_microphysics_derivatives(
+            BMT.Microphysics1Moment(),
+            mp, tps, ρ, T, q_tot, q_lcl, q_icl, q_rai, q_sno,
+        )
+        # Derivative should be negative (autoconversion + accretion sinks)
+        @test derivs.∂tendency_∂q_icl < FT(0)
+    end
+
+    @testset "BulkMicrophysicsDerivatives 1M - Rain sink derivatives (cold)" begin
+        # In cold conditions with both rain and snow, rain is lost to snow
+        ρ = FT(1.2)
+        T = T_freeze - FT(5)  # below freezing
+        q_lcl = FT(0)
+        q_icl = FT(1e-4)  # ice present for accretion_rain_sink
+        q_rai = FT(5e-4)
+        q_sno = FT(5e-4)
+        q_tot = FT(0.012)
+
+        derivs = BMT.bulk_microphysics_derivatives(
+            BMT.Microphysics1Moment(),
+            mp, tps, ρ, T, q_tot, q_lcl, q_icl, q_rai, q_sno,
+        )
+        # Rain derivative should be negative (evaporation not active in subsaturated,
+        # plus accretion sinks from ice and snow collisions)
+        @test derivs.∂tendency_∂q_rai < FT(0)
+    end
+
+    @testset "BulkMicrophysicsDerivatives 1M - Snow sink derivatives (warm)" begin
+        # In warm conditions with snow, melting and accretion sinks should make
+        # the derivative strongly negative
+        ρ = FT(1.2)
+        T = T_freeze + FT(5)  # above freezing
+        q_lcl = FT(5e-4)  # liquid present for accretion melt
+        q_icl = FT(0)
+        q_rai = FT(5e-4)  # rain present for snow-rain accretion
+        q_sno = FT(1e-3)
+        q_tot = FT(0.015)
+
+        derivs = BMT.bulk_microphysics_derivatives(
+            BMT.Microphysics1Moment(),
+            mp, tps, ρ, T, q_tot, q_lcl, q_icl, q_rai, q_sno,
+        )
+        # Snow derivative should be strongly negative (melting + accretion sinks)
+        @test derivs.∂tendency_∂q_sno < FT(0)
+
+        # Compare: snow derivative should be more negative in warm than cold
+        T_cold = T_freeze - FT(10)
+        derivs_cold = BMT.bulk_microphysics_derivatives(
+            BMT.Microphysics1Moment(),
+            mp, tps, ρ, T_cold, q_tot, q_lcl, q_icl, q_rai, q_sno,
+        )
+        @test derivs.∂tendency_∂q_sno < derivs_cold.∂tendency_∂q_sno
+    end
+
+    @testset "BulkMicrophysicsDerivatives 1M - Finite difference validation" begin
+        # Validate derivatives against finite differences for a representative state
+        ρ = FT(1.2)
+        T = T_freeze + FT(3)  # warm, snow processes active
+        q_lcl = FT(5e-4)
+        q_icl = FT(0)
+        q_rai = FT(5e-4)
+        q_sno = FT(5e-4)
+        q_tot = FT(0.015)
+
+        derivs = BMT.bulk_microphysics_derivatives(
+            BMT.Microphysics1Moment(),
+            mp, tps, ρ, T, q_tot, q_lcl, q_icl, q_rai, q_sno,
+        )
+
+        # Finite difference for ∂(dq_sno_dt)/∂q_sno
+        δ = FT(1e-8)
+        tend_plus = BMT.bulk_microphysics_tendencies(
+            BMT.Microphysics1Moment(),
+            mp, tps, ρ, T, q_tot, q_lcl, q_icl, q_rai, q_sno + δ,
+        )
+        tend_minus = BMT.bulk_microphysics_tendencies(
+            BMT.Microphysics1Moment(),
+            mp, tps, ρ, T, q_tot, q_lcl, q_icl, q_rai, q_sno - δ,
+        )
+        fd_sno = (tend_plus.dq_sno_dt - tend_minus.dq_sno_dt) / (2δ)
+
+        # The analytic derivative should be within a factor of 3 of the finite difference
+        # (the rate/q approximation is not exact, so we use a loose tolerance)
+        @test derivs.∂tendency_∂q_sno < FT(0)
+        @test fd_sno < FT(0)
+        ratio = derivs.∂tendency_∂q_sno / fd_sno
+        @test FT(0.1) < ratio < FT(10)
+    end
 end
 
 ###