From 01ca4732d49af6fbf6467c0ce93750acecc4512e Mon Sep 17 00:00:00 2001 From: Anastasia Popova Date: Mon, 22 Apr 2024 15:51:15 -0700 Subject: [PATCH] Add P3 terminal velocity functions --- docs/src/API.md | 7 +- docs/src/plots/P3TerminalVelocityPlots.jl | 2 +- src/CloudMicrophysics.jl | 3 +- src/Microphysics2M.jl | 22 ++ src/P3.jl | 6 +- src/P3_terminal_velocity.jl | 121 ++++++++++- src/TerminalVelocity.jl | 79 ------- test/p3_tests.jl | 246 +++++++++++++++++++++- test/performance_tests.jl | 2 +- test/runtests.jl | 1 - test/terminal_velocity_tests.jl | 55 ----- 11 files changed, 383 insertions(+), 161 deletions(-) delete mode 100644 src/TerminalVelocity.jl delete mode 100644 test/terminal_velocity_tests.jl diff --git a/docs/src/API.md b/docs/src/API.md index facea6e98..61ad28e5b 100644 --- a/docs/src/API.md +++ b/docs/src/API.md @@ -66,12 +66,7 @@ Microphysics2M.conv_q_liq_to_q_rai P3Scheme P3Scheme.thresholds P3Scheme.distribution_parameter_solver -``` - -# Terminal Velocity -```@docs -TerminalVelocity -TerminalVelocity.velocity_chen +P3Scheme.ice_terminal_velocity ``` # Aerosol model diff --git a/docs/src/plots/P3TerminalVelocityPlots.jl b/docs/src/plots/P3TerminalVelocityPlots.jl index ce3f2c3ce..3ff20f9dc 100644 --- a/docs/src/plots/P3TerminalVelocityPlots.jl +++ b/docs/src/plots/P3TerminalVelocityPlots.jl @@ -31,7 +31,7 @@ function get_values( ρ_r = ρ_rs[j] V_m[i, j] = - P3.terminal_velocity(p3, Chen2022, q, N, ρ_r, F_r, ρ_a)[2] + P3.ice_terminal_velocity(p3, Chen2022, q, N, ρ_r, F_r, ρ_a)[2] # get D_m in mm for plots D_m[i, j] = 1e3 * P3.D_m(p3, q, N, ρ_r, F_r) end diff --git a/src/CloudMicrophysics.jl b/src/CloudMicrophysics.jl index b1a25d2c2..5688c8157 100644 --- a/src/CloudMicrophysics.jl +++ b/src/CloudMicrophysics.jl @@ -14,13 +14,12 @@ include("Common.jl") include("Microphysics0M.jl") include("Microphysics1M.jl") include("Microphysics2M.jl") -include("TerminalVelocity.jl") +include("IceNucleation.jl") include("P3.jl") include("MicrophysicsNonEq.jl") include("AerosolModel.jl") include("AerosolActivation.jl") include("Nucleation.jl") -include("IceNucleation.jl") include("PrecipitationSusceptibility.jl") include("ArtifactCalling.jl") diff --git a/src/Microphysics2M.jl b/src/Microphysics2M.jl index dd8f248f1..6605a6186 100644 --- a/src/Microphysics2M.jl +++ b/src/Microphysics2M.jl @@ -504,6 +504,28 @@ function rain_self_collection_and_breakup( return (; sc, br) end +""" + rain_particle_terminal_velocity(D, Chen2022, ρₐ) + + - D - maximum particle dimension + - Chen2022 - a struct with terminal velocity parameters from Chen 2022 + - ρₐ - air density + +Returns the terminal velocity of an individual rain drop as a function +of its size (maximum dimension) following Chen 2022 velocity parametrization. +Needed for numerical integrals in the P3 scheme. +""" +function rain_particle_terminal_velocity( + D::FT, + Chen2022::CMP.Chen2022VelTypeRain, + ρₐ::FT, +) where {FT} + (ai, bi, ci) = CO.Chen2022_vel_coeffs_small(Chen2022, ρₐ) + + v = sum(@. sum(ai * D^bi * exp(-ci * D))) + return v +end + """ rain_terminal_velocity(SB2006, vel, q_rai, ρ, N_rai) diff --git a/src/P3.jl b/src/P3.jl index 40fab104b..38b8a9cdd 100644 --- a/src/P3.jl +++ b/src/P3.jl @@ -14,7 +14,11 @@ import HCubature as HC import ClimaParams as CP import CloudMicrophysics.Parameters as CMP import CloudMicrophysics.Common as CO -import CloudMicrophysics.TerminalVelocity as TV +import CloudMicrophysics.HetIceNucleation as CM_HetIce +import CloudMicrophysics.Microphysics2M as CM2 + +import Thermodynamics as TD +import Thermodynamics.Parameters as TDP const PSP3 = CMP.ParametersP3 diff --git a/src/P3_terminal_velocity.jl b/src/P3_terminal_velocity.jl index cab19c71b..000702809 100644 --- a/src/P3_terminal_velocity.jl +++ b/src/P3_terminal_velocity.jl @@ -1,5 +1,53 @@ """ - terminal_velocity(p3, Chen2022, q, N, ρ_r, F_r, ρ_a) + ϕᵢ(mᵢ, aᵢ, ρᵢ) + + - mᵢ - particle mass + - aᵢ - particle area + - ρᵢ - ice density + +Returns the aspect ratio (ϕ) for an ice particle with given mass, area, and ice density. +""" +function ϕᵢ(mᵢ::FT, aᵢ::FT, ρᵢ::FT) where {FT} + # TODO - add some notes on how we derived it + # TODO - should we make use of other P3 properties like rimed fraction and volume? + return 16 * ρᵢ^2 * aᵢ^3 / (9 * π * mᵢ^2) +end + +""" + ice_particle_terminal_velocity(D, Chen2022, ρₐ, mᵢ, aᵢ, ρᵢ) + + - D - maximum particle dimension + - Chen2022 - a struct with terminal velocity parameters from Chen 2022 + - ρₐ - air density + - mᵢ - particle mass + - aᵢ - particle area + - ρᵢ - ice density + +Returns the terminal velocity of a single ice particle as a function +of its size (maximum dimension) using Chen 2022 parametrization. +Needed for numerical integrals in the P3 scheme. +""" +function ice_particle_terminal_velocity( + D::FT, + Chen2022::CMP.Chen2022VelTypeSnowIce, + ρₐ::FT, + mᵢ::FT, + aᵢ::FT, + ρᵢ::FT, +) where {FT} + if D <= Chen2022.cutoff + (ai, bi, ci) = CO.Chen2022_vel_coeffs_small(Chen2022, ρₐ) + else + (ai, bi, ci) = CO.Chen2022_vel_coeffs_large(Chen2022, ρₐ) + end + v = sum(@. sum(ai * D^bi * exp(-ci * D))) + + κ = FT(-1 / 6) + return ϕᵢ(mᵢ, aᵢ, ρᵢ)^κ * v +end + +""" + ice_terminal_velocity(p3, Chen2022, q, N, ρ_r, F_r, ρ_a) - p3 - a struct with P3 scheme parameters - Chen2022 - a struch with terminal velocity parameters as in Chen(2022) @@ -9,10 +57,10 @@ - F_r - rime mass fraction (q_rim/q_i) - ρ_a - density of air - Returns the mass and number weighted fall speeds - Eq C10 of Morrison and Milbrandt (2015) +Returns the mass and number weighted fall speeds for ice following +eq C10 of Morrison and Milbrandt (2015). """ -function terminal_velocity( +function ice_terminal_velocity( p3::PSP3, Chen2022::CMP.Chen2022VelTypeSnowIce, q::FT, @@ -52,6 +100,7 @@ function terminal_velocity( spheres_n(a, b, c) = (aₛ(a), bₛ(b), cₛ(c)) spheres_m(a, b, c) = (aₛ_m(a), bₛ_m(b), cₛ(c)) + # TODO - combine different aspect ratio computations aₙₛ(a) = aₛ(a) * (16 * p3.ρ_i^2 * p3.γ^3 / (9 * FT(π) * α_va^2))^κ bₙₛ(b) = bₛ(b) + κ * (3 * p3.σ - 2 * p3.β_va) @@ -183,3 +232,67 @@ function terminal_velocity( end return (v_n / N, v_m / q) end + +""" + velocity_difference(type, Dₗ, Dᵢ, p3, Chen2022, ρ_a, F_r, th) + + - type - a struct containing the size distribution parameters of the particle colliding with ice + - Dₗ - maximum dimension of the particle colliding with ice + - Dᵢ - maximum dimension of ice particle + - p3 - a struct containing P3 parameters + - Chen2022 - a struct containing Chen 2022 velocity parameters + - ρ_a - density of air + - F_r - rime mass fraction (q_rim/ q_i) + - th - P3 particle properties thresholds + +Returns the absolute value of the velocity difference between an ice particle and +cloud or rain drop as a function of their sizes. It uses Chen 2022 velocity +parameterization for ice and rain and assumes no sedimentation of cloud droplets. +Needed for P3 numerical integrals +""" +function velocity_difference( + pdf_r::Union{ + CMP.RainParticlePDF_SB2006{FT}, + CMP.RainParticlePDF_SB2006_limited{FT}, + }, + Dₗ::FT, + Dᵢ::FT, + p3::PSP3, + Chen2022::CMP.Chen2022VelType, + ρ_a::FT, + F_r::FT, + th, +) where {FT} + return abs( + ice_particle_terminal_velocity( + Dᵢ, + Chen2022.snow_ice, + ρ_a, + p3_mass(p3, Dᵢ, F_r, th), + p3_area(p3, Dᵢ, F_r, th), + p3.ρ_i, + ) - CM2.rain_particle_terminal_velocity(Dₗ, Chen2022.rain, ρ_a), + ) +end +function velocity_difference( + pdf_c::CMP.CloudParticlePDF_SB2006{FT}, + Dₗ::FT, + Dᵢ::FT, + p3::PSP3, + Chen2022::CMP.Chen2022VelType, + ρ_a::FT, + F_r::FT, + th, +) where {FT} + # velocity difference for cloud collisions + return abs( + ice_particle_terminal_velocity( + Dᵢ, + Chen2022.snow_ice, + ρ_a, + p3_mass(p3, Dᵢ, F_r, th), + p3_area(p3, Dᵢ, F_r, th), + p3.ρ_i, + ), + ) +end diff --git a/src/TerminalVelocity.jl b/src/TerminalVelocity.jl deleted file mode 100644 index 7da45e153..000000000 --- a/src/TerminalVelocity.jl +++ /dev/null @@ -1,79 +0,0 @@ -""" -Terminal velocity calculations including: -- Chen 2022 parametrization for ice and rain -""" -module TerminalVelocity - -import CloudMicrophysics.Common as CO -import CloudMicrophysics.Parameters as CMP - -export velocity_chen - -""" - ϕ_ice(mass, area, ρ_i) - - - mass - mass of particle - - area - area of particle - - ρ_i - density of ice - - Returns the aspect ratio (ϕ) for a particle with given mass, area, and ice density -""" -function ϕ_ice(mass::FT, area::FT, ρ_i::FT) where {FT} - return 16 * ρ_i^2 * area^3 / (9 * π * mass^2) -end - -""" - velocity_chen(D, Chen2022, ρ_a, mass, area, ρ_i) - - - D - maximum particle dimension - - Chen2022 - a struct with terminal velocity parameters as in Chen(2022) - - ρ_a - density of air - - mass - mass of particle - - area - area of particle - - ρ_i - density of ice - - Returns the terminal velocity of ice at given particle dimension using Chen 2022 parametrizations -""" -function velocity_chen( - D::FT, - Chen2022::CMP.Chen2022VelTypeSnowIce, - ρ_a::FT, - mass::FT, - area::FT, - ρ_i::FT, -) where {FT} - if D <= Chen2022.cutoff - (ai, bi, ci) = CO.Chen2022_vel_coeffs_small(Chen2022, ρ_a) - else - (ai, bi, ci) = CO.Chen2022_vel_coeffs_large(Chen2022, ρ_a) - end - - κ = FT(-1 / 6) - - v = sum(@. sum(ai * D^bi * exp(-ci * D))) - - return ϕ_ice(mass, area, ρ_i)^κ * v -end - -""" - velocity_chen(p3, Chen2022, ρ_a, D) - - - D - maximum particle dimension - - Chen2022 - a struct with terminal velocity parameters as in Chen(2022) - - ρ_a - density of air - - Returns the terminal velocity of rain given Chen 2022 velocity parametrizations -""" -function velocity_chen( - D::FT, - Chen2022::CMP.Chen2022VelTypeRain, - ρ_a::FT, -) where {FT} - (ai, bi, ci) = CO.Chen2022_vel_coeffs_small(Chen2022, ρ_a) - - v = sum(@. sum(ai * D^bi * exp(-ci * D))) - - return v -end - -end diff --git a/test/p3_tests.jl b/test/p3_tests.jl index 15f470ef9..8645afd1e 100644 --- a/test/p3_tests.jl +++ b/test/p3_tests.jl @@ -3,7 +3,8 @@ import ClimaParams import CloudMicrophysics as CM import CloudMicrophysics.P3Scheme as P3 import CloudMicrophysics.Parameters as CMP -import CloudMicrophysics.TerminalVelocity as TV +import CloudMicrophysics.Microphysics2M as CM2 +import Thermodynamics as TD import QuadGK as QGK @@ -88,11 +89,11 @@ function test_p3_thresholds(FT) end TT.@testset "mass and area tests" begin - # values + # values ρ_r = FT(500) F_r = FT(0.5) - # get thresholds + # get thresholds D_th = P3.D_th_helper(p3) th = P3.thresholds(p3, ρ_r, F_r) (; D_gr, D_cr) = th @@ -102,13 +103,13 @@ function test_p3_thresholds(FT) D_2 = (D_th + D_gr) / 2 D_3 = (D_gr + D_cr) / 2 - # test area + # test area TT.@test P3.p3_area(p3, D_1, F_r, th) == P3.A_s(D_1) TT.@test P3.p3_area(p3, D_2, F_r, th) == P3.A_ns(p3, D_2) TT.@test P3.p3_area(p3, D_3, F_r, th) == P3.A_s(D_3) TT.@test P3.p3_area(p3, D_cr, F_r, th) == P3.A_r(p3, F_r, D_cr) - # test mass + # test mass TT.@test P3.p3_mass(p3, D_1, F_r, th) == P3.mass_s(D_1, p3.ρ_i) TT.@test P3.p3_mass(p3, D_2, F_r, th) == P3.mass_nl(p3, D_2) TT.@test P3.p3_mass(p3, D_3, F_r, th) == P3.mass_s(D_3, th.ρ_g) @@ -172,7 +173,46 @@ function test_p3_shape_solver(FT) end end -function test_velocities(FT) +function test_particle_terminal_velocities(FT) + + p3 = CMP.ParametersP3(FT) + Chen2022 = CMP.Chen2022VelType(FT) + ρ_a = FT(1.2) + + TT.@testset "Chen 2022 - Rain" begin + Ds = range(FT(1e-6), stop = FT(1e-5), length = 5) + expected = [0.002508, 0.009156, 0.01632, 0.02377, 0.03144] + for i in axes(Ds, 1) + vel = CM2.rain_particle_terminal_velocity(Ds[i], Chen2022.rain, ρ_a) + TT.@test vel >= 0 + TT.@test vel ≈ expected[i] rtol = 1e-3 + end + end + + TT.@testset "Chen 2022 - Ice" begin + F_r = FT(0.5) + ρ_r = FT(500) + th = P3.thresholds(p3, ρ_r, F_r) + # Allow for a D falling into every regime of the P3 Scheme + Ds = range(FT(0.5e-4), stop = FT(4.5e-4), length = 5) + expected = [0.08109, 0.3838, 0.5917, 0.8637, 1.148] + for i in axes(Ds, 1) + D = Ds[i] + vel = P3.ice_particle_terminal_velocity( + D, + Chen2022.snow_ice, + ρ_a, + P3.p3_mass(p3, D, F_r, th), + P3.p3_area(p3, D, F_r, th), + p3.ρ_i, + ) + TT.@test vel >= 0 + TT.@test vel ≈ expected[i] rtol = 1e-3 + end + end +end + +function test_bulk_terminal_velocities(FT) Chen2022 = CMP.Chen2022VelType(FT) p3 = CMP.ParametersP3(FT) q = FT(0.22) @@ -199,7 +239,7 @@ function test_velocities(FT) ρ_r = ρ_rs[i] F_r = F_rs[j] - calculated_vel = P3.terminal_velocity( + calculated_vel = P3.ice_terminal_velocity( p3, Chen2022.snow_ice, q, @@ -242,6 +282,187 @@ function test_velocities(FT) end end +#function test_tendencies(FT) +# +# tps = TD.Parameters.ThermodynamicsParameters(FT) +# +# p3 = CMP.ParametersP3(FT) +# Chen2022 = CMP.Chen2022VelType(FT) +# aps = CMP.AirProperties(FT) +# +# SB2006 = CMP.SB2006(FT, false) # no limiters +# pdf_r = SB2006.pdf_r +# pdf_c = SB2006.pdf_c +# +# TT.@testset "Collision Tendencies Smoke Test" begin +# N = FT(1e8) +# ρ_a = FT(1.2) +# ρ_r = FT(500) +# F_r = FT(0.5) +# T_warm = FT(300) +# T_cold = FT(200) +# +# qs = range(0.001, stop = 0.005, length = 5) +# q_const = FT(0.05) +# +# cloud_expected_warm_previous = +# [6.341e-5, 0.0002099, 0.0004258, 0.0007047, 0.001042] +# cloud_expected_cold_previous = +# [0.002196, 0.003525, 0.004687, 0.005761, 0.006781] +# cloud_expected_warm = +# [8.043e-27, 3.641e-26, 8.773e-26, 1.63e-25, 2.625e-25] +# cloud_expected_cold = +# [8.197e-33, 1.865e-32, 3.012e-32, 4.233e-32, 5.523e-32] +# rain_expected_warm = [0.000402, 0.0008436, 0.001304, 0.001777, 0.00226] +# rain_expected_cold = [0.4156, 0.4260, 0.433, 0.4383, 0.4427] +# +# for i in axes(qs, 1) +# cloud_warm = P3.ice_collisions( +# pdf_c, +# p3, +# Chen2022, +# qs[i], +# N, +# q_const, +# N, +# ρ_a, +# F_r, +# ρ_r, +# T_warm, +# ) +# cloud_cold = P3.ice_collisions( +# pdf_c, +# p3, +# Chen2022, +# qs[i], +# N, +# q_const, +# N, +# ρ_a, +# F_r, +# ρ_r, +# T_cold, +# ) +# +# TT.@test cloud_warm >= 0 +# TT.@test cloud_warm ≈ cloud_expected_warm[i] rtol = 1e-3 +# TT.@test cloud_cold >= 0 +# TT.@test cloud_cold ≈ cloud_expected_cold[i] rtol = 1e-3 +# +# rain_warm = P3.ice_collisions( +# pdf_r, +# p3, +# Chen2022, +# qs[i], +# N, +# q_const, +# N, +# ρ_a, +# F_r, +# ρ_r, +# T_warm, +# ) +# rain_cold = P3.ice_collisions( +# pdf_r, +# p3, +# Chen2022, +# qs[i], +# N, +# q_const, +# N, +# ρ_a, +# F_r, +# ρ_r, +# T_cold, +# ) +# +# TT.@test rain_warm >= 0 +# TT.@test rain_warm ≈ rain_expected_warm[i] rtol = 1e-3 +# TT.@test rain_cold >= 0 +# TT.@test rain_cold ≈ rain_expected_cold[i] rtol = 1e-3 +# end +# end +# +# TT.@testset "Melting Tendencies Smoke Test" begin +# N = FT(1e8) +# ρ_a = FT(1.2) +# ρ_r = FT(500) +# F_r = FT(0.5) +# T_freeze = FT(273.15) +# +# qs = range(0.001, stop = 0.005, length = 5) +# +# expected_melt = [0.0006982, 0.0009034, 0.001054, 0.001177, 0.001283] +# +# for i in axes(qs, 1) +# rate = P3.p3_melt( +# p3, +# Chen2022, +# aps, +# tps, +# qs[i], +# N, +# T_freeze + 2, +# ρ_a, +# F_r, +# ρ_r, +# ) +# +# TT.@test rate >= 0 +# TT.@test rate ≈ expected_melt[i] rtol = 1e-3 +# end +# end +# +# TT.@testset "Heterogeneous Freezing Smoke Test" begin +# T = FT(250) +# N = FT(1e8) +# ρ_a = FT(1.2) +# qᵥ = FT(8.1e-4) +# aero_type = CMP.Illite(FT) +# +# qs = range(0.001, stop = 0.005, length = 5) +# +# expected_freeze_q = +# [2.036e-61, 6.463e-61, 1.270e-60, 2.052e-60, 2.976e-60] +# expected_freeze_N = +# [1.414e-51, 2.244e-51, 2.941e-51, 3.562e-51, 4.134e-51] +# +# for i in axes(qs, 1) +# rate_mass = P3.p3_rain_het_freezing( +# true, +# pdf_r, +# p3, +# tps, +# qs[i], +# N, +# T, +# ρ_a, +# qᵥ, +# aero_type, +# ) +# rate_num = P3.p3_rain_het_freezing( +# false, +# pdf_r, +# p3, +# tps, +# qs[i], +# N, +# T, +# ρ_a, +# qᵥ, +# aero_type, +# ) +# +# TT.@test rate_mass >= 0 +# TT.@test rate_mass ≈ expected_freeze_q[i] rtol = 1e-3 +# +# TT.@test rate_num >= 0 +# TT.@test rate_num ≈ expected_freeze_N[i] rtol = 1e-3 +# end +# end +# +#end + function test_integrals(FT) p3 = CMP.ParametersP3(FT) Chen2022 = CMP.Chen2022VelType(FT) @@ -259,7 +480,7 @@ function test_integrals(FT) q = qs[i] # Velocity comparisons - vel_N, vel_m = P3.terminal_velocity( + vel_N, vel_m = P3.ice_terminal_velocity( p3, Chen2022.snow_ice, q, @@ -272,7 +493,7 @@ function test_integrals(FT) λ, N_0 = P3.distribution_parameter_solver(p3, q, N, ρ_r, F_r) th = P3.thresholds(p3, ρ_r, F_r) ice_bound = P3.get_ice_bound(p3, λ, tolerance) - vel(d) = TV.velocity_chen( + vel(d) = P3.ice_particle_terminal_velocity( d, Chen2022.snow_ice, ρ_a, @@ -292,7 +513,7 @@ function test_integrals(FT) TT.@test vel_N ≈ qgk_vel_N rtol = 1e-7 TT.@test vel_m ≈ qgk_vel_m rtol = 1e-7 - # Dₘ comparisons + # Dₘ comparisons D_m = P3.D_m(p3, q, N, ρ_r, F_r) f_d(d) = d * P3.p3_mass(p3, d, F_r, th) * P3.N′ice(p3, d, λ, N_0) @@ -306,6 +527,7 @@ end println("Testing Float32") test_p3_thresholds(Float32) +test_particle_terminal_velocities(Float64) #TODO - only works for Float64 now. We should switch the units inside the solver # from SI base to something more managable #test_p3_shape_solver(Float32) @@ -313,5 +535,7 @@ test_p3_thresholds(Float32) println("Testing Float64") test_p3_thresholds(Float64) test_p3_shape_solver(Float64) -test_velocities(Float64) +test_particle_terminal_velocities(Float64) +test_bulk_terminal_velocities(Float64) +#test_tendencies(Float64) test_integrals(Float64) diff --git a/test/performance_tests.jl b/test/performance_tests.jl index 156007cb0..92bd6e90d 100644 --- a/test/performance_tests.jl +++ b/test/performance_tests.jl @@ -156,7 +156,7 @@ function benchmark_test(FT) 1e5, ) bench_press( - P3.terminal_velocity, + P3.ice_terminal_velocity, (p3, ch2022.snow_ice, q_ice, N, ρ_r, F_r, ρ_air), 2e5, 3e4, diff --git a/test/runtests.jl b/test/runtests.jl index cb78fe194..cc7d85e3d 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -11,7 +11,6 @@ include("common_functions_tests.jl") include("common_types_tests.jl") include("nucleation_unit_tests.jl") include("precipitation_susceptibility_tests.jl") -include("terminal_velocity_tests.jl") include("p3_tests.jl") include("aqua.jl") include("performance_tests.jl") diff --git a/test/terminal_velocity_tests.jl b/test/terminal_velocity_tests.jl deleted file mode 100644 index ce1c1e15c..000000000 --- a/test/terminal_velocity_tests.jl +++ /dev/null @@ -1,55 +0,0 @@ -import Test as TT -import ClimaParams -import CloudMicrophysics as CM -import CloudMicrophysics.P3Scheme as P3 -import CloudMicrophysics.Parameters as CMP -import CloudMicrophysics.TerminalVelocity as TV - -@info "Terminal Velocity Tests" - -function test_chen_velocities(FT) - - p3 = CMP.ParametersP3(FT) - Chen2022 = CMP.Chen2022VelType(FT) - - ρ_a = FT(1.2) - - TT.@testset "Chen 2022 - Rain" begin - Ds = range(FT(1e-6), stop = FT(1e-5), length = 5) - expected = [0.002508, 0.009156, 0.01632, 0.02377, 0.03144] - for i in axes(Ds, 1) - vel = TV.velocity_chen(Ds[i], Chen2022.rain, ρ_a) - TT.@test vel >= 0 - TT.@test vel ≈ expected[i] rtol = 1e-3 - end - end - - TT.@testset "Chen 2022 - Ice" begin - F_r = FT(0.5) - ρ_r = FT(500) - th = P3.thresholds(p3, ρ_r, F_r) - # Allow for a D falling into every regime of the P3 Scheme - Ds = range(FT(0.5e-4), stop = FT(4.5e-4), length = 5) - expected = [0.08109, 0.3838, 0.5917, 0.8637, 1.148] - for i in axes(Ds, 1) - D = Ds[i] - vel = TV.velocity_chen( - D, - Chen2022.snow_ice, - ρ_a, - P3.p3_mass(p3, D, F_r, th), - P3.p3_area(p3, D, F_r, th), - p3.ρ_i, - ) - TT.@test vel >= 0 - TT.@test vel ≈ expected[i] rtol = 1e-3 - end - - end -end - -println("Testing Float32") -test_chen_velocities(Float32) - -println("Testing Float64") -test_chen_velocities(Float64)