From 4b68eeeb555a4a657e5e8d4aaf8ccf2ffc5e2988 Mon Sep 17 00:00:00 2001 From: Haakon Ludvig Langeland Ervik <45243236+haakon-e@users.noreply.github.com> Date: Thu, 27 Feb 2025 17:44:16 -0800 Subject: [PATCH] formatter --- docs/src/plots/P3AspectRatioPlot.jl | 24 +- docs/src/plots/P3Melting.jl | 7 +- docs/src/plots/P3SchemePlots.jl | 24 +- docs/src/plots/P3SlopeParameterizations.jl | 19 +- docs/src/plots/P3TerminalVelocityPlots.jl | 266 +++++++++++++-------- docs/src/plots/P3Thresholds.jl | 10 +- src/P3.jl | 3 +- src/P3_integral_properties.jl | 12 +- src/P3_particle_properties.jl | 24 +- src/P3_size_distribution.jl | 17 +- src/P3_terminal_velocity.jl | 20 +- src/parameters/MicrophysicsP3.jl | 21 +- test/p3_tests.jl | 48 ++-- 13 files changed, 304 insertions(+), 191 deletions(-) diff --git a/docs/src/plots/P3AspectRatioPlot.jl b/docs/src/plots/P3AspectRatioPlot.jl index 5f8e56cff..24a8808c1 100644 --- a/docs/src/plots/P3AspectRatioPlot.jl +++ b/docs/src/plots/P3AspectRatioPlot.jl @@ -5,18 +5,18 @@ import CloudMicrophysics.Parameters as CMP import CloudMicrophysics.P3Scheme as P3 axis_theme = Theme( - Axis = ( - xscale = log10, - xminorticksvisible = true, xminorticks = IntervalsBetween(5), - xticks = [0.01, 0.1, 1, 10], - limits = ((0.01, 10.0), (0, 1.05)), - xgridvisible = false, ygridvisible = false, - xlabel = "D (mm)", - ), - linewidth = 3, - VLines = ( - linewidth = 1.5, - ) + Axis = ( + xscale = log10, + xminorticksvisible = true, + xminorticks = IntervalsBetween(5), + xticks = [0.01, 0.1, 1, 10], + limits = ((0.01, 10.0), (0, 1.05)), + xgridvisible = false, + ygridvisible = false, + xlabel = "D (mm)", + ), + linewidth = 3, + VLines = (linewidth = 1.5,), ) logocolors = Makie.Colors.JULIA_LOGO_COLORS diff --git a/docs/src/plots/P3Melting.jl b/docs/src/plots/P3Melting.jl index 08c63994b..d40b7ee4f 100644 --- a/docs/src/plots/P3Melting.jl +++ b/docs/src/plots/P3Melting.jl @@ -10,8 +10,11 @@ FT = Float64 # parameters tps = TD.Parameters.ThermodynamicsParameters(FT) params = CMP.ParametersP3( - CP.create_toml_dict(FT; override_file=Dict("Heymsfield_mu_coeff1" => Dict( "value" => 3.0))); - slope_law=:constant + CP.create_toml_dict( + FT; + override_file = Dict("Heymsfield_mu_coeff1" => Dict("value" => 3.0)), + ); + slope_law = :constant, ) vel = CMP.Chen2022VelType(FT) aps = CMP.AirProperties(FT) diff --git a/docs/src/plots/P3SchemePlots.jl b/docs/src/plots/P3SchemePlots.jl index 88e3ef606..0b2317e0a 100644 --- a/docs/src/plots/P3SchemePlots.jl +++ b/docs/src/plots/P3SchemePlots.jl @@ -5,17 +5,17 @@ import CloudMicrophysics.Parameters as CMP import CloudMicrophysics.P3Scheme as P3 axis_theme = Theme( - Axis = ( - xscale = log10, - xminorticksvisible = true, xminorticks = IntervalsBetween(5), - xticks = [0.01, 0.1, 1, 10], - limits = ((0.01, 10.0), nothing), - xgridvisible = false, ygridvisible = false, - ), - linewidth = 3, - VLines = ( - linewidth = 1.5, - ) + Axis = ( + xscale = log10, + xminorticksvisible = true, + xminorticks = IntervalsBetween(5), + xticks = [0.01, 0.1, 1, 10], + limits = ((0.01, 10.0), nothing), + xgridvisible = false, + ygridvisible = false, + ), + linewidth = 3, + VLines = (linewidth = 1.5,), ) logocolors = Makie.Colors.JULIA_LOGO_COLORS @@ -96,4 +96,4 @@ end fig = with_theme(p3_relations_plot, axis_theme) save("P3Scheme_relations.svg", fig) -fig \ No newline at end of file +fig diff --git a/docs/src/plots/P3SlopeParameterizations.jl b/docs/src/plots/P3SlopeParameterizations.jl index 46520658a..38f4defcd 100644 --- a/docs/src/plots/P3SlopeParameterizations.jl +++ b/docs/src/plots/P3SlopeParameterizations.jl @@ -10,30 +10,29 @@ FT = Float64 slope_power_law = CMP.ParametersP3(FT).slope # Constant parameterization -slope_constant = CMP.ParametersP3( - CP.create_toml_dict(FT; override_file=Dict("Heymsfield_mu_coeff1" => Dict( "value" => 3.0))); - slope_law=:constant -).slope +override_file = Dict("Heymsfield_mu_coeff1" => Dict("value" => 3.0)) +params = CMP.ParametersP3(CP.create_toml_dict(FT; override_file)) +slope_constant = params.slope -logλs = @. log(10.0 ^ FT(2:0.01:6)) + +logλs = @. log(10.0^FT(2:0.01:6)) function make_slope_plot(slope_law, title) - fig = Figure(size=(400, 300), figure_padding = 20) + fig = Figure(size = (400, 300), figure_padding = 20) - ax = Axis(fig[1,1]; title, xscale = log10, xlabel = "λ", ylabel = "μ") + ax = Axis(fig[1, 1]; title, xscale = log10, xlabel = "λ", ylabel = "μ") lines!(ax, exp.(logλs), P3.get_μ.(slope_law, logλs)) return fig end -fig = with_theme(theme_minimal()) do +fig = with_theme(theme_minimal()) do make_slope_plot(slope_power_law, "μ as a function of λ (power law)") end save("P3SlopeParameterizations_power_law.svg", fig) -fig = with_theme(theme_minimal()) do +fig = with_theme(theme_minimal()) do make_slope_plot(slope_constant, "μ as a constant") end save("P3SlopeParameterizations_constant.svg", fig) - diff --git a/docs/src/plots/P3TerminalVelocityPlots.jl b/docs/src/plots/P3TerminalVelocityPlots.jl index 9e09f7ca1..a28857d9d 100644 --- a/docs/src/plots/P3TerminalVelocityPlots.jl +++ b/docs/src/plots/P3TerminalVelocityPlots.jl @@ -8,87 +8,87 @@ const PSP3 = CMP.ParametersP3 FT = Float64 -params = CMP.ParametersP3( - CP.create_toml_dict(FT; override_file=Dict("Heymsfield_mu_coeff1" => Dict( "value" => 3.0))); - slope_law=:constant -) +override_file = Dict("Heymsfield_mu_coeff1" => Dict("value" => 3.0)) +toml_dict = CP.create_toml_dict(FT; override_file) +params = CMP.ParametersP3(toml_dict; slope_law = :constant) function get_values( - params::CMP.ParametersP3, - Chen2022::CMP.Chen2022VelType, - L::FT, - N::FT, - ρ_a::FT, - x_resolution::Int, - y_resolution::Int, + params::CMP.ParametersP3, + Chen2022::CMP.Chen2022VelType, + L::FT, + N::FT, + ρ_a::FT, + x_resolution::Int, + y_resolution::Int, ) where {FT} - F_rims = range(FT(0), stop = FT(1 - eps(FT)), length = x_resolution) - ρ_rs = range(FT(25), stop = FT(975), length = y_resolution) - - V_m = zeros(x_resolution, y_resolution) - V_m_ϕ = zeros(x_resolution, y_resolution) - D_m = zeros(x_resolution, y_resolution) - D_m_regimes = zeros(x_resolution, y_resolution) - ϕᵢ = zeros(x_resolution, y_resolution) - - for i in 1:x_resolution - for j in 1:y_resolution - F_rim = F_rims[i] - ρ_r = ρ_rs[j] - use_aspect_ratio = true - do_not_use_aspect = false - state = P3.get_state(params; F_rim, ρ_r) - dist = P3.get_distribution_parameters(state; L, N) - - try - V_m[i, j] = P3.ice_terminal_velocity(dist, Chen2022, ρ_a, do_not_use_aspect)[2] - catch e - @show F_rim, ρ_r, L, N, ρ_a, do_not_use_aspect - rethrow() - end - V_m_ϕ[i, j] = P3.ice_terminal_velocity(dist, Chen2022, ρ_a, use_aspect_ratio)[2] - - D_m[i, j] = P3.D_m(dist) - D_m_regimes[i, j] = D_m[i, j] - ϕᵢ[i, j] = P3.ϕᵢ(state, D_m[i, j]) - - # plot the regimes - if state.D_th > D_m[i, j] - # small spherical ice - D_m_regimes[i, j] = 0 - elseif state.F_rim == 0 - # large nonspherical unrimed ice - D_m_regimes[i, j] = 0 - elseif state.D_gr > D_m[i, j] >= state.D_th - # dense nonspherical ice - D_m_regimes[i, j] = 1 - elseif state.D_cr > D_m[i, j] >= state.D_gr - # graupel - D_m_regimes[i, j] = 2 - else #elseif D >= state.D_cr - # partially rimed ice - D_m_regimes[i, j] = 3 - end - end - end - D_m *= 1e3 - return (; F_rims, ρ_rs, D_m_regimes, D_m, ϕᵢ, V_m, V_m_ϕ) + F_rims = range(FT(0), stop = FT(1 - eps(FT)), length = x_resolution) + ρ_rs = range(FT(25), stop = FT(975), length = y_resolution) + + V_m = zeros(x_resolution, y_resolution) + V_m_ϕ = zeros(x_resolution, y_resolution) + D_m = zeros(x_resolution, y_resolution) + D_m_regimes = zeros(x_resolution, y_resolution) + ϕᵢ = zeros(x_resolution, y_resolution) + + for i in 1:x_resolution + for j in 1:y_resolution + F_rim = F_rims[i] + ρ_r = ρ_rs[j] + use_aspect_ratio = true + do_not_use_aspect = false + state = P3.get_state(params; F_rim, ρ_r) + dist = P3.get_distribution_parameters(state; L, N) + + V_m[i, j] = P3.ice_terminal_velocity( + dist, + Chen2022, + ρ_a, + do_not_use_aspect, + )[2] + V_m_ϕ[i, j] = + P3.ice_terminal_velocity(dist, Chen2022, ρ_a, use_aspect_ratio)[2] + + D_m[i, j] = P3.D_m(dist) + D_m_regimes[i, j] = D_m[i, j] + ϕᵢ[i, j] = P3.ϕᵢ(state, D_m[i, j]) + + # plot the regimes + if state.D_th > D_m[i, j] + # small spherical ice + D_m_regimes[i, j] = 0 + elseif state.F_rim == 0 + # large nonspherical unrimed ice + D_m_regimes[i, j] = 0 + elseif state.D_gr > D_m[i, j] >= state.D_th + # dense nonspherical ice + D_m_regimes[i, j] = 1 + elseif state.D_cr > D_m[i, j] >= state.D_gr + # graupel + D_m_regimes[i, j] = 2 + else #elseif D >= state.D_cr + # partially rimed ice + D_m_regimes[i, j] = 3 + end + end + end + D_m *= 1e3 + return (; F_rims, ρ_rs, D_m_regimes, D_m, ϕᵢ, V_m, V_m_ϕ) end theme = Theme( Axis = (; width = 350, - height = 350, - limits = ((0, 1.0), (25, 975)), - xticks = [0.0, 0.2, 0.4, 0.6, 0.8, 1.0], - yticks = [200, 400, 600, 800], - ), + height = 350, + limits = ((0, 1.0), (25, 975)), + xticks = [0.0, 0.2, 0.4, 0.6, 0.8, 1.0], + yticks = [200, 400, 600, 800], + ), Contour = (; color = :black, - labels = true, - levels = 3, - linewidth = 1.5, - labelsize = 18, + labels = true, + levels = 3, + linewidth = 1.5, + labelsize = 18, ), ) @@ -112,71 +112,125 @@ function figure_2() xres = 100 yres = 100 - (F_rims, ρ_rs, D_m_regimes_s, D_m_s, ϕᵢ_s, V_m_s, V_m_ϕ_s) = get_values(params, Chen2022, L_s, N_s, ρ_a, xres, yres) - (F_rimm, ρ_rm, D_m_regimes_m, D_m_m, ϕᵢ_m, V_m_m, V_m_ϕ_m) = get_values(params, Chen2022, L_m, N_m, ρ_a, xres, yres) - (F_riml, ρ_rl, D_m_regimes_l, D_m_l, ϕᵢ_l, V_m_l, V_m_ϕ_l) = get_values(params, Chen2022, L_l, N_l, ρ_a, xres, yres) + (F_rims, ρ_rs, D_m_regimes_s, D_m_s, ϕᵢ_s, V_m_s, V_m_ϕ_s) = + get_values(params, Chen2022, L_s, N_s, ρ_a, xres, yres) + (F_rimm, ρ_rm, D_m_regimes_m, D_m_m, ϕᵢ_m, V_m_m, V_m_ϕ_m) = + get_values(params, Chen2022, L_m, N_m, ρ_a, xres, yres) + (F_riml, ρ_rl, D_m_regimes_l, D_m_l, ϕᵢ_l, V_m_l, V_m_ϕ_l) = + get_values(params, Chen2022, L_l, N_l, ρ_a, xres, yres) ### PLOT ### - fig = Figure() + fig = Figure() - # Plot velocities as in Fig 2 in Morrison and Milbrandt 2015 + # Plot velocities as in Fig 2 in Morrison and Milbrandt 2015 colormap = cgrad(:PuBuGn_3, 3, categorical = true) regime_contour_kwargs = (; levels = 3, colormap) row = 1 - ax1 = Axis(fig[row,1]; title = "Particle regimes with small Dₘ") - hm = contourf!(ax1, F_rims, ρ_rs, D_m_regimes_s; regime_contour_kwargs...) + ax1 = Axis(fig[row, 1]; title = "Particle regimes with small Dₘ") + hm = contourf!(ax1, F_rims, ρ_rs, D_m_regimes_s; regime_contour_kwargs...) - ax2 = Axis(fig[row,2]; title = "Particle regimes with medium Dₘ") - hm = contourf!(ax2, F_rimm, ρ_rm, D_m_regimes_m; regime_contour_kwargs...) + ax2 = Axis(fig[row, 2]; title = "Particle regimes with medium Dₘ") + hm = contourf!(ax2, F_rimm, ρ_rm, D_m_regimes_m; regime_contour_kwargs...) - ax3 = Axis(fig[row,3]; title = "Particle regimes with large Dₘ") - hm = contourf!(ax3, F_riml, ρ_rl, D_m_regimes_l; regime_contour_kwargs...) + ax3 = Axis(fig[row, 3]; title = "Particle regimes with large Dₘ") + hm = contourf!(ax3, F_riml, ρ_rl, D_m_regimes_l; regime_contour_kwargs...) map(1:3) do col - Colorbar(fig[row,col]; - colormap, ticks = ([1, 2, 3], ["dense\n nonspherical ice", "graupel", "partially\n rimed ice"]), - vertical = false, width = Relative(0.95), height = 10, halign = 0.5, valign = 0.02, tellheight = false, - colorrange = (0.5, 3.5), ticklabelpad = 0, + ticks = ( + [1, 2, 3], + ["dense\n nonspherical ice", "graupel", "partially\n rimed ice"], + ) + Colorbar( + fig[row, col]; + colormap, + ticks, + vertical = false, + width = Relative(0.95), + height = 10, + halign = 0.5, + valign = 0.02, + tellheight = false, + colorrange = (0.5, 3.5), + ticklabelpad = 0, ) end - function make_plots(row, col, F_rim, ρ_r; cfvals, cvals = nothing, title = "") + function make_plots( + row, + col, + F_rim, + ρ_r; + cfvals, + cvals = nothing, + title = "", + ) gp = fig[row, col] ax = Axis(gp; title) - row3_opts = row == 3 ? (; ticklabelcolor = :white, ticks = 0:0.25:1, leftspinecolor = :white, rightspinecolor = :white, bottomspinecolor = :white, topspinecolor = :white) : (;) + row3_opts = + row == 3 ? + (; + ticklabelcolor = :white, + ticks = 0:0.25:1, + leftspinecolor = :white, + rightspinecolor = :white, + bottomspinecolor = :white, + topspinecolor = :white, + ) : (;) hm = contourf!(ax, F_rim, ρ_r, cfvals) - Colorbar(gp, hm; halign = 0.05, valign = 0.05, height = Relative(0.60), width = 10, tellwidth = false, ticklabelpad = 0, row3_opts...) + Colorbar( + gp, + hm; + halign = 0.05, + valign = 0.05, + height = Relative(0.60), + width = 10, + tellwidth = false, + ticklabelpad = 0, + row3_opts..., + ) !isnothing(cvals) && contour!(ax, F_rim, ρ_r, cvals) end row += 1 - make_plots(row, 1, F_rims, ρ_rs; cfvals = D_m_s, title = "Dₘ (mm) (L = 8e-4 kgm⁻³, N = 1e6 m⁻³)") - make_plots(row, 2, F_rimm, ρ_rm; cfvals = D_m_m, title = "Dₘ (mm) (L = 0.22 kgm⁻³, N = 1e6 m⁻³)") - make_plots(row, 3, F_riml, ρ_rl; cfvals = D_m_l, title = "Dₘ (mm) (L = 0.7 kgm⁻³, N = 1e6 m⁻³)") + title = "Dₘ (mm) (L = 8e-4 kgm⁻³, N = 1e6 m⁻³)" + make_plots(row, 1, F_rims, ρ_rs; cfvals = D_m_s, title) + title = "Dₘ (mm) (L = 0.22 kgm⁻³, N = 1e6 m⁻³)" + make_plots(row, 2, F_rimm, ρ_rm; cfvals = D_m_m, title) + title = "Dₘ (mm) (L = 0.7 kgm⁻³, N = 1e6 m⁻³)" + make_plots(row, 3, F_riml, ρ_rl; cfvals = D_m_l, title) row += 1 - make_plots(row, 1, F_rims, ρ_rs; cfvals = ϕᵢ_s, cvals = D_m_s, title = "ϕᵢ with small Dₘ") - make_plots(row, 2, F_rimm, ρ_rm; cfvals = ϕᵢ_m, cvals = D_m_m, title = "ϕᵢ with medium Dₘ") - make_plots(row, 3, F_riml, ρ_rl; cfvals = ϕᵢ_l, cvals = D_m_l, title = "ϕᵢ with large Dₘ") - + title = "ϕᵢ with small Dₘ" + make_plots(row, 1, F_rims, ρ_rs; cfvals = ϕᵢ_s, cvals = D_m_s, title) + title = "ϕᵢ with medium Dₘ" + make_plots(row, 2, F_rimm, ρ_rm; cfvals = ϕᵢ_m, cvals = D_m_m, title) + title = "ϕᵢ with large Dₘ" + make_plots(row, 3, F_riml, ρ_rl; cfvals = ϕᵢ_l, cvals = D_m_l, title) + row += 1 - make_plots(row, 1, F_rims, ρ_rs; cfvals = V_m_s, cvals = D_m_s, title = "Vₘ (ϕᵢ = 1) with small Dₘ") - make_plots(row, 2, F_rimm, ρ_rm; cfvals = V_m_m, cvals = D_m_m, title = "Vₘ (ϕᵢ = 1) with medium Dₘ") - make_plots(row, 3, F_riml, ρ_rl; cfvals = V_m_l, cvals = D_m_l, title = "Vₘ (ϕᵢ = 1) with large Dₘ") + title = "Vₘ (ϕᵢ = 1) with small Dₘ" + make_plots(row, 1, F_rims, ρ_rs; cfvals = V_m_s, cvals = D_m_s, title) + title = "Vₘ (ϕᵢ = 1) with medium Dₘ", + make_plots(row, 2, F_rimm, ρ_rm; cfvals = V_m_m, cvals = D_m_m, title) + title = "Vₘ (ϕᵢ = 1) with large Dₘ" + make_plots(row, 3, F_riml, ρ_rl; cfvals = V_m_l, cvals = D_m_l, title) row += 1 - make_plots(row, 1, F_rims, ρ_rs; cfvals = V_m_ϕ_s, cvals = D_m_s, title = "Vₘ (using ϕᵢ) with small Dₘ") - make_plots(row, 2, F_rimm, ρ_rm; cfvals = V_m_ϕ_m, cvals = D_m_m, title = "Vₘ (using ϕᵢ) with medium Dₘ") - make_plots(row, 3, F_riml, ρ_rl; cfvals = V_m_ϕ_l, cvals = D_m_l, title = "Vₘ (using ϕᵢ) with large Dₘ") + title = "Vₘ (using ϕᵢ) with small Dₘ" + make_plots(row, 1, F_rims, ρ_rs; cfvals = V_m_ϕ_s, cvals = D_m_s, title) + title = "Vₘ (using ϕᵢ) with medium Dₘ" + make_plots(row, 2, F_rimm, ρ_rm; cfvals = V_m_ϕ_m, cvals = D_m_m, title) + title = "Vₘ (using ϕᵢ) with large Dₘ" + make_plots(row, 3, F_riml, ρ_rl; cfvals = V_m_ϕ_l, cvals = D_m_l, title) axs = filter(ax -> ax isa Axis, fig.content) - linkaxes!(axs...) + linkaxes!(axs...) - resize_to_layout!(fig) - save("MorrisonandMilbrandtFig2.svg", fig) - fig + resize_to_layout!(fig) + save("MorrisonandMilbrandtFig2.svg", fig) + fig end #! format: on diff --git a/docs/src/plots/P3Thresholds.jl b/docs/src/plots/P3Thresholds.jl index 4a28868e5..a52393131 100644 --- a/docs/src/plots/P3Thresholds.jl +++ b/docs/src/plots/P3Thresholds.jl @@ -15,7 +15,7 @@ get_state2(F_rim, ρ_r) = P3.get_state(params; F_rim, ρ_r) states = get_state2.(F_rims, ρ_rs') -fig = Figure(size=(1200, 900), figure_padding = 20) +fig = Figure(size = (1200, 900), figure_padding = 20) D_th = round(states[1].D_th * 1e3, digits = 4) # mm # Make a plot for each threshold @@ -26,23 +26,23 @@ for (i, key) in enumerate(thresh_keys) else "" end - ax = Axis(fig[i,1], xlabel = "F_rim", ylabel = "ρ_r", title = title) + ax = Axis(fig[i, 1], xlabel = "F_rim", ylabel = "ρ_r", title = title) threshold = getproperty.(states, key) * 1e3 hm = heatmap!(ax, F_rims, ρ_rs, threshold) - Colorbar(fig[i,2], hm, label = "mm") + Colorbar(fig[i, 2], hm, label = "mm") end # Make a plot for each density ρ_keys = (:ρ_g, :ρ_d) for (i, key) in enumerate(ρ_keys) - ax = Axis(fig[i,3], xlabel = "F_rim", ylabel = "ρ_r", title = string(key)) + ax = Axis(fig[i, 3], xlabel = "F_rim", ylabel = "ρ_r", title = string(key)) density = if key == :ρ_g getproperty.(states, key) else # :ρ_d P3.get_ρ_d_exact.(params.mass, F_rims, ρ_rs') end hm = heatmap!(ax, F_rims, ρ_rs, density) - Colorbar(fig[i,4], hm, label = "kg/m³") + Colorbar(fig[i, 4], hm, label = "kg/m³") end save("P3Thresholds.svg", fig) diff --git a/src/P3.jl b/src/P3.jl index 0f7a0b7ab..abf9ee72a 100644 --- a/src/P3.jl +++ b/src/P3.jl @@ -24,7 +24,8 @@ import Thermodynamics as TD import Thermodynamics.Parameters as TDP const PSP3 = CMP.ParametersP3 -using CloudMicrophysics.Parameters: MassPowerLaw, AreaPowerLaw, SlopeLaw, SlopePowerLaw, SlopeConstant +using CloudMicrophysics.Parameters: + MassPowerLaw, AreaPowerLaw, SlopeLaw, SlopePowerLaw, SlopeConstant include("P3_particle_properties.jl") include("P3_size_distribution.jl") diff --git a/src/P3_integral_properties.jl b/src/P3_integral_properties.jl index 5db5f3a8d..f937444a5 100644 --- a/src/P3_integral_properties.jl +++ b/src/P3_integral_properties.jl @@ -115,12 +115,18 @@ function D_m(dist::P3Distribution{FT}) where {FT} else G_dense_nonspherical = G(D_th, D_gr, α_va, β_va + 1, μ, log_λ) G_graupel = G(D_gr, D_cr, ρ_g * π / 6, 4, μ, log_λ) - G_partially_rimed = G(D_cr, ∞, α_va / (1 - F_rim), β_va + 1, μ, log_λ) - logsumexp((G_small_spherical, G_dense_nonspherical, G_graupel, G_partially_rimed)) + G_partially_rimed = + G(D_cr, ∞, α_va / (1 - F_rim), β_va + 1, μ, log_λ) + logsumexp(( + G_small_spherical, + G_dense_nonspherical, + G_graupel, + G_partially_rimed, + )) end # compute `F_liq`-weighted average and normalize by `L` return exp(G_summed) * exp(log_N₀) / L # TODO: Implement `F_liq`-weighted average # return weighted_average(state.F_liq, exp(G_liqfrac), exp(G_summed)) * exp(log_N₀) / L -end \ No newline at end of file +end diff --git a/src/P3_particle_properties.jl b/src/P3_particle_properties.jl index b5f7afb8c..65f252d32 100644 --- a/src/P3_particle_properties.jl +++ b/src/P3_particle_properties.jl @@ -138,10 +138,10 @@ F_rim, ρ_r = FT(0.5), FT(916.7) ``` """ function get_ρ_d_exact((; β_va)::MassPowerLaw, F_rim, ρ_r) - k = (1 - F_rim) ^ (-1/(3 - β_va)) - num = ρ_r * F_rim - den = (β_va - 2) * (k - 1) / ( (1 - F_rim) * k - 1 ) - (1 - F_rim) - return num / den + k = (1 - F_rim)^(-1 / (3 - β_va)) + num = ρ_r * F_rim + den = (β_va - 2) * (k - 1) / ((1 - F_rim) * k - 1) - (1 - F_rim) + return num / den end """ @@ -188,7 +188,8 @@ where for the different thresholds, `ρ` is: - `params`: [`CMP.MassPowerLaw`](@ref) parameters - `ρ`: density [kg/m³] """ -get_threshold((; α_va, β_va)::MassPowerLaw, ρ) = (6α_va / (π * ρ))^(1 / (3 - β_va)) +get_threshold((; α_va, β_va)::MassPowerLaw, ρ) = + (6α_va / (π * ρ))^(1 / (3 - β_va)) """ get_D_th(mass::MassPowerLaw, ρ_i) @@ -215,7 +216,8 @@ Return the size of equal mass for graupel and partially rimed ice [meters] See Eq. 14 in [MorrisonMilbrandt2015](@cite). """ -get_D_cr(mass::MassPowerLaw, ρ_g, F_rim) = get_threshold(mass, ρ_g * (1 - F_rim)) +get_D_cr(mass::MassPowerLaw, ρ_g, F_rim) = + get_threshold(mass, ρ_g * (1 - F_rim)) """ weighted_average(f_a, a, b) @@ -354,9 +356,11 @@ function ∂ice_mass_∂D(state::P3State, D) (; mass, ρ_i) = get_parameters(state) ∂mass_spherical_∂D(ρ, D) = ρ * D^2 * π / 2 - ∂mass_nonspherical_∂D((; α_va, β_va)::MassPowerLaw, D) = β_va * α_va * D^(β_va - 1) - ∂mass_rimed_∂D((; α_va, β_va)::MassPowerLaw, D, F_rim) = β_va * α_va * D^(β_va - 1) / (1 - F_rim) - + ∂mass_nonspherical_∂D((; α_va, β_va)::MassPowerLaw, D) = + β_va * α_va * D^(β_va - 1) + ∂mass_rimed_∂D((; α_va, β_va)::MassPowerLaw, D, F_rim) = + β_va * α_va * D^(β_va - 1) / (1 - F_rim) + return if D < D_th ∂mass_spherical_∂D(ρ_i, D) # small spherical ice elseif isunrimed(state) @@ -379,7 +383,7 @@ Calculate the cross-sectional area of a spherical particle. # Arguments - `D`: maximum particle dimension [m] """ -area_spherical(D) = D^2 * π / 4 +area_spherical(D) = D^2 * π / 4 """ area_nonspherical(area, D) diff --git a/src/P3_size_distribution.jl b/src/P3_size_distribution.jl index 12bb3f7db..613ef3f23 100644 --- a/src/P3_size_distribution.jl +++ b/src/P3_size_distribution.jl @@ -117,7 +117,8 @@ get_μ(dist::P3Distribution) = get_μ(get_state(dist), dist.log_λ) get_μ(state::P3State, log_λ) = get_μ(get_parameters(state), log_λ) get_μ(params::PSP3, log_λ) = get_μ(params.slope, log_λ) -get_μ((; a, b, c, μ_max)::SlopePowerLaw, log_λ) = clamp(a * exp(log_λ)^b - c, 0, μ_max) +get_μ((; a, b, c, μ_max)::SlopePowerLaw, log_λ) = + clamp(a * exp(log_λ)^b - c, 0, μ_max) get_μ((; μ)::SlopeConstant, _) = μ """ @@ -295,7 +296,11 @@ Find all solutions for `log_λ` given the `state` ([`P3State`](@ref)), `L`, and [`SlopePowerLaw`](@ref) parameterization, which can have multiple solutions for `log_λ` for a given `log_L` and `log_N`. """ -function get_distribution_parameters_all_solutions(state::P3State{FT}; L, N) where {FT} +function get_distribution_parameters_all_solutions( + state::P3State{FT}; + L, + N, +) where {FT} # Find bounds by evaluating function incrementally, then apply root finding with bounds above and below zero-point target_log_LdN = log(L) - log(N) @@ -305,9 +310,9 @@ function get_distribution_parameters_all_solutions(state::P3State{FT}; L, N) whe λs = 10.0 .^ (2.0:Δλ:6.0) λ_bnds = Tuple[] # Loop over λs and find where shape_problem changes sign - for i = 1:length(λs)-1 - if shape_problem(log(λs[i])) * shape_problem(log(λs[i+1])) < 0 - push!(λ_bnds, (λs[i], λs[i+1])) + for i in 1:(length(λs) - 1) + if shape_problem(log(λs[i])) * shape_problem(log(λs[i + 1])) < 0 + push!(λ_bnds, (λs[i], λs[i + 1])) end end @@ -331,4 +336,4 @@ function Base.show(io::IO, p3s::P3Distribution{FT}) where {FT} println(io, "├── state: is $rimed") println(io, "├── log_λ = $(p3s.log_λ) [log(1/m)]") println(io, "└── log_N₀ = $(p3s.log_N₀) [log(1/m^3)]") -end \ No newline at end of file +end diff --git a/src/P3_terminal_velocity.jl b/src/P3_terminal_velocity.jl index d15fe8c73..8df4c3186 100644 --- a/src/P3_terminal_velocity.jl +++ b/src/P3_terminal_velocity.jl @@ -56,7 +56,8 @@ function p3_particle_terminal_velocity( use_aspect_ratio = true, ) where {FT} # TODO: Refactor to be a parameterization for use with `F_liq` - v_i = ice_particle_terminal_velocity(state, D, Chen2022, ρₐ, use_aspect_ratio) + v_i = + ice_particle_terminal_velocity(state, D, Chen2022, ρₐ, use_aspect_ratio) # v_r = CM2.rain_particle_terminal_velocity(D, Chen2022.rain, ρₐ) return v_i # return weighted_average(state.F_liq, v_r, v_i) @@ -145,7 +146,13 @@ See Eq. C10 of [MorrisonMilbrandt2015](@cite) and use the Chen 2022 terminal vel - `v_n`: The number weighted fall speed - `v_m`: The mass weighted fall speed """ -function ice_terminal_velocity(dist::P3Distribution{FT}, Chen2022::CMP.Chen2022VelType, ρₐ::FT, use_aspect_ratio = true; accurate = false) where {FT} +function ice_terminal_velocity( + dist::P3Distribution{FT}, + Chen2022::CMP.Chen2022VelType, + ρₐ::FT, + use_aspect_ratio = true; + accurate = false, +) where {FT} state = get_state(dist) (; L, N) = dist if N < eps(FT) || L < eps(FT) @@ -161,8 +168,13 @@ function ice_terminal_velocity(dist::P3Distribution{FT}, Chen2022::CMP.Chen2022V # ∫N(D) v(D) dD v_n = ∫fdD(state; accurate) do D - N′ice(dist, D) * - p3_particle_terminal_velocity(state, D, Chen2022, ρₐ, use_aspect_ratio) + N′ice(dist, D) * p3_particle_terminal_velocity( + state, + D, + Chen2022, + ρₐ, + use_aspect_ratio, + ) end return (v_n / N, v_m / L) diff --git a/src/parameters/MicrophysicsP3.jl b/src/parameters/MicrophysicsP3.jl index a69c37454..9aea052f9 100644 --- a/src/parameters/MicrophysicsP3.jl +++ b/src/parameters/MicrophysicsP3.jl @@ -1,5 +1,6 @@ export ParametersP3 -export MassPowerLaw, AreaPowerLaw, SlopePowerLaw, SlopeConstant, VentilationSB2005 +export MassPowerLaw, + AreaPowerLaw, SlopePowerLaw, SlopeConstant, VentilationSB2005 ### ----------------------------- ### ### --- SUB-PARAMETERIZATIONS --- ### @@ -42,7 +43,10 @@ Base.@kwdef struct MassPowerLaw{FT} <: ParametersType{FT} β_va::FT end function MassPowerLaw(td::CP.AbstractTOMLDict) - name_map = (; :BF1995_mass_coeff_alpha => :α_va, :BF1995_mass_exponent_beta => :β_va) + name_map = (; + :BF1995_mass_coeff_alpha => :α_va, + :BF1995_mass_exponent_beta => :β_va, + ) (; β_va) = p = CP.get_parameter_values(td, name_map, "CloudMicrophysics") α_va = p.α_va * 10^(6 * β_va - 3) FT = CP.float_type(td) @@ -79,7 +83,8 @@ Base.@kwdef struct AreaPowerLaw{FT} <: ParametersType{FT} σ::FT end function AreaPowerLaw(td::CP.AbstractTOMLDict) - name_map = (; :M1996_area_coeff_gamma => :γ, :M1996_area_exponent_sigma => :σ) + name_map = + (; :M1996_area_coeff_gamma => :γ, :M1996_area_exponent_sigma => :σ) params = CP.get_parameter_values(td, name_map, "CloudMicrophysics") FT = CP.float_type(td) return AreaPowerLaw{FT}(; params...) @@ -175,7 +180,7 @@ Base.@kwdef struct SlopeConstant{FT} <: SlopeLaw{FT} μ::FT end function SlopeConstant(td::CP.AbstractTOMLDict) - name_map = (; :Heymsfield_mu_coeff1 => :μ,) + name_map = (; :Heymsfield_mu_coeff1 => :μ) params = CP.get_parameter_values(td, name_map, "CloudMicrophysics") FT = CP.float_type(td) return SlopeConstant{FT}(; params...) @@ -288,7 +293,8 @@ ParametersP3 └── T_freeze = 273.15 [K] ``` """ -ParametersP3(::Type{FT}; kw...) where {FT<:AbstractFloat} = ParametersP3(CP.create_toml_dict(FT); kw...) +ParametersP3(::Type{FT}; kw...) where {FT <: AbstractFloat} = + ParametersP3(CP.create_toml_dict(FT); kw...) """ ParametersP3(td::CP.AbstractTOMLDict; slope_law = :powerlaw) @@ -380,7 +386,10 @@ function Base.show( if typeof(value) <: Number # Simple value - print directly with unit unit = _get_parameter_unit(field) - println(io, "$(indent)$(prefix_char)── $(field) = $(value) [$(unit)]") + println( + io, + "$(indent)$(prefix_char)── $(field) = $(value) [$(unit)]", + ) else # Nested struct - recursively show with proper context nested_io = IOContext( diff --git a/test/p3_tests.jl b/test/p3_tests.jl index d45bd9719..aa62320fe 100644 --- a/test/p3_tests.jl +++ b/test/p3_tests.jl @@ -30,12 +30,24 @@ function test_p3_state_creation(FT) # Test thresholds for rimed state @test length(P3.threshold_tuple(state_rimed)) == 3 @test P3.threshold_tuple(state_rimed) == - (state_rimed.D_th, state_rimed.D_gr, state_rimed.D_cr) + (state_rimed.D_th, state_rimed.D_gr, state_rimed.D_cr) # Test parameter boundary validation - @test_throws AssertionError P3.get_state(params; F_rim = FT(-0.1), ρ_r = FT(400)) - @test_throws AssertionError P3.get_state(params; F_rim = FT(1), ρ_r = FT(400)) - @test_throws AssertionError P3.get_state(params; F_rim = FT(0.5), ρ_r = FT(-400)) + @test_throws AssertionError P3.get_state( + params; + F_rim = FT(-0.1), + ρ_r = FT(400), + ) + @test_throws AssertionError P3.get_state( + params; + F_rim = FT(1), + ρ_r = FT(400), + ) + @test_throws AssertionError P3.get_state( + params; + F_rim = FT(0.5), + ρ_r = FT(-400), + ) end end @@ -81,7 +93,11 @@ function test_thresholds_solver(FT) state = P3.get_state(params; F_rim, ρ_r) ρ_d = P3.get_ρ_d_exact(params.mass, F_rim, ρ_r) - @test get_ρ_d(params.mass; D_cr = state.D_cr, D_gr = state.D_gr) ≈ ρ_d + @test get_ρ_d( + params.mass; + D_cr = state.D_cr, + D_gr = state.D_gr, + ) ≈ ρ_d @test P3.get_ρ_g(ρ_r, F_rim, ρ_d) ≈ state.ρ_g @test P3.get_D_cr(params.mass, state.ρ_g, F_rim) ≈ state.D_cr @test P3.get_D_gr(params.mass, state.ρ_g) ≈ state.D_gr @@ -130,13 +146,15 @@ function test_thresholds_solver(FT) @test P3.ice_area(state, D_1) == P3.area_spherical(D_1) @test P3.ice_area(state, D_2) == P3.area_nonspherical(params.area, D_2) @test P3.ice_area(state, D_3) == P3.area_spherical(D_3) - @test P3.ice_area(state, D_cr) == P3.area_rimed(params.area, F_rim, D_cr) + @test P3.ice_area(state, D_cr) == + P3.area_rimed(params.area, F_rim, D_cr) # test mass @test P3.ice_mass(state, D_1) == P3.mass_spherical(params.ρ_i, D_1) @test P3.ice_mass(state, D_2) == P3.mass_nonspherical(params.mass, D_2) @test P3.ice_mass(state, D_3) == P3.mass_spherical(state.ρ_g, D_3) - @test P3.ice_mass(state, D_cr) == P3.mass_rimed(params.mass, D_cr, F_rim) + @test P3.ice_mass(state, D_cr) == + P3.mass_rimed(params.mass, D_cr, F_rim) # test density @test P3.ice_density(state, D_1) ≈ params.ρ_i @@ -199,11 +217,12 @@ function test_shape_solver(FT) if L_calc < FT(1) # Solve for shape parameters - (; log_λ, log_N₀) = P3.get_distribution_parameters( - state; - L = L_calc, - N, - ) + (; log_λ, log_N₀) = + P3.get_distribution_parameters( + state; + L = L_calc, + N, + ) # Compare solved values with the input expected values @test log_λ ≈ logλ_ex rtol = ep @@ -336,12 +355,13 @@ function test_bulk_terminal_velocities(FT) ref_v_m = [7.790036674136336, 5.799151695989574] ref_v_m_ϕ = [7.790036674136336, 5.799151695989574] - for k = 1:length(F_rims) + for k in 1:length(F_rims) F_rim = F_rims[k] state = P3.get_state(params; F_rim, ρ_r) dist = P3.get_distribution_parameters(state; L, N) vel_n, vel_m = P3.ice_terminal_velocity(dist, Chen2022, ρ_a, false) - vel_n_ϕ, vel_m_ϕ = P3.ice_terminal_velocity(dist, Chen2022, ρ_a, true) + vel_n_ϕ, vel_m_ϕ = + P3.ice_terminal_velocity(dist, Chen2022, ρ_a, true) # number weighted @test vel_n > 0