From 8cae7dc8525ffb7b40c2bb96a3c0d1a85d3b90da Mon Sep 17 00:00:00 2001 From: Anna Jaruga Date: Wed, 12 Jun 2024 11:14:35 -0700 Subject: [PATCH] Cleanup in limiters --- src/Microphysics2M.jl | 81 ++++++++++++++++++------------------ test/Project.toml | 3 +- test/microphysics2M_tests.jl | 71 ++++++++++++++----------------- 3 files changed, 73 insertions(+), 82 deletions(-) diff --git a/src/Microphysics2M.jl b/src/Microphysics2M.jl index 6f5b794b8..7d81200b0 100644 --- a/src/Microphysics2M.jl +++ b/src/Microphysics2M.jl @@ -66,11 +66,11 @@ function pdf_cloud( ) where {FT} (; νc, μc, ρw) = pdf_c - ϵ_N = FT(10) # TODO - should it be a parameter + #ϵ_N = FT(10) # TODO - should it be a parameter χ = 9 # convert mean droplet mass from kg to μg - if qₗ < eps(FT) || Nₗ < ϵ_N + if qₗ < eps(FT) || Nₗ < eps(FT) return ( xc = FT(0), χ, @@ -116,9 +116,9 @@ function pdf_rain( Nᵣ::FT, ) where {FT} (; νr, μr, ρw) = pdf_r - ϵ_N = FT(10) # TODO - should it be a parameter? + #ϵ_N = FT(10) # TODO - should it be a parameter? - if qᵣ < eps(FT) || Nᵣ < ϵ_N + if qᵣ < eps(FT) || Nᵣ < eps(FT) return (λr = FT(0), xr = FT(0), Ar = FT(0), Br = FT(0)) else λr = (Nᵣ * FT(π) * ρw / ρₐ / qᵣ)^FT(1 / 3) @@ -133,24 +133,24 @@ function pdf_rain( end function pdf_rain( pdf_r::CMP.RainParticlePDF_SB2006_limited{FT}, - q_rai::FT, - ρ::FT, - N_rai::FT, + qᵣ::FT, + ρₐ::FT, + Nᵣ::FT, ) where {FT} (; νr, μr, xr_min, xr_max, N0_min, N0_max, λ_min, λ_max, ρw) = pdf_r - L_rai = ρ * q_rai - xr_0 = L_rai / N_rai + qᵣ = qᵣ > FT(0) ? qᵣ : FT(0) + + Lᵣ = ρₐ * qᵣ + xr_0 = Nᵣ > eps(FT) ? Lᵣ / Nᵣ : FT(0) xr_hat = max(xr_min, min(xr_max, xr_0)) - N0 = max(N0_min, min(N0_max, N_rai * (FT(π) * ρw / xr_hat)^FT(1 / 3))) - λr = max(λ_min, min(λ_max, (FT(π) * ρw * N0 / L_rai)^FT(1 / 4))) - xr = max(xr_min, min(xr_max, L_rai * λr / N0)) - αr = N_rai * λr + N0 = max(N0_min, min(N0_max, Nᵣ * (FT(π) * ρw / xr_hat)^FT(1 / 3))) + λr = max(λ_min, min(λ_max, (FT(π) * ρw * N0 / Lᵣ)^FT(1 / 4))) + xr = max(xr_min, min(xr_max, Lᵣ * λr / N0)) + αr = Nᵣ * λr - Br = - (xr == 0) ? FT(0) : - (SF.gamma(FT(νr + 1) / μr) * xr / SF.gamma(FT(νr + 2) / μr))^(-μr) - Ar = μr * N_rai * Br^(FT(νr + 1) / μr) / SF.gamma(FT(νr + 1) / μr) + Br = (SF.gamma(FT(νr + 1) / μr) * xr / SF.gamma(FT(νr + 2) / μr))^(-μr) + Ar = μr * Nᵣ * Br^(FT(νr + 1) / μr) / SF.gamma(FT(νr + 1) / μr) return (; λr, αr, xr, Ar, Br) end @@ -176,7 +176,7 @@ function autoconversion( N_liq, ) where {FT} - if q_liq < eps(FT) + if q_liq < eps(FT) || N_liq < eps(FT) return LiqRaiRates{FT}() end @@ -221,7 +221,7 @@ collisions between raindrops and cloud droplets (accretion) for `scheme == SB200 """ function accretion((; accr)::CMP.SB2006{FT}, q_liq, q_rai, ρ, N_liq) where {FT} - if q_liq < eps(FT) || q_rai < eps(FT) + if q_liq < eps(FT) || q_rai < eps(FT) || N_liq < eps(FT) return LiqRaiRates{FT}() end @@ -326,7 +326,7 @@ function rain_self_collection( N_rai::FT, ) where {FT} - if q_rai < eps(FT) + if q_rai < eps(FT) || N_rai < eps(FT) return FT(0) end @@ -334,11 +334,7 @@ function rain_self_collection( (; ρ0, ρw) = pdf L_rai = ρ * q_rai - #λr = - # pdf_rain(pdf, q_rai, ρ, N_rai).λr * - # (SF.gamma(FT(4)) / FT(π) / ρw)^FT(1 / 3) λr = pdf_rain(pdf, q_rai, ρ, N_rai).Br - dN_rai_dt_sc = -krr * N_rai * L_rai * sqrt(ρ0 / ρ) * (1 + κrr / λr)^d return dN_rai_dt_sc @@ -368,7 +364,7 @@ function rain_breakup( dN_rai_dt_sc::FT, ) where {FT} - if q_rai < eps(FT) + if q_rai < eps(FT) || N_rai < eps(FT) return FT(0) end (; Deq, Dr_th, kbr, κbr) = brek @@ -431,12 +427,15 @@ function rain_terminal_velocity( N_rai::FT, ) where {FT} # TODO: Input argument list needs to be redesigned - if q_rai < eps(FT) - return (FT(0), FT(0)) - end + λr = pdf_rain(pdf_r, q_rai, ρ, N_rai).λr - vt0 = max(FT(0), sqrt(ρ0 / ρ) * (aR - bR / (1 + cR / λr))) - vt1 = max(FT(0), sqrt(ρ0 / ρ) * (aR - bR / (1 + cR / λr)^FT(4))) + + vt0 = + N_rai < eps(FT) ? FT(0) : + max(FT(0), sqrt(ρ0 / ρ) * (aR - bR / (1 + cR / λr))) + vt1 = + q_rai < eps(FT) ? FT(0) : + max(FT(0), sqrt(ρ0 / ρ) * (aR - bR / (1 + cR / λr)^FT(4))) return (vt0, vt1) end function rain_terminal_velocity( @@ -446,9 +445,6 @@ function rain_terminal_velocity( ρ::FT, N_rai::FT, ) where {FT} - if q_rai < eps(FT) - return (FT(0), FT(0)) - end # coefficients from Table B1 from Chen et. al. 2022 aiu, bi, ciu = CO.Chen2022_vel_coeffs_small(vel, ρ) # size distribution parameter @@ -458,8 +454,8 @@ function rain_terminal_velocity( vt0 = sum(CO.Chen2022_vel_add.(aiu, bi, ciu, λ, 0)) vt3 = sum(CO.Chen2022_vel_add.(aiu, bi, ciu, λ, 3)) - vt0 = max(FT(0), vt0) - vt3 = max(FT(0), vt3) + vt0 = N_rai < eps(FT) ? FT(0) : max(FT(0), vt0) + vt3 = q_rai < eps(FT) ? FT(0) : max(FT(0), vt3) # It should be (ϕ^κ * vt0, ϕ^κ * vt3), but for rain drops ϕ = 1 and κ = 0 return (vt0, vt3) end @@ -507,7 +503,7 @@ function rain_evaporation( evap_rate_1 = FT(0) S = TD.supersaturation(tps, q, ρ, T, TD.Liquid()) - if (q_rai > FT(0) && S < FT(0)) + if ((q_rai > eps(FT) || N_rai > eps(FT)) && S < FT(0)) (; ν_air, D_vapor) = aps (; av, bv, α, β, ρ0) = evap @@ -534,6 +530,9 @@ function rain_evaporation( evap_rate_0 = min(FT(0), FT(2) * FT(π) * G * S * N_rai * Dr * Fv0 / xr) evap_rate_1 = min(FT(0), FT(2) * FT(π) * G * S * N_rai * Dr * Fv1 / ρ) + + evap_rate_0 = N_rai < eps(FT) ? FT(0) : evap_rate_0 + evap_rate_1 = q_rai < eps(FT) ? FT(0) : evap_rate_1 end return (; evap_rate_0, evap_rate_1) @@ -578,19 +577,19 @@ function radar_reflectivity( χ = pdf_cloud(pdf_c, q_liq, ρ_air, N_liq).χ Zc = - Bc == 0 ? FT(0) : + Bc < eps(FT) ? FT(0) : Ac * C^(νc + 1) * (Bc * C^μc)^(-(3 + νc) / μc) * SF.gamma((3 + νc) / μc) / μc * FT(10)^(-2 * χ) Zr = - Br == 0 ? FT(0) : + Br < eps(FT) ? FT(0) : Ar * C^(νr + 1) * (Br * C^μr)^(-(3 + νr) / μr) * SF.gamma((3 + νr) / μr) / μr - return Zc + Zr == 0 ? FT(-135) : 10 * (log10(Zc + Zr) - log_10_Z₀) + return max(FT(-150), 10 * (log10(Zc + Zr) - log_10_Z₀)) end """ @@ -655,7 +654,7 @@ function effective_radius( (Br * C^μr)^(-(5 + 3 * νr) / (3 * μr)) * SF.gamma((5 + 3 * νr) / (3 * μr)) / μr - return M2_c + M2_r == 0 ? FT(0) : (M3_c + M3_r) / (M2_c + M2_r) + return M2_c + M2_r <= eps(FT) ? FT(0) : (M3_c + M3_r) / (M2_c + M2_r) end """ @@ -681,7 +680,7 @@ function effective_radius_Liu_Hallet_97( k = FT(0.8) r_vol = - ((N_liq + N_rai) == 0) ? FT(0) : + ((N_liq + N_rai) < eps(FT)) ? FT(0) : ( (FT(3) * (q_liq + q_rai) * ρ_air) / (FT(4) * π * ρ_w * (N_liq + N_rai)) diff --git a/test/Project.toml b/test/Project.toml index 242ca35af..5b76de58a 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -13,6 +13,7 @@ EnsembleKalmanProcesses = "aa8a2aa5-91d8-4396-bcef-d4f2ec43552d" EvoTrees = "f6006082-12f8-11e9-0c9c-0d5d367ab1e5" Flux = "587475ba-b771-5e3f-ad9e-33799f191a9c" GaussianProcesses = "891a1506-143c-57d2-908e-e1f8e92e6de9" +HCubature = "19dc6840-f33b-545b-b366-655c7e3ffd49" JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b" KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" @@ -21,8 +22,8 @@ MLJFlux = "094fc8d1-fd35-5302-93ea-dabda2abf845" MLJModels = "d491faf4-2d78-11e9-2867-c94bc002c0b7" Optim = "429524aa-4258-5aef-a3af-852621145aeb" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" -Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" RootSolvers = "7181ea78-2dcb-4de3-ab41-2b8ab5a31e74" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" diff --git a/test/microphysics2M_tests.jl b/test/microphysics2M_tests.jl index 710f8a9f7..6853abfdc 100644 --- a/test/microphysics2M_tests.jl +++ b/test/microphysics2M_tests.jl @@ -150,8 +150,8 @@ function test_microphysics2M(FT) # 2M_microphysics - Seifert and Beheng 2006 double moment scheme tests TT.@testset "limiting lambda_r and x_r - Seifert and Beheng 2006" begin #setup - q_rai = [FT(0), FT(1e-4), FT(1e-2)] - N_rai = [FT(1e1), FT(1e3), FT(1e5)] + q_rai = [FT(0), FT(1e-3), FT(1e-4), FT(1e-2)] + N_rai = [FT(1e1), FT(1e1), FT(1e3), FT(1e5)] ρ = FT(1) (; xr_min, xr_max, λ_min, λ_max) = SB2006.pdf_r @@ -162,7 +162,6 @@ function test_microphysics2M(FT) λ = CM2.pdf_rain(SB2006.pdf_r, qr, ρ, Nr).λr xr = CM2.pdf_rain(SB2006.pdf_r, qr, ρ, Nr).xr - #test TT.@test λ_min <= λ <= λ_max TT.@test xr_min <= xr <= xr_max end @@ -322,9 +321,8 @@ function test_microphysics2M(FT) sc_br_rai = CM2.rain_self_collection_and_breakup(SB, q_rai, ρ, N_rai) - λr = - CM2.pdf_rain(SB.pdf_r, q_rai, ρ, N_rai).λr * - FT(6 / π / 1000)^FT(1 / 3) + λr = CM2.pdf_rain(SB.pdf_r, q_rai, ρ, N_rai).Br + dNrdt_sc = -krr * N_rai * ρ * q_rai * (1 + κrr / λr)^-5 * sqrt(ρ0 / ρ) @@ -394,12 +392,13 @@ function test_microphysics2M(FT) TT.@test vt_rai isa Tuple TT.@test vt_rai[1] ≈ vt0 rtol = 1e-6 TT.@test vt_rai[2] ≈ vt1 rtol = 1e-6 + TT.@test CM2.rain_terminal_velocity( SB, SB2006Vel, - FT(0), + q_rai, ρ, - N_rai, + FT(0), )[1] ≈ 0 atol = eps(FT) TT.@test CM2.rain_terminal_velocity( SB, @@ -432,9 +431,9 @@ function test_microphysics2M(FT) TT.@test CM2.rain_terminal_velocity( SB, Chen2022Vel, - FT(0), + q_rai, ρ, - N_rai, + FT(0), )[1] ≈ 0 atol = eps(FT) TT.@test CM2.rain_terminal_velocity( SB, @@ -492,9 +491,9 @@ function test_microphysics2M(FT) aps, tps, q, - FT(0), + q_rai, ρ, - N_rai, + FT(0), T, ).evap_rate_0 ≈ 0 atol = eps(FT) TT.@test CM2.rain_evaporation( @@ -510,38 +509,30 @@ function test_microphysics2M(FT) end end - TT.@testset "2M_microphysics - Seifert and Beheng 2006 radar reflectivity" begin + TT.@testset "2M_microphysics - Seifert and Beheng 2006 effective radius and reflectivity" begin #setup - ρ_air = FT(1) - q_liq = FT(2.128e-4) - N_liq = FT(15053529) - q_rai = FT(1.573e-4) - N_rai = FT(510859) + ρₐ = FT(1) - for SB in [SB2006, SB2006_no_limiters] - #action - rr = CM2.radar_reflectivity(SB, q_liq, q_rai, N_liq, N_rai, ρ_air) + q_liq = [FT(2.128e-4), FT(2.128e-20), FT(1.6e-12), FT(0), FT(1.037e-25)] + N_liq = [FT(15053529), FT(3), FT(5512), FT(0), FT(5.225e-12)] + q_rai = [FT(1.573e-4), FT(1.573e-4), FT(1.9e-15), FT(0), FT(2.448e-27)] + N_rai = [FT(510859), FT(510859), FT(0), FT(0), FT(5.136e-18)] - # test - TT.@test rr ≈ FT(-13) atol = FT(0.5) - end - end + # reference values + rr = [FT(-12.561951), FT(-12.579899), FT(-150), FT(-150), FT(-150)] + reff = [FT(2.319383e-5), FT(6.91594e-5), FT(0), FT(0), FT(0)] - TT.@testset "2M_microphysics - Seifert and Beheng 2006 effective radius" begin - #setup - ρ_air = FT(1) - ρ_w = FT(1000) - q_liq = FT(2.128e-4) - N_liq = FT(15053529) - q_rai = FT(1.573e-4) - N_rai = FT(510859) + for (qₗ, Nₗ, qᵣ, Nᵣ, rₑ, Z) in zip(q_liq, N_liq, q_rai, N_rai, reff, rr) + for SB in [SB2006, SB2006_no_limiters] - for SB in [SB2006, SB2006_no_limiters] - #action - reff = CM2.effective_radius(SB, q_liq, q_rai, N_liq, N_rai, ρ_air) + #action + Z_val = CM2.radar_reflectivity(SB, qₗ, qᵣ, Nₗ, Nᵣ, ρₐ) + rₑ_val = CM2.effective_radius(SB, qₗ, qᵣ, Nₗ, Nᵣ, ρₐ) - #test - TT.@test reff ≈ FT(2.66e-05) atol = FT(8e-6) + #test + TT.@test rₑ_val ≈ rₑ atol = FT(1e-6) + TT.@test Z_val ≈ Z atol = FT(1e-4) + end end end @@ -624,7 +615,7 @@ function test_microphysics2M(FT) TT.@test qD ≈ qᵣ atol = FT(1e-6) TT.@test qx ≈ qᵣ atol = FT(1e-6) - # The mass integrals don't work with limiters + # The mass integrals don't work with limiters TT.@test qx_lim ≈ qᵣ atol = FT(1e-6) TT.@test qD_lim ≈ qx_lim atol = FT(1e-6) end @@ -665,7 +656,7 @@ function test_microphysics2M(FT) TT.@test qx ≈ qₗ rtol = FT(1e-6) TT.@test Nx ≈ Nₗ rtol = FT(1e-6) - # mass of liquid droplets as a function of its diameter in mm and μg + # mass of liquid droplets as a function of its diameter in mm and μg _m(D) = FT(π / 6) * ρₗ * FT(10)^(χ - 9) * D^3 # integral bounds in millimiters