Skip to content

Commit

Permalink
clean up
Browse files Browse the repository at this point in the history
  • Loading branch information
@anastasia-popova
anastasia-popova committed Apr 10, 2024
1 parent 6e275d7 commit 001b718
Show file tree
Hide file tree
Showing 4 changed files with 80 additions and 109 deletions.
74 changes: 28 additions & 46 deletions docs/src/plots/P3TerminalVelocityPlots.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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",
Expand All @@ -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],
Expand All @@ -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],
Expand All @@ -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],
Expand Down Expand Up @@ -270,4 +252,4 @@ function figure_2()
end

# Terminal Velocity figure
figure_2()
figure_2()
81 changes: 27 additions & 54 deletions src/P3Scheme.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand All @@ -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] + λ,
)
Expand All @@ -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] + λ,
)
Expand All @@ -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] + λ,
)
Expand All @@ -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] + λ,
)
Expand All @@ -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

Expand All @@ -550,57 +525,64 @@ 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,
th.D_gr,
α_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
if !large && th.D_cr > cutoff
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
Expand All @@ -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

Expand All @@ -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 ->
Expand All @@ -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
Expand Down
18 changes: 9 additions & 9 deletions test/p3_tests.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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)
Expand Down
16 changes: 16 additions & 0 deletions test/performance_tests.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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(
Expand Down

0 comments on commit 001b718

Please sign in to comment.