From 001b7182d2ac6ccf13923616d64cb5b92a136ff4 Mon Sep 17 00:00:00 2001 From: Anastasia Popova Date: Tue, 9 Apr 2024 13:52:14 -0400 Subject: [PATCH] clean up --- docs/src/plots/P3TerminalVelocityPlots.jl | 74 ++++++++------------- src/P3Scheme.jl | 81 ++++++++--------------- test/p3_tests.jl | 18 ++--- test/performance_tests.jl | 16 +++++ 4 files changed, 80 insertions(+), 109 deletions(-) diff --git a/docs/src/plots/P3TerminalVelocityPlots.jl b/docs/src/plots/P3TerminalVelocityPlots.jl index f99dbe352..71fcabf15 100644 --- a/docs/src/plots/P3TerminalVelocityPlots.jl +++ b/docs/src/plots/P3TerminalVelocityPlots.jl @@ -70,34 +70,6 @@ function figure_2() (F_rl, ρ_rl, V_ml, D_ml) = get_values(p3, Chen2022.snow_ice, q_l, N_l, ρ_a, xres, yres) - println( - "small q = ", - q_s, - ": min: ", - minimum(d for d in D_ms if !isnan(d)), - " max: ", - maximum(d for d in D_ms if !isnan(d)), - ) - println( - "medium q = ", - q_m, - ": min: ", - minimum(d for d in D_mm if !isnan(d)), - " max: ", - maximum(d for d in D_mm if !isnan(d)), - ) - println( - "large q = ", - q_l, - ": min: ", - minimum(d for d in D_ml if !isnan(d)), - " max: ", - maximum(d for d in D_ml if !isnan(d)), - ) - - points = [0, 0.1, 0.2, 0.5, 1, 2, 3, 4, 5, 7.5, 10] - labels = string.(points) - ax1 = Plt.Axis( f[1, 1], xlabel = "F_r", @@ -115,11 +87,18 @@ function figure_2() ρ_rs, V_ms, colormap = Plt.reverse(Plt.cgrad(:Spectral)), - colorrange = (min, max), - levels = points, - extendhigh = :darkred, ) Plt.contour!(F_rs, ρ_rs, D_ms, color = :black, labels = true, levels = 3) + Plt.Colorbar( + f[2, 1], + limits = ( + minimum(v for v in V_ms if !isnan(v)), + maximum(v for v in V_ms if !isnan(v)), + ), + colormap = Plt.reverse(Plt.cgrad(:Spectral)), + flipaxis = true, + vertical = false, + ) ax2 = Plt.Axis( f[1, 2], @@ -138,11 +117,18 @@ function figure_2() ρ_rm, V_mm, colormap = Plt.reverse(Plt.cgrad(:Spectral)), - colorrange = (min, max), - levels = points, - extendhigh = :darkred, ) Plt.contour!(F_rm, ρ_rm, D_mm, color = :black, labels = true, levels = 3) + Plt.Colorbar( + f[2, 2], + limits = ( + minimum(v for v in V_mm if !isnan(v)), + maximum(v for v in V_mm if !isnan(v)), + ), + colormap = Plt.reverse(Plt.cgrad(:Spectral)), + flipaxis = true, + vertical = false, + ) ax3 = Plt.Axis( f[1, 3], @@ -161,29 +147,25 @@ function figure_2() ρ_rl, V_ml, colormap = Plt.reverse(Plt.cgrad(:Spectral)), - colorrange = (min, max), - levels = points, - extendhigh = :darkred, ) Plt.contour!(F_rl, ρ_rl, D_ml, color = :black, labels = true, levels = 3) - Plt.Colorbar( - f[2, 2], - limits = (min, max), - colormap = Plt.reverse( - Plt.cgrad(:Spectral, length(points), categorical = true), + f[2, 3], + limits = ( + minimum(v for v in V_ml if !isnan(v)), + maximum(v for v in V_ml if !isnan(v)), ), + colormap = Plt.reverse(Plt.cgrad(:Spectral)), flipaxis = true, vertical = false, - ticks = (points, labels), - highclip = :darkred, ) Plt.linkxaxes!(ax1, ax2) Plt.linkxaxes!(ax2, ax3) Plt.linkyaxes!(ax1, ax2) Plt.linkyaxes!(ax2, ax3) - # Attempt to plot D_m for clearer information than the contour plots and debugging + + # Plot D_m as second row of comparisons ax4 = Plt.Axis( f[3, 1], @@ -270,4 +252,4 @@ function figure_2() end # Terminal Velocity figure -figure_2() \ No newline at end of file +figure_2() diff --git a/src/P3Scheme.jl b/src/P3Scheme.jl index 5a5f3a48e..6bca4ea95 100644 --- a/src/P3Scheme.jl +++ b/src/P3Scheme.jl @@ -251,31 +251,6 @@ function integrate(a::FT, b::FT, c1::FT, c2::FT, c3::FT) where {FT} end end -""" - get_coeffs(p3, th) - - - p3 - a struct with P3 scheme parameters - - th - thresholds tuple as returned by thresholds() - - F_r - rime mass fraction [q_rim/q_i] - -Returns the coefficients for m(D), a(D), and the respective powers of D - Where the indices are as follows: - 1 - small, spherical ice - 2 - large, unrimed ice - 3 - dense, nonspherical ice - 4 - graupel - 5 - partially rimed ice - 6 - second half of partially rimed ice (only for a) -""" -function get_coeffs(p3::PSP3, th, F_r::FT) where {FT} - α_va = α_va_si(p3) - m = [π / 6 * p3.ρ_i, α_va, α_va, π / 6 * th.ρ_g, α_va / (1 - F_r)] - m_power = [FT(3), p3.β_va, p3.β_va, 3, p3.β_va] - a = [π / 4, p3.γ, p3.γ, π / 4, F_r * π / 4, (1 - F_r) * p3.γ] - a_power = [FT(2), p3.σ, p3.σ, FT(2), FT(2), p3.σ] - return (; m, m_power, a, a_power) -end - """ q_(p3, ρ, F_r, λ, μ, D_min, D_max) @@ -473,6 +448,7 @@ function terminal_velocity_mass( (λ, N_0) = distribution_parameter_solver(p3, q, N, ρ_r, F_r) μ = DSD_μ(p3, λ) + # TO DO: Change when each value used depending on type of particle κ = FT(-1 / 6) #FT(1/3) # Redefine α_va to be in si units @@ -482,11 +458,11 @@ function terminal_velocity_mass( for i in 1:2 if F_r == 0 # Velocity coefficients for small particles - (ai, bi, ci) = CO.Chen2022_vel_coeffs_small(Chen2022, ρ_a) + (ai, bi, ci) = CO.Chen2022_vel_coeffs_small(Chen2022, ρ_a) v += integrate( FT(0), D_th, - π / 6 * p3.ρ_i * ai[i] * N_0, + FT(π) / 6 * p3.ρ_i * ai[i] * N_0, bi[i] + μ + 3, ci[i] + λ, ) @@ -496,7 +472,7 @@ function terminal_velocity_mass( α_va * ai[i] * N_0 * - (16 * p3.ρ_i^2 * p3.γ^3 / (9 * π * α_va^2))^κ, + (16 * p3.ρ_i^2 * p3.γ^3 / (9 * FT(π) * α_va^2))^κ, bi[i] + μ + p3.β_va + κ * (3 * p3.σ - 2 * p3.β_va), ci[i] + λ, ) @@ -509,7 +485,7 @@ function terminal_velocity_mass( α_va * ai[i] * N_0 * - (16 * p3.ρ_i^2 * p3.γ^3 / (9 * π * α_va^2))^κ, + (16 * p3.ρ_i^2 * p3.γ^3 / (9 * FT(π) * α_va^2))^κ, bi[i] + μ + p3.β_va + κ * (3 * p3.σ - 2 * p3.β_va), ci[i] + λ, ) @@ -521,7 +497,7 @@ function terminal_velocity_mass( v += integrate( FT(0), D_th, - π / 6 * p3.ρ_i * ai[i] * N_0, + FT(π) / 6 * p3.ρ_i * ai[i] * N_0, bi[i] + μ + 3, ci[i] + λ, ) @@ -534,13 +510,12 @@ function terminal_velocity_mass( α_va * ai[i] * N_0 * - (16 * p3.ρ_i^2 * p3.γ^3 / (9 * π * α_va^2))^κ, + (16 * p3.ρ_i^2 * p3.γ^3 / (9 * FT(π) * α_va^2))^κ, bi[i] + μ + p3.β_va + κ * (3 * p3.σ - 2 * p3.β_va), ci[i] + λ, ) - #= println("D_th to cutoff") =# - # large particles + # Switch to large particles (ai, bi, ci) = CO.Chen2022_vel_coeffs_large(Chen2022, ρ_a) large = true @@ -550,11 +525,10 @@ function terminal_velocity_mass( α_va * ai[i] * N_0 * - (16 * p3.ρ_i^2 * p3.γ^3 / (9 * π * α_va^2))^κ, + (16 * p3.ρ_i^2 * p3.γ^3 / (9 * FT(π) * α_va^2))^κ, bi[i] + μ + p3.β_va + κ * (3 * p3.σ - 2 * p3.β_va), ci[i] + λ, ) - #= println("large, cutoff to D_gr") =# else v += integrate( D_th, @@ -562,11 +536,10 @@ function terminal_velocity_mass( α_va * ai[i] * N_0 * - (16 * p3.ρ_i^2 * p3.γ^3 / (9 * π * α_va^2))^κ, + (16 * p3.ρ_i^2 * p3.γ^3 / (9 * FT(π) * α_va^2))^κ, bi[i] + μ + p3.β_va + κ * (3 * p3.σ - 2 * p3.β_va), ci[i] + λ, ) - #= println("D_th to D_gr") =# end # D_gr to D_cr @@ -574,33 +547,42 @@ function terminal_velocity_mass( v += integrate( th.D_gr, cutoff, - π / 6 * th.ρ_g * ai[i] * N_0 * (p3.ρ_i / th.ρ_g)^(2 * κ), + FT(π) / 6 * + th.ρ_g * + ai[i] * + N_0 * + (p3.ρ_i / th.ρ_g)^(2 * κ), bi[i] + μ + 3, ci[i] + λ, ) - #= println("D_gr to cutoff") =# - # large particles + # Switch to large particles (ai, bi, ci) = CO.Chen2022_vel_coeffs_large(Chen2022, ρ_a) large = true v += integrate( cutoff, th.D_cr, - π / 6 * th.ρ_g * ai[i] * N_0 * (p3.ρ_i / th.ρ_g)^(2 * κ), + FT(π) / 6 * + th.ρ_g * + ai[i] * + N_0 * + (p3.ρ_i / th.ρ_g)^(2 * κ), bi[i] + μ + 3, ci[i] + λ, ) - #= println("large, cutoff to D_cr") =# else v += integrate( th.D_gr, th.D_cr, - π / 6 * th.ρ_g * ai[i] * N_0 * (p3.ρ_i / th.ρ_g)^(2 * κ), + FT(π) / 6 * + th.ρ_g * + ai[i] * + N_0 * + (p3.ρ_i / th.ρ_g)^(2 * κ), bi[i] + μ + 3, ci[i] + λ, ) - #= printlln("D_gr to D_cr") =# end # D_cr to Infinity @@ -624,12 +606,8 @@ function terminal_velocity_mass( cutoff, ) v += I - #= println("D_cr to cutoff") =# - - # approximating sigma as 2 (for closed form integration) (this overestimates A LOT) - #v += integrate(th.D_cr, cutoff, 1/(1 - F_r) * α_va * ai[i] * N_0 * (16 * p3.ρ_i ^ 2 * (F_r * π / 4 + (1 - F_r) * p3.γ)^3 / (9 * π * (1/(1 - F_r) * α_va) ^ 2)) ^ (κ), bi[i] + μ + p3.β_va + κ*(6 - 2 * p3.β_va), ci[i] + λ) - # large particles + # Switch to large particles (ai, bi, ci) = CO.Chen2022_vel_coeffs_large(Chen2022, ρ_a) large = true @@ -652,8 +630,6 @@ function terminal_velocity_mass( Inf, ) v += I - #= println("large, cutoff to Infinity") =# - #v += integrate(cutoff, Inf, 1/(1 - F_r) * α_va * ai[i] * N_0 * (16 * p3.ρ_i ^ 2 * (F_r * π / 4 + (1 - F_r) * p3.γ)^3 / (9 * π * (1/(1 - F_r) * α_va) ^ 2)) ^ (κ), bi[i] + μ + p3.β_va + κ*(6 - 2 * p3.β_va), ci[i] + λ) else (I, e) = QGK.quadgk( D -> @@ -674,9 +650,6 @@ function terminal_velocity_mass( Inf, ) v += I - #= println("D_cr to Infinity") =# - # approximating sigma as 2 (for closed form integration) (this overestimates A LOT) - #v += integrate(th.D_cr, Inf, 1/(1 - F_r) * α_va * ai[i] * N_0 * (16 * p3.ρ_i ^ 2 * (F_r * π / 4 + (1 - F_r) * p3.γ)^3 / (9 * π * (1/(1 - F_r) * α_va) ^ 2)) ^ (κ), bi[i] + μ + p3.β_va + κ*(6 - 2 * p3.β_va), ci[i] + λ) end end end diff --git a/test/p3_tests.jl b/test/p3_tests.jl index ac3a09614..dcc0f29a0 100644 --- a/test/p3_tests.jl +++ b/test/p3_tests.jl @@ -142,14 +142,14 @@ function test_velocities(FT) N = FT(1e6) ρ_a = FT(1.2) ρ_rs = [FT(200), FT(400), FT(600), FT(800)] - F_rs = [FT(0.2), FT(0.4), FT(0.6), FT(0.8)] + F_rs = [FT(0), FT(0.2), FT(0.4), FT(0.6), FT(0.8)] TT.@testset "Mass-weighted terminal velocities" begin paper_vals = [ - [1.5, 1.5, 1.5, 1.5], - [1.5, 2.5, 2.5, 2.5], - [2.5, 2.5, 2.5, 2.5], - [2.5, 3.5, 3.5, 3.5], + [1.5, 1.5, 1.5, 1.5, 1.5], + [1.5, 1.5, 2.5, 2.5, 2.5], + [1.5, 2.5, 2.5, 2.5, 2.5], + [1.5, 2.5, 3.5, 3.5, 3.5], ] for i in 1:length(ρ_rs) for j in 1:length(F_rs) @@ -175,10 +175,10 @@ function test_velocities(FT) TT.@testset "Mass-weighted mean diameters" begin paper_vals = [ - [5, 5, 5, 5], - [4.5, 4.5, 4.5, 4.5], - [3.5, 3.5, 3.5, 3.5], - [3.5, 2.5, 2.5, 2.5], + [5, 5, 5, 5, 5], + [4.5, 4.5, 4.5, 4.5, 4.5], + [3.5, 3.5, 3.5, 3.5, 3.5], + [3.5, 3.5, 2.5, 2.5, 2.5], ] for i in 1:length(ρ_rs) for j in 1:length(F_rs) diff --git a/test/performance_tests.jl b/test/performance_tests.jl index cbe34b7e3..25c13af00 100644 --- a/test/performance_tests.jl +++ b/test/performance_tests.jl @@ -86,6 +86,7 @@ function benchmark_test(FT) ρ_r = FT(400.0) F_r = FT(0.95) + N = FT(1e8) T_air_2 = FT(250) T_air_cold = FT(230) @@ -121,6 +122,21 @@ function benchmark_test(FT) # P3 scheme bench_press(P3.thresholds, (p3, ρ_r, F_r), 12e6, 2048, 80) + if FT == Float64 + bench_press( + P3.distribution_parameter_solver, + (p3, q_ice, N, ρ_r, F_r), + 1e5, + ) + bench_press( + P3.terminal_velocity_mass, + (p3, ch2022.snow_ice, q_ice, N, ρ_r, F_r, ρ_air), + 2e6, + 4e4, + 2e3, + ) + bench_press(P3.D_m, (p3, q_ice, N, ρ_r, F_r), 1e5) + end # aerosol activation bench_press(