From e6581c31d84ee97fd47f87c94d1a2f14aaf07616 Mon Sep 17 00:00:00 2001 From: Anna Jaruga Date: Wed, 25 Mar 2026 11:05:42 -0700 Subject: [PATCH] Use condensate in eaporation/sublimation --- Project.toml | 2 +- docs/src/API.md | 3 - docs/src/MicrophysicsNonEq.md | 30 ----- docs/src/plots/NonEqCondEvapRate.jl | 195 ---------------------------- src/BulkMicrophysicsTendencies.jl | 6 +- src/MicrophysicsNonEq.jl | 122 +++-------------- test/gpu_tests.jl | 15 +-- test/microphysics_noneq_tests.jl | 104 +-------------- test/performance_tests.jl | 2 +- 9 files changed, 29 insertions(+), 450 deletions(-) diff --git a/Project.toml b/Project.toml index ac88928d8..75ec69253 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "CloudMicrophysics" uuid = "6a9e3e04-43cd-43ba-94b9-e8782df3c71b" -version = "0.32" +version = "0.33" authors = ["Climate Modeling Alliance"] [deps] diff --git a/docs/src/API.md b/docs/src/API.md index a36279ed5..17a424250 100644 --- a/docs/src/API.md +++ b/docs/src/API.md @@ -11,13 +11,10 @@ MicrophysicsNonEq MicrophysicsNonEq.τ_relax MicrophysicsNonEq.conv_q_vap_to_q_lcl_icl MicrophysicsNonEq.conv_q_vap_to_q_lcl_icl_MM2015 -MicrophysicsNonEq.∂conv_q_vap_to_q_lcl_icl_MM2015_∂q_cld -MicrophysicsNonEq.limit_MM2015_sinks MicrophysicsNonEq.INP_limiter MicrophysicsNonEq.terminal_velocity MicrophysicsNonEq.dqcld_dT MicrophysicsNonEq.gamma_helper -MicrophysicsNonEq.d2qcld_dT2 ``` # 0-moment precipitation microphysics diff --git a/docs/src/MicrophysicsNonEq.md b/docs/src/MicrophysicsNonEq.md index 524c78e5b..f627adc2a 100644 --- a/docs/src/MicrophysicsNonEq.md +++ b/docs/src/MicrophysicsNonEq.md @@ -124,38 +124,8 @@ include("plots/NonEqCondEvapRate.jl") ![](condensation_evaporation_ql_z.svg) ![](condensation_evaporation_ql_T.svg) -### Derivative of the tendency -For implicit time integration, we also compute the total derivative of the tendency - with respect to the cloud condensate, accounting for the implicit feedback - between the condensate and temperature through latent heat release - (``dT/dq = L/c_p``): -```math -\begin{equation} - \frac{d\dot{q}_{lcl}}{dq_{lcl}} = -\frac{1}{\tau_l} - - \frac{\dot{q}_{lcl}}{\Gamma_l} \left(\frac{L_v}{c_p}\right)^2 \frac{d^2 q_{sl}}{dT^2}; - \;\;\;\;\;\;\;\; - \frac{d\dot{q}_{icl}}{dq_{icl}} = -\frac{1}{\tau_i} - - \frac{\dot{q}_{icl}}{\Gamma_i} \left(\frac{L_s}{c_p}\right)^2 \frac{d^2 q_{si}}{dT^2} -\end{equation} -``` -where ``\dot{q}_{lcl}`` and ``\dot{q}_{icl}`` denote the condensation/evaporation - and deposition/sublimation tendencies respectively, - and the second derivative of saturation specific humidity is: -```math -\begin{equation} - \frac{d^2 q_{sl}}{dT^2} = q_{sl} \left[\left(\frac{L_v}{R_v T^2} - \frac{1}{T}\right)^2 - + \frac{1}{T^2} - \frac{2 L_v}{R_v T^3}\right] -\end{equation} -``` - -The leading-order term ``-1/\tau`` provides a simple approximation, - while the second term captures the higher-order correction from the - thermodynamic adjustment factor ``\Gamma``. - -![](condensation_evaporation_deriv.svg) -![](condensation_evaporation_deriv_vs_qvap.svg) ## Cloud condensate sedimentation diff --git a/docs/src/plots/NonEqCondEvapRate.jl b/docs/src/plots/NonEqCondEvapRate.jl index fded88c0a..308aa557a 100644 --- a/docs/src/plots/NonEqCondEvapRate.jl +++ b/docs/src/plots/NonEqCondEvapRate.jl @@ -98,7 +98,6 @@ function generate_cond_evap_rate(::Range_qₗ_T) ) end -# --- Derivative line plots --- minussign(s) = replace(s, "-" => "−") @@ -122,200 +121,6 @@ function make_figure(case) save(file_name, fig) end -function make_deriv_line_plot() - # Temperature range - T_range = collect(220:0.5:310) - - # Fixed parameters - qₜ = FT(10e-3) # 10 g/kg total water - qₗ = FT(1e-3) # 1 g/kg cloud liquid - qᵢ = FT(1e-3) # 1 g/kg cloud ice - qᵣ = FT(0) - qₛ = FT(0) - ρ = FT(1) # kg/m³ - - τ_relax = FT(1) # 1 s relaxation time - _liq = CM.Parameters.CloudLiquid(FT) - cm_liq = CM.Parameters.CloudLiquid(τ_relax, _liq.ρw, _liq.r_eff, _liq.N_0) - _ice = CM.Parameters.CloudIce(FT) - cm_ice = CM.Parameters.CloudIce(_ice.pdf, _ice.mass, _ice.r0, _ice.r_ice_snow, τ_relax, _ice.ρᵢ, _ice.r_eff, _ice.N_0) - - # Compute tendency and derivatives for each temperature - S_liq = zeros(FT, length(T_range)) - dS_liq_simple = zeros(FT, length(T_range)) - dS_liq_full = zeros(FT, length(T_range)) - - S_ice = zeros(FT, length(T_range)) - dS_ice_simple = zeros(FT, length(T_range)) - dS_ice_full = zeros(FT, length(T_range)) - - for (i, T) in enumerate(T_range) - T_ft = FT(T) - s = CMNe.conv_q_vap_to_q_lcl_icl_MM2015( - cm_liq, thp, qₜ, qₗ, qᵢ, qᵣ, qₛ, ρ, T_ft, - ) - ds_s = CMNe.∂conv_q_vap_to_q_lcl_icl_MM2015_∂q_cld( - cm_liq, thp, qₜ, qₗ, qᵢ, qᵣ, qₛ, ρ, T_ft, - ) - ds_f = CMNe.∂conv_q_vap_to_q_lcl_icl_MM2015_∂q_cld( - cm_liq, thp, qₜ, qₗ, qᵢ, qᵣ, qₛ, ρ, T_ft; - simplified = false, - ) - S_liq[i] = s - dS_liq_simple[i] = ds_s - dS_liq_full[i] = ds_f - - s = CMNe.conv_q_vap_to_q_lcl_icl_MM2015( - cm_ice, thp, qₜ, qₗ, qᵢ, qᵣ, qₛ, ρ, T_ft, - ) - ds_s = CMNe.∂conv_q_vap_to_q_lcl_icl_MM2015_∂q_cld( - cm_ice, thp, qₜ, qₗ, qᵢ, qᵣ, qₛ, ρ, T_ft, - ) - ds_f = CMNe.∂conv_q_vap_to_q_lcl_icl_MM2015_∂q_cld( - cm_ice, thp, qₜ, qₗ, qᵢ, qᵣ, qₛ, ρ, T_ft; - simplified = false, - ) - S_ice[i] = s - dS_ice_simple[i] = ds_s - dS_ice_full[i] = ds_f - end - - T_freeze = FT(TD.Parameters.T_freeze(thp)) - - # --- Figure: 2 panels --- - fig = Figure(size = (900, 450)) - - # Left panel: tendency vs T - ax1 = Axis(fig[1, 1]; - xlabel = "Temperature (K)", - ylabel = "Tendency (1/s)", - title = "Condensation/Evaporation & Deposition/Sublimation\n(qₜ=10, qₗ=1, qᵢ=1 g/kg, ρ=1 kg/m³, τ=1 s)", - ) - lines!(ax1, T_range, S_liq; color = :blue, label = "Liquid (cond/evap)") - lines!(ax1, T_range, S_ice; color = :magenta, label = "Ice (dep/subl)") - vlines!(ax1, [T_freeze]; color = :gray, linestyle = :dash, label = "T_freeze") - hlines!(ax1, [0]; color = :black, linewidth = 0.5) - text!(ax1, T_freeze + 2, minimum(S_ice) * 0.5; - text = "INP limiter\n(ice dep → 0)", - fontsize = 10, color = :gray, - ) - axislegend(ax1; position = :lb, framevisible = false) - - # Right panel: derivative vs T - ax2 = Axis(fig[1, 2]; - xlabel = "Temperature (K)", - ylabel = "∂(tendency)/∂q (1/s)", - title = "Derivative of tendency w.r.t. cloud condensate", - ) - lines!(ax2, T_range, dS_liq_full; color = :blue, label = "Liquid (full)") - lines!(ax2, T_range, dS_ice_full; color = :magenta, label = "Ice (full)") - # Simple derivative is -1/τ for both phases (identical when τ_liq == τ_ice) - lines!(ax2, T_range, dS_liq_simple; color = :black, linestyle = :dash, label = "Simple = −1/τ (both)") - vlines!(ax2, [T_freeze]; color = :gray, linestyle = :dash, label = "T_freeze") - text!(ax2, T_freeze + 2, minimum(dS_ice_full) * 0.5; - text = "INP limiter\n(∂ice → 0)", - fontsize = 10, color = :gray, - ) - axislegend(ax2; position = :lt, framevisible = false) - - save("condensation_evaporation_deriv.svg", fig) -end - -function make_deriv_vs_qvap_plot() - # Available water vapor range (q_vap = q_tot - q_lcl - q_icl - q_rai - q_sno) - # We vary q_tot while keeping q_lcl, q_icl fixed, so q_vap changes - qₗ = FT(1e-3) # 1 g/kg cloud liquid (fixed) - qᵢ = FT(1e-3) # 1 g/kg cloud ice (fixed) - qᵣ = FT(0) - qₛ = FT(0) - ρ = FT(1) # kg/m³ - - # q_tot range: from 2 g/kg (all condensate, no vapor) to 20 g/kg - qₜ_range = collect(range(FT(2e-3), stop = FT(20e-3), length = 181)) - qᵥ_range = qₜ_range .- qₗ .- qᵢ # available water vapor - - τ_relax = FT(1) - _liq = CM.Parameters.CloudLiquid(FT) - cm_liq = CM.Parameters.CloudLiquid(τ_relax, _liq.ρw, _liq.r_eff, _liq.N_0) - _ice = CM.Parameters.CloudIce(FT) - cm_ice = CM.Parameters.CloudIce(_ice.pdf, _ice.mass, _ice.r0, _ice.r_ice_snow, τ_relax, _ice.ρᵢ, _ice.r_eff, _ice.N_0) - - # Evaluate at two representative temperatures - T_vals = [FT(250), FT(280)] - T_labels = ["T=250K", "T=280K"] - T_colors_liq = [:royalblue, :cornflowerblue] - T_colors_ice = [:magenta, :hotpink] - - fig = Figure(size = (900, 450)) - - ax1 = Axis(fig[1, 1]; - xlabel = "Water vapor, qᵥ (g/kg)", - ylabel = "Tendency (1/s)", - title = "Tendency vs available water vapor\n(qₗ=1, qᵢ=1 g/kg, ρ=1 kg/m³, τ=1 s)", - ) - ax2 = Axis(fig[1, 2]; - xlabel = "Water vapor, qᵥ (g/kg)", - ylabel = "∂(tendency)/∂q (1/s)", - title = "Derivative vs available water vapor", - ) - - for (j, T) in enumerate(T_vals) - S_liq = zeros(FT, length(qₜ_range)) - dS_liq_f = zeros(FT, length(qₜ_range)) - dS_liq_s = zeros(FT, length(qₜ_range)) - S_ice = zeros(FT, length(qₜ_range)) - dS_ice_f = zeros(FT, length(qₜ_range)) - - for (i, qₜ) in enumerate(qₜ_range) - s = CMNe.conv_q_vap_to_q_lcl_icl_MM2015( - cm_liq, thp, qₜ, qₗ, qᵢ, qᵣ, qₛ, ρ, T, - ) - ds_s = CMNe.∂conv_q_vap_to_q_lcl_icl_MM2015_∂q_cld( - cm_liq, thp, qₜ, qₗ, qᵢ, qᵣ, qₛ, ρ, T, - ) - ds_f = CMNe.∂conv_q_vap_to_q_lcl_icl_MM2015_∂q_cld( - cm_liq, thp, qₜ, qₗ, qᵢ, qᵣ, qₛ, ρ, T; - simplified = false, - ) - S_liq[i] = s - dS_liq_f[i] = ds_f - dS_liq_s[i] = ds_s - - s = CMNe.conv_q_vap_to_q_lcl_icl_MM2015( - cm_ice, thp, qₜ, qₗ, qᵢ, qᵣ, qₛ, ρ, T, - ) - ds_f = CMNe.∂conv_q_vap_to_q_lcl_icl_MM2015_∂q_cld( - cm_ice, thp, qₜ, qₗ, qᵢ, qᵣ, qₛ, ρ, T; - simplified = false, - ) - S_ice[i] = s - dS_ice_f[i] = ds_f - end - - qᵥ_gkg = qᵥ_range * 1e3 - - # Left panel: tendency - lines!(ax1, qᵥ_gkg, S_liq; color = T_colors_liq[j], label = "Liquid $(T_labels[j])") - lines!(ax1, qᵥ_gkg, S_ice; color = T_colors_ice[j], label = "Ice $(T_labels[j])") - - # Right panel: derivative - lines!(ax2, qᵥ_gkg, dS_liq_f; color = T_colors_liq[j], label = "Liquid $(T_labels[j])") - lines!(ax2, qᵥ_gkg, dS_ice_f; color = T_colors_ice[j], label = "Ice $(T_labels[j])") - - if j == 1 - lines!(ax2, qᵥ_gkg, dS_liq_s; color = :black, linestyle = :dash, label = "Simple = −1/τ") - end - end - - hlines!(ax1, [0]; color = :black, linewidth = 0.5) - axislegend(ax1; position = :lt, framevisible = false) - axislegend(ax2; position = :lb, framevisible = false) - - save("condensation_evaporation_deriv_vs_qvap.svg", fig) -end - make_figure(HydrostaticBalance_qₗ_z()) make_figure(Range_qₗ_T()) -make_deriv_line_plot() -make_deriv_vs_qvap_plot() diff --git a/src/BulkMicrophysicsTendencies.jl b/src/BulkMicrophysicsTendencies.jl index ae4218f62..ecc97c959 100644 --- a/src/BulkMicrophysicsTendencies.jl +++ b/src/BulkMicrophysicsTendencies.jl @@ -434,7 +434,7 @@ and avoids quadratures over the derivatives. vel = mp.terminal_velocity # --- 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) + ∂tendency_∂q_lcl = -1 / CMNonEq.τ_relax(lcl) # Autoconversion q_lcl → q_rai (sink) 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) @@ -446,7 +446,7 @@ and avoids quadratures over the derivatives. max(q_lcl, UT.ϵ_numerics(typeof(q_lcl))) # --- 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) + ∂tendency_∂q_icl = -1 / CMNonEq.τ_relax(icl) # 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) @@ -519,7 +519,7 @@ derivatives; snow and cloud formation derivatives are zero for now. N_rai = ρ * n_rai # TODO: Cloud formation — 2M bulk_microphysics_tendencies does not call cloud formation yet; - # once it does, set these from MM2015 (same as 1M) via ∂conv_q_vap_to_q_lcl_icl_MM2015_∂q_cld. + # once it does, set these from MM2015 (same as 1M) via -1/τ_relax. ∂tendency_∂q_lcl = zero(ρ) ∂tendency_∂q_icl = zero(ρ) diff --git a/src/MicrophysicsNonEq.jl b/src/MicrophysicsNonEq.jl index 829cea931..231eae470 100644 --- a/src/MicrophysicsNonEq.jl +++ b/src/MicrophysicsNonEq.jl @@ -17,12 +17,9 @@ import ..Utilities as UT export τ_relax export conv_q_vap_to_q_lcl_icl export conv_q_vap_to_q_lcl_icl_MM2015 -export ∂conv_q_vap_to_q_lcl_icl_MM2015_∂q_cld export INP_limiter -export limit_MM2015_sinks export dqcld_dT export gamma_helper -export d2qcld_dT2 const CondEvap = Union{CMP.CloudLiquid, CMP.CondEvap2M} const SubDep = Union{CMP.CloudIce, CMP.SubDep2M} @@ -72,24 +69,6 @@ end return (q_sat_ice - q_icl) / τ_relax end -""" - limit_MM2015_sinks(tendency, q_cld) - -Checks if the cloud evaporation or sublimation tendency need to be -limited to prevent negative cloud condensate specific content. -The tendency is limited if it is negative and the cloud condensate specific -content is zero or negative. - -# Arguments -- `tendency` - cloud condensate tendency [kg/kg/s] -- `q_cld` - cloud liquid water or ice specific content [kg/kg] - -# Returns -- `true` if the tendency needs to be limited, `false` otherwise -""" -@inline function limit_MM2015_sinks(tendency, q_cld) - return tendency <= zero(tendency) && q_cld <= zero(q_cld) -end """ INP_limiter(tendency, tps, T) @@ -131,22 +110,6 @@ Computes the thermodynamic adjustment factor Γ. return 1 + (L / cₚ_air) * dqcld_dT end -""" - d2qcld_dT2(qᵥ_sat, L, Rᵥ, T) - -Computes the second derivative of the saturation specific humidity with respect to -temperature for a given phase of water. - -# Arguments -- `qᵥ_sat` - saturation specific humidity [kg/kg] -- `L` - latent heat [J/kg] -- `Rᵥ` - gas constant for water vapor [J/kg/K] -- `T` - temperature [K] -""" -@inline function d2qcld_dT2(qᵥ_sat, L, Rᵥ, T) - return qᵥ_sat * ((L / Rᵥ / T^2 - 1 / T)^2 + (1 / T^2 - 2 * L / Rᵥ / T^3)) -end - """ conv_q_vap_to_q_lcl_icl_MM2015(params, tps, q_tot, q_lcl, q_icl, q_rai, q_sno, ρ, T) @@ -179,6 +142,7 @@ function conv_q_vap_to_q_lcl_icl_MM2015( (; τ_relax)::CondEvap, tps::TDI.PS, q_tot, q_lcl, q_icl, q_rai, q_sno, ρ, T, ) + FT = eltype(tps) Rᵥ = TDI.Rᵥ(tps) Lᵥ = TDI.Lᵥ(tps, T) cₚ_air = TDI.cpₘ(tps, q_tot, q_lcl + q_rai, q_icl + q_sno) @@ -189,15 +153,21 @@ function conv_q_vap_to_q_lcl_icl_MM2015( dqsl_dT = dqcld_dT(qᵥ_sat_liq, Lᵥ, Rᵥ, T) Γₗ = gamma_helper(Lᵥ, cₚ_air, dqsl_dT) - # compute the tendency - tendency = (qᵥ - qᵥ_sat_liq) / (τ_relax * Γₗ) + sat_excess = qᵥ - qᵥ_sat_liq + timescale = τ_relax * Γₗ - return ifelse(limit_MM2015_sinks(tendency, q_lcl), zero(tendency), tendency) + # compute the tendency + return ifelse( + sat_excess < 0, + -max(FT(0), q_lcl) / timescale, + sat_excess / timescale, + ) end function conv_q_vap_to_q_lcl_icl_MM2015( (; τ_relax)::SubDep, tps::TDI.PS, q_tot, q_lcl, q_icl, q_rai, q_sno, ρ, T, ) + FT = eltype(tps) Rᵥ = TDI.Rᵥ(tps) Lₛ = TDI.Lₛ(tps, T) cₚ_air = TDI.cpₘ(tps, q_tot, q_lcl + q_rai, q_icl + q_sno) @@ -208,73 +178,19 @@ function conv_q_vap_to_q_lcl_icl_MM2015( dqsi_dT = dqcld_dT(qᵥ_sat_ice, Lₛ, Rᵥ, T) Γᵢ = gamma_helper(Lₛ, cₚ_air, dqsi_dT) - # compute the tendency - tendency = (qᵥ - qᵥ_sat_ice) / (τ_relax * Γᵢ) + sat_excess = qᵥ - qᵥ_sat_ice + timescale = τ_relax * Γᵢ - limiter = limit_MM2015_sinks(tendency, q_icl) || INP_limiter(tendency, tps, T) + # compute the tendency + tendency = ifelse( + sat_excess < 0, + -max(FT(0), q_icl) / timescale, + sat_excess / timescale, + ) + limiter = INP_limiter(tendency, tps, T) return ifelse(limiter, zero(tendency), tendency) end -""" - ∂conv_q_vap_to_q_lcl_icl_MM2015_∂q_cld(params::CloudLiquid, tps, q_tot, q_lcl, q_icl, q_rai, q_sno, ρ, T; simplified = true) - ∂conv_q_vap_to_q_lcl_icl_MM2015_∂q_cld(params::CloudIce, tps, q_tot, q_lcl, q_icl, q_rai, q_sno, ρ, T; simplified = true) - -Returns the derivative of the cloud condensate tendency with respect to the -corresponding cloud condensate species: - - `∂(tendency)/∂q_lcl` for the `CloudLiquid` dispatch - - `∂(tendency)/∂q_icl` for the `CloudIce` dispatch - -# Keyword Arguments -- `simplified` — if `true` (default), returns the leading-order approximation (`-1/τ_relax`); - if `false`, returns the total derivative accounting for implicit temperature feedback. -""" -function ∂conv_q_vap_to_q_lcl_icl_MM2015_∂q_cld( - lcl::CondEvap, tps::TDI.PS, - q_tot, q_lcl, q_icl, q_rai, q_sno, ρ, T; - simplified = true, -) - tendency = conv_q_vap_to_q_lcl_icl_MM2015(lcl, tps, q_tot, q_lcl, q_icl, q_rai, q_sno, ρ, T) - ∂tendency_∂q_lcl = -1 / lcl.τ_relax - if !simplified - Rᵥ = TDI.Rᵥ(tps) - Lᵥ = TDI.Lᵥ(tps, T) - cₚ_air = TDI.cpₘ(tps, q_tot, q_lcl + q_rai, q_icl + q_sno) - qᵥ_sat_liq = TDI.saturation_vapor_specific_content_over_liquid(tps, T, ρ) - - dqsl_dT = dqcld_dT(qᵥ_sat_liq, Lᵥ, Rᵥ, T) - Γₗ = gamma_helper(Lᵥ, cₚ_air, dqsl_dT) - d²qᵥ_sat_liq_dT² = d2qcld_dT2(qᵥ_sat_liq, Lᵥ, Rᵥ, T) - - ∂tendency_∂q_lcl -= tendency / Γₗ * (Lᵥ / cₚ_air)^2 * d²qᵥ_sat_liq_dT² - end - - return ifelse(limit_MM2015_sinks(tendency, q_lcl), zero(∂tendency_∂q_lcl), ∂tendency_∂q_lcl) -end -function ∂conv_q_vap_to_q_lcl_icl_MM2015_∂q_cld( - icl::SubDep, tps::TDI.PS, - q_tot, q_lcl, q_icl, q_rai, q_sno, ρ, T; - simplified = true, -) - tendency = conv_q_vap_to_q_lcl_icl_MM2015(icl, tps, q_tot, q_lcl, q_icl, q_rai, q_sno, ρ, T) - ∂tendency_∂q_icl = -1 / icl.τ_relax - if !simplified - Rᵥ = TDI.Rᵥ(tps) - Lₛ = TDI.Lₛ(tps, T) - cₚ_air = TDI.cpₘ(tps, q_tot, q_lcl + q_rai, q_icl + q_sno) - qᵥ_sat_ice = TDI.saturation_vapor_specific_content_over_ice(tps, T, ρ) - - dqsi_dT = dqcld_dT(qᵥ_sat_ice, Lₛ, Rᵥ, T) - Γᵢ = gamma_helper(Lₛ, cₚ_air, dqsi_dT) - d²qᵥ_sat_ice_dT² = d2qcld_dT2(qᵥ_sat_ice, Lₛ, Rᵥ, T) - - ∂tendency_∂q_icl -= tendency / Γᵢ * (Lₛ / cₚ_air)^2 * d²qᵥ_sat_ice_dT² - end - - limiter = limit_MM2015_sinks(tendency, q_icl) || INP_limiter(tendency, tps, T) - return ifelse(limiter, zero(∂tendency_∂q_icl), ∂tendency_∂q_icl) -end - - ### -------------------------- ### ### 1-moment terminal velocity ### ### -------------------------- ### diff --git a/test/gpu_tests.jl b/test/gpu_tests.jl index 39bc50675..3bea1c00f 100644 --- a/test/gpu_tests.jl +++ b/test/gpu_tests.jl @@ -83,15 +83,8 @@ end S_cond_MM2015 = CMN.conv_q_vap_to_q_lcl_icl_MM2015( lcl, tps, qᵥ_sl[i], q_lcl, q_icl, q_rai, q_sno, ρ[i], T[i], ) - S_cond_deriv_simple = CMN.∂conv_q_vap_to_q_lcl_icl_MM2015_∂q_cld( - lcl, tps, qᵥ_sl[i], q_lcl, q_icl, q_rai, q_sno, ρ[i], T[i], - ) - S_cond_deriv = CMN.∂conv_q_vap_to_q_lcl_icl_MM2015_∂q_cld( - lcl, tps, qᵥ_sl[i], q_lcl, q_icl, q_rai, q_sno, ρ[i], T[i]; - simplified = false, - ) S_cond = CMN.conv_q_vap_to_q_lcl_icl(icl, qᵢ_s[i], qᵢ[i]) - output[i] = (; S_cond_MM2015, S_cond_deriv_simple, S_cond_deriv, S_cond) + output[i] = (; S_cond_MM2015, S_cond) end @kernel inbounds = true function test_chen2022_terminal_velocity_kernel!( @@ -515,7 +508,7 @@ function test_gpu(FT) end TT.@testset "non-equilibrium microphysics kernels" begin - DT = @NamedTuple{S_cond_MM2015::FT, S_cond_deriv_simple::FT, S_cond_deriv::FT, S_cond::FT} + DT = @NamedTuple{S_cond_MM2015::FT, S_cond::FT} (; ndrange, output) = setup_output(1, DT) ρ = ArrayType([FT(0.8)]) @@ -526,11 +519,9 @@ function test_gpu(FT) kernel! = test_noneq_micro_kernel!(backend, work_groups) kernel!(lcl, icl, tps, output, ρ, T, qᵥ_sl, qᵢ, qᵢ_s; ndrange) - (; S_cond_MM2015, S_cond_deriv_simple, S_cond_deriv, S_cond) = Array(output)[1] + (; S_cond_MM2015, S_cond) = Array(output)[1] # test that nonequilibrium cloud formation is callable and returns a reasonable value TT.@test S_cond_MM2015 ≈ FT(3.76347635339803e-5) - TT.@test S_cond_deriv < FT(0) # derivative should be negative - TT.@test S_cond_deriv_simple ≈ FT(-0.1) # -1/τ_relax = -1/10 TT.@test S_cond ≈ FT(-1e-4) end diff --git a/test/microphysics_noneq_tests.jl b/test/microphysics_noneq_tests.jl index ca6898726..0b1f0266a 100644 --- a/test/microphysics_noneq_tests.jl +++ b/test/microphysics_noneq_tests.jl @@ -89,114 +89,14 @@ function test_microphysics_noneq(FT) # Liquid condensation above freezing should be unaffected TT.@test CMNe.conv_q_vap_to_q_lcl_icl_MM2015(liquid, tps, FT(1.5 * qᵥ_sl_w), FT(0), FT(0), FT(0), FT(0), ρ, T_warm) > FT(0) # Ice sublimation (negative tendency) above freezing should NOT be limited + # When subsaturated w.r.t. ice and q_icl > 0, the tendency is negative (sublimation). + # The INP limiter only zeroes positive tendencies above freezing, so sublimation passes through. TT.@test CMNe.conv_q_vap_to_q_lcl_icl_MM2015(ice, tps, FT(0.5 * qᵥ_si_w), FT(0), FT(0.001), FT(0), FT(0), ρ, T_warm) < FT(0) #! format: on end - TT.@testset "CondEvap_DepSub_MM2015_derivatives" begin - ρ = FT(0.8) - T_freeze = TDI.T_freeze(tps) - T = FT(T_freeze - 10) - q_tot = FT(0.01) # 10 g/kg total water - q_lcl = FT(0.001) # 1 g/kg cloud liquid - q_icl = FT(0.0005) # 0.5 g/kg cloud ice - q_rai = FT(0) - q_sno = FT(0) - τ = CMNe.τ_relax(liquid) - - # --- Liquid condensation derivative --- - S_liq = CMNe.conv_q_vap_to_q_lcl_icl_MM2015( - liquid, tps, q_tot, q_lcl, q_icl, q_rai, q_sno, ρ, T, - ) - dS_liq_simple = CMNe.∂conv_q_vap_to_q_lcl_icl_MM2015_∂q_cld( - liquid, tps, q_tot, q_lcl, q_icl, q_rai, q_sno, ρ, T, - ) - dS_liq = CMNe.∂conv_q_vap_to_q_lcl_icl_MM2015_∂q_cld( - liquid, tps, q_tot, q_lcl, q_icl, q_rai, q_sno, ρ, T; - simplified = false, - ) - - # derivative should be negative (increasing q_lcl decreases tendency) - TT.@test dS_liq < FT(0) - TT.@test dS_liq_simple < FT(0) - - # simple derivative should be -1/τ - TT.@test dS_liq_simple ≈ FT(-1) / τ - - # full derivative should be close to but not equal to simple - TT.@test dS_liq ≈ dS_liq_simple rtol = 0.5 - TT.@test dS_liq != dS_liq_simple - - # --- Ice deposition derivative --- - S_ice = CMNe.conv_q_vap_to_q_lcl_icl_MM2015( - ice, tps, q_tot, q_lcl, q_icl, q_rai, q_sno, ρ, T, - ) - dS_ice_simple = CMNe.∂conv_q_vap_to_q_lcl_icl_MM2015_∂q_cld( - ice, tps, q_tot, q_lcl, q_icl, q_rai, q_sno, ρ, T, - ) - dS_ice = CMNe.∂conv_q_vap_to_q_lcl_icl_MM2015_∂q_cld( - ice, tps, q_tot, q_lcl, q_icl, q_rai, q_sno, ρ, T; - simplified = false, - ) - - TT.@test dS_ice < FT(0) - TT.@test dS_ice_simple ≈ FT(-1) / CMNe.τ_relax(ice) - - # --- Finite difference validation at constant e_tot, ρ, q_tot --- - # Perturb q_cld, recover T from energy conservation, - # then check analytical derivative vs finite difference. - if FT == Float64 - ε = FT(1e-8) - q_liq = q_lcl + q_rai - q_ice = q_icl + q_sno - e_int = TDI.TD.internal_energy(tps, T, q_tot, q_liq, q_ice) - - # Validate ∂S/∂q_lcl for liquid at constant e_tot - q_liq_p = (q_lcl + ε) + q_rai - q_liq_m = (q_lcl - ε) + q_rai - T_p = TDI.TD.air_temperature(tps, e_int, q_tot, q_liq_p, q_ice) - T_m = TDI.TD.air_temperature(tps, e_int, q_tot, q_liq_m, q_ice) - S_plus = CMNe.conv_q_vap_to_q_lcl_icl_MM2015( - liquid, tps, q_tot, q_lcl + ε, q_icl, q_rai, q_sno, ρ, T_p, - ) - S_minus = CMNe.conv_q_vap_to_q_lcl_icl_MM2015( - liquid, tps, q_tot, q_lcl - ε, q_icl, q_rai, q_sno, ρ, T_m, - ) - dS_fd_liq = (S_plus - S_minus) / (2ε) - TT.@test dS_liq ≈ dS_fd_liq rtol = 0.2 - - # Validate ∂S/∂q_icl for ice at constant e_tot - q_ice_p = q_icl + ε + q_sno - q_ice_m = q_icl - ε + q_sno - T_p_i = TDI.TD.air_temperature(tps, e_int, q_tot, q_liq, q_ice_p) - T_m_i = TDI.TD.air_temperature(tps, e_int, q_tot, q_liq, q_ice_m) - S_plus_i = CMNe.conv_q_vap_to_q_lcl_icl_MM2015( - ice, tps, q_tot, q_lcl, q_icl + ε, q_rai, q_sno, ρ, T_p_i, - ) - S_minus_i = CMNe.conv_q_vap_to_q_lcl_icl_MM2015( - ice, tps, q_tot, q_lcl, q_icl - ε, q_rai, q_sno, ρ, T_m_i, - ) - dS_fd_ice = (S_plus_i - S_minus_i) / (2ε) - TT.@test dS_ice ≈ dS_fd_ice rtol = 0.2 - end - - # --- Limiter: when tendency < 0 and q <= 0, all outputs should be zero --- - S_z = CMNe.conv_q_vap_to_q_lcl_icl_MM2015( - liquid, tps, FT(0.5 * 0.002), FT(0), FT(0), FT(0), FT(0), ρ, T, - ) - dS_z = CMNe.∂conv_q_vap_to_q_lcl_icl_MM2015_∂q_cld( - liquid, tps, FT(0.5 * 0.002), FT(0), FT(0), FT(0), FT(0), ρ, T, - ) - dS_z_f = CMNe.∂conv_q_vap_to_q_lcl_icl_MM2015_∂q_cld( - liquid, tps, FT(0.5 * 0.002), FT(0), FT(0), FT(0), FT(0), ρ, T; - simplified = false, - ) - TT.@test S_z == FT(0) - TT.@test dS_z == FT(0) - TT.@test dS_z_f == FT(0) - end TT.@testset "Cloud condensate sedimentation - liquid" begin #setup diff --git a/test/performance_tests.jl b/test/performance_tests.jl index 74dac8b19..1d4d77f1e 100644 --- a/test/performance_tests.jl +++ b/test/performance_tests.jl @@ -176,7 +176,7 @@ function benchmark_test(FT) bench_press(FT, P3.get_distribution_logλ, (params_P3, L_ice, N_ice, F_rim, ρ_rim), 30_000) bench_press(FT, P3.ice_terminal_velocity_number_weighted, (ch2022, ρ_air, state, logλ), 120_000) bench_press(FT, P3.ice_terminal_velocity_mass_weighted, (ch2022, ρ_air, state, logλ), 135_000) - bench_press(FT, P3.integrate, (x -> x^4, FT(0), FT(1)), 3_300) + bench_press(FT, P3.integrate, (x -> x^4, FT(0), FT(1)), 3_700) bench_press(FT, P3.D_m, (state, logλ), 3_000) @info "P3 Ice Nucleation"