P3 Terminal Velocity

docs/src/API.md

@@ -58,7 +58,6 @@ Microphysics2M.conv_q_liq_to_q_rai
 ```@docs
 P3Scheme
 P3Scheme.thresholds
-P3Scheme.p3_mass
 P3Scheme.distribution_parameter_solver
 ```
 

docs/src/plots/P3SchemePlots.jl

@@ -9,6 +9,66 @@ FT = Float64
 const PSP3 = CMP.ParametersP3
 p3 = CMP.ParametersP3(FT)
 
+"""
+    mass_(p3, D, ρ, F_r)
+
+ - p3 - a struct with P3 scheme parameters
+ - D - maximum particle dimension [m]
+ - ρ - bulk ice density (ρ_i for small ice, ρ_g for graupel) [kg/m3]
+ - F_r - rime mass fraction [q_rim/q_i]
+
+Returns mass as a function of size for differen particle regimes [kg]
+"""
+# for spherical ice (small ice or completely rimed ice)
+mass_s(D::FT, ρ::FT) where {FT <: Real} = FT(π) / 6 * ρ * D^3
+# for large nonspherical ice (used for unrimed and dense types)
+mass_nl(p3::PSP3, D::FT) where {FT <: Real} = P3.α_va_si(p3) * D^p3.β_va
+# for partially rimed ice
+mass_r(p3::PSP3, D::FT, F_r::FT) where {FT <: Real} =
+    P3.α_va_si(p3) / (1 - F_r) * D^p3.β_va
+
+"""
+    p3_mass(p3, D, F_r, th)
+
+ - p3 - a struct with P3 scheme parameters
+ - D - maximum particle dimension
+ - F_r - rime mass fraction (q_rim/q_i)
+ - th - P3Scheme nonlinear solve output tuple (D_cr, D_gr, ρ_g, ρ_d)
+
+Returns mass(D) regime, used to create figures for the docs page.
+"""
+function p3_mass(
+    p3::PSP3,
+    D::FT,
+    F_r::FT,
+    th = (; D_cr = FT(0), D_gr = FT(0), ρ_g = FT(0), ρ_d = FT(0)),
+) where {FT <: Real}
+    D_th = P3.D_th_helper(p3)
+    if D_th > D
+        return mass_s(D, p3.ρ_i)          # small spherical ice
+    end
+    if F_r == 0
+        return mass_nl(p3, D)             # large nonspherical unrimed ice
+    end
+    if th.D_gr > D >= D_th
+        return mass_nl(p3, D)             # dense nonspherical ice
+    end
+    if th.D_cr > D >= th.D_gr
+        return mass_s(D, th.ρ_g)          # graupel
+    end
+    if D >= th.D_cr
+        return mass_r(p3, D, F_r)         # partially rimed ice
+    end
+
+    # TODO - would something like this be better?
+    #return ifelse(D_th_helper(p3) > D, mass_s(D, p3.ρ_i),
+    #       ifelse(F_r == 0, mass_nl(p3, D),
+    #       ifelse(th.D_gr > D >= D_th_helper(p3), mass_nl(p3, D),
+    #       ifelse(th.D_cr > D >= th.D_gr, mass_s(D, th.ρ_g),
+    #           mass_r(p3, D, F_r)))))
+
+end
+
 """
     A_(p3, D)
 
@@ -88,9 +148,9 @@ function p3_mass_plot()
     sol_8 = P3.thresholds(p3, 400.0, 0.8)
 
     #! format: off
-    fig1_a_0 = Plt.lines!(ax1_a, D_range * 1e3, [P3.p3_mass(p3, D, 0.0       ) for D in D_range], color = cl[1], linewidth = lw)
-    fig1_a_5 = Plt.lines!(ax1_a, D_range * 1e3, [P3.p3_mass(p3, D, 0.5, sol_5) for D in D_range], color = cl[2], linewidth = lw)
-    fig1_a_8 = Plt.lines!(ax1_a, D_range * 1e3, [P3.p3_mass(p3, D, 0.8, sol_8) for D in D_range], color = cl[3], linewidth = lw)
+    fig1_a_0 = Plt.lines!(ax1_a, D_range * 1e3, [p3_mass(p3, D, 0.0       ) for D in D_range], color = cl[1], linewidth = lw)
+    fig1_a_5 = Plt.lines!(ax1_a, D_range * 1e3, [p3_mass(p3, D, 0.5, sol_5) for D in D_range], color = cl[2], linewidth = lw)
+    fig1_a_8 = Plt.lines!(ax1_a, D_range * 1e3, [p3_mass(p3, D, 0.8, sol_8) for D in D_range], color = cl[3], linewidth = lw)
 
     d_tha  = Plt.vlines!(ax1_a, P3.D_th_helper(p3) * 1e3, linestyle = :dash, color = cl[4], linewidth = lw)
     d_cr_5 = Plt.vlines!(ax1_a, sol_5[1]           * 1e3, linestyle = :dot,  color = cl[2], linewidth = lw)
@@ -128,9 +188,9 @@ function p3_mass_plot()
     sol_8 = P3.thresholds(p3, 800.0, 0.95)
 
     #! format: off
-    fig1_b200 = Plt.lines!(ax1_b, D_range * 1e3, [P3.p3_mass(p3, D, 0.95, sol_2) for D in D_range], color = cl[1], linewidth = lw)
-    fig1_b400 = Plt.lines!(ax1_b, D_range * 1e3, [P3.p3_mass(p3, D, 0.95, sol_4) for D in D_range], color = cl[2], linewidth = lw)
-    fig1_b800 = Plt.lines!(ax1_b, D_range * 1e3, [P3.p3_mass(p3, D, 0.95, sol_8) for D in D_range], color = cl[3], linewidth = lw)
+    fig1_b200 = Plt.lines!(ax1_b, D_range * 1e3, [p3_mass(p3, D, 0.95, sol_2) for D in D_range], color = cl[1], linewidth = lw)
+    fig1_b400 = Plt.lines!(ax1_b, D_range * 1e3, [p3_mass(p3, D, 0.95, sol_4) for D in D_range], color = cl[2], linewidth = lw)
+    fig1_b800 = Plt.lines!(ax1_b, D_range * 1e3, [p3_mass(p3, D, 0.95, sol_8) for D in D_range], color = cl[3], linewidth = lw)
 
     d_thb    = Plt.vlines!(ax1_b, P3.D_th_helper(p3) * 1e3, linestyle = :dash, color = cl[4], linewidth = lw)
     d_cr_200 = Plt.vlines!(ax1_b, sol_2[1] * 1e3,           linestyle = :dot,  color = cl[1], linewidth = lw)

docs/src/plots/P3TerminalVelocityPlots.jl

@@ -0,0 +1,208 @@
+import CLIMAParameters
+import CloudMicrophysics as CM
+import CloudMicrophysics.P3Scheme as P3
+import CloudMicrophysics.Parameters as CMP
+import CairoMakie as Plt
+
+const PSP3 = CMP.ParametersP3
+
+FT = Float64 
+
+p3 = CMP.ParametersP3(FT)
+
+function get_values(p3::PSP3, Chen2022::CMP.Chen2022VelTypeSnowIce, q::FT, N::FT, ρ_a::FT, x_resolution::Int, y_resolution::Int) where {FT} 
+    F_rs = 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)
+    D_m = zeros(x_resolution, y_resolution)
+
+    for i in 1:x_resolution 
+        for j in 1:y_resolution 
+            F_r = F_rs[i]
+            ρ_r = ρ_rs[j]
+
+            V_m[i, j] = P3.terminal_velocity_mass(p3, Chen2022, q, N, ρ_r, F_r, ρ_a)
+            # get D_m in mm for plots 
+            D_m[i, j] = 1e3 * P3.D_m(p3, q, N, ρ_r, F_r)
+        end 
+    end
+    return (; F_rs, ρ_rs, V_m, D_m)
+end
+
+function figure_2()
+    Chen2022 = CMP.Chen2022VelType(FT)
+    # density of air in kg/m^3
+    ρ_a = FT(1.293)
+
+    f = Plt.Figure()
+    xres = 100
+    yres = 100
+    min = FT(0)
+    max = FT(10)
+
+    # small D_m
+    q_s = FT(0.0008)
+    N_s = FT(1e6)    
+
+    # medium D_m 
+    q_m = FT(0.22) 
+    N_m = FT(1e6)
+
+    # large D_m
+    q_l = FT(0.7) 
+    N_l = FT(1e6)
+
+    # get V_m and D_m
+    (F_rs, ρ_rs, V_ms, D_ms) = get_values(p3, Chen2022.snow_ice, q_s, N_s, ρ_a, xres, yres)
+    (F_rm, ρ_rm, V_mm, D_mm) = get_values(p3, Chen2022.snow_ice, q_m, N_m, ρ_a, xres, yres)
+    (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", 
+        ylabel = "ρ_r", 
+        title = "Small Dₘ", 
+        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],
+    )
+
+    Plt.contourf!(F_rs, ρ_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,)
+
+    ax2 = Plt.Axis(
+        f[1, 2], 
+        xlabel = "F_r", 
+        ylabel = "ρ_r", 
+        title = "Medium Dₘ", 
+        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],
+    )
+
+    Plt.contourf!(F_rm, ρ_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)
+
+    ax3 = Plt.Axis(
+        f[1, 3], 
+        xlabel = "F_r", 
+        ylabel = "ρ_r", 
+        title = "Large Dₘ", 
+        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],
+    )
+
+    Plt.contourf!(F_rl, ρ_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)),
+        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 
+
+    ax4 = Plt.Axis(f[3, 1], height = 350, width = 350, xlabel = "F_r", ylabel = "ρ_r", title = "Small Dₘ vs F_r and ρ_r")
+    Plt.contourf!(
+        F_rs, 
+        ρ_rs, 
+        D_ms,
+        colormap = Plt.reverse(Plt.cgrad(:Spectral))
+    )
+    Plt.Colorbar(
+        f[4, 1], 
+        limits = (minimum(d for d in D_ms if !isnan(d)), maximum(d for d in D_ms if !isnan(d))), 
+        colormap = Plt.reverse(Plt.cgrad(:Spectral)),
+        flipaxis = true,
+        vertical = false, 
+    )
+    ax5 = Plt.Axis(f[3, 2], height = 350, width = 350, xlabel = "F_r", ylabel = "ρ_r", title = "Medium Dₘ vs F_r and ρ_r")
+    Plt.contourf!(
+        F_rm, 
+        ρ_rm, 
+        D_mm,
+        colormap = Plt.reverse(Plt.cgrad(:Spectral))
+    )
+    Plt.Colorbar(
+        f[4, 2], 
+        limits = (minimum(d for d in D_mm if !isnan(d)), maximum(d for d in D_mm if !isnan(d))), 
+        colormap = Plt.reverse(Plt.cgrad(:Spectral)),
+        flipaxis = true,
+        vertical = false, 
+    )
+    ax6 = Plt.Axis(f[3, 3], height = 350, width = 350, xlabel = "F_r", ylabel = "ρ_r", title = "Large Dₘ vs F_r and ρ_r")
+    Plt.contourf!(
+        F_rl, 
+        ρ_rl, 
+        D_ml,
+        colormap = Plt.reverse(Plt.cgrad(:Spectral))
+    )
+    Plt.Colorbar(
+        f[4, 3], 
+        limits = (minimum(d for d in D_ml if !isnan(d)), maximum(d for d in D_ml if !isnan(d))), 
+        colormap = Plt.reverse(Plt.cgrad(:Spectral)),
+        flipaxis = true,
+        vertical = false, 
+    )
+
+    Plt.linkxaxes!(ax1, ax4)
+    Plt.linkxaxes!(ax4, ax5)
+    Plt.linkxaxes!(ax5, ax6)
+    Plt.linkyaxes!(ax1, ax4)
+    Plt.linkyaxes!(ax4, ax5)
+    Plt.linkyaxes!(ax5, ax6)
+
+    Plt.resize_to_layout!(f) 
+    Plt.save("MorrisonandMilbrandtFig2.svg", f)
+end
+
+println("start") 
+figure_2()
+println("done")
+#println(P3.D_th_helper(p3))
+#println(P3.thresholds(p3, FT(500), FT(0.5)))
+#println("")
+#println("small = ", P3.q_gamma(p3, FT(0.5), FT(1e7), FT(log(4.9 * 10^2)), P3.thresholds(p3, FT(500), FT(0.5))))
+
+
+import RootSolvers as RS
+ρ_r = FT(500) 
+F_r = FT(0.8)
+N = FT(1e8) 
+Dₘ = 0.006
+println("started")
+shape_problem(x) = Dₘ - P3.D_m(p3, exp(x), N, ρ_r, F_r)
+x =
+        RS.find_zero(
+            shape_problem,
+            RS.SecantMethod(log(0.01), log(0.015)),
+            RS.CompactSolution(),
+            RS.RelativeSolutionTolerance(eps(FT)),
+            5,
+        ).root
+
+println("q_solved = ", exp(x))