Skip to content

Commit

Permalink
lambda testing with colorful plots
Browse files Browse the repository at this point in the history
  • Loading branch information
@anastasia-popova
anastasia-popova committed Feb 6, 2024
1 parent fc187c8 commit cbbb4a1
Show file tree
Hide file tree
Showing 2 changed files with 235 additions and 25 deletions.
25 changes: 18 additions & 7 deletions src/P3Scheme.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -228,7 +228,7 @@ Returns the shape parameter μ for a given λ value
Eq. 3 in Morrison and Milbrandt (2015).
"""
function DSD_μ(p3::PSP3, λ::FT) where {FT}
@assert λ > FT(0)
# @assert λ > FT(0)
return min(p3.μ_max, max(FT(0), p3.a * λ^p3.b - p3.c))
end

Expand Down Expand Up @@ -306,6 +306,7 @@ function q_gamma(

D_th = D_th_helper(p3)
λ = exp(log_λ)
# println(" log λ = ", log_λ, ", λ = ", λ)
N_0 = DSD_N₀(p3, N, λ)

return ifelse(
Expand Down Expand Up @@ -340,25 +341,35 @@ function distribution_parameter_solver(
F_r::FT,
) where {FT}

N̂ = FT(1e30)
N̂ = FT(1e20)

# Get the thresholds for different particles regimes
th = thresholds(p3, ρ_r, F_r)

# To ensure that λ is positive solve for x such that λ = exp(x)
shape_problem(x) = q/N̂ - q_gamma(p3, F_r, N/N̂, x, th)

# Find slope parameter
x =
sol =
RS.find_zero(
shape_problem,
RS.SecantMethod(FT(log(4 * 1e4)), FT(log(5 * 1e4))),
RS.CompactSolution(),
RS.RelativeSolutionTolerance(eps(FT)),
10,
)
x=sol.root
converge = sol.converged
#= x =
RS.find_zero(
shape_problem,
RS.SecantMethod(FT(log(20000)), FT(log(50000))),
RS.SecantMethod(FT(log(4 * 1e4)), FT(log(5 * 1e4))),
RS.CompactSolution(),
RS.RelativeSolutionTolerance(eps(FT)),
10,
).root

return (; λ = exp(x), N_0 = DSD_N₀(p3, N, exp(x)))
=#
return (; λ = exp(x), N_0 = DSD_N₀(p3, N, exp(x)), converged = converge)
end

end
235 changes: 217 additions & 18 deletions test/p3_LambdaTesting.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,22 +19,24 @@ function λ_diff(F_r::FT, ρ_r::FT, N::FT, λ_ex::FT, p3::PSP3)
x = log(λ_ex)
# Compute mass density based on input shape parameters
q_calc = P3.q_gamma(p3, F_r, N, x, th)
(λ_calculated, ) = P3.distribution_parameter_solver(p3, q_calc, N, ρ_r, F_r)
return abs(λ_ex - λ_calculated)
(λ_calculated, N_0, converged) = P3.distribution_parameter_solver(p3, q_calc, N, ρ_r, F_r,)
return (diff = abs(λ_ex - λ_calculated), converged = converged)
end

function find_gaps_in_solver(F_r::FT, ρ_r::FT, λ_min::FT, λ_max::FT, p3::PSP3)
function find_gaps_in_solver(F_r::FT, ρ_r::FT, λ_min::FT, λ_max::FT, p3::PSP3, numSteps)
N = FT(1e8)

prev = FT(0)

λ_range = range(λ_min, stop = λ_max, length = Int(round(abs((λ_min - λ_max)/100))))
λ_range = range(λ_min, stop = λ_max, length = numSteps)

xs = Vector{FT}()
ysFr = Vector{FT}()
ysρr = Vector{FT}()

for λ in λ_range
# println("λ = ", λ, " F_r = ", F_r, " ρ_r = ", ρ_r)

next = λ_diff(F_r, ρ_r, N, λ, p3)
p = isequal(prev, NaN)
n = isequal(next, NaN)
Expand All @@ -51,6 +53,9 @@ function find_gaps_in_solver(F_r::FT, ρ_r::FT, λ_min::FT, λ_max::FT, p3::PSP3
append!(ysρr, ρ_r)
end
prev = next

#println("done, λ_difference = ", next)
#println(" ")
end
if isequal(prev, NaN)
append!(xs, λ_max)
Expand All @@ -60,7 +65,7 @@ function find_gaps_in_solver(F_r::FT, ρ_r::FT, λ_min::FT, λ_max::FT, p3::PSP3
return (xs, ysFr, ysρr)
end

function plot_gaps(F_r_input::FT, ρ_r::FT, λ_min::FT, λ_max::FT, p3::PSP3)
function plot_gaps(F_r_input::FT, ρ_r::FT, λ_min::FT, λ_max::FT, p3::PSP3, numSteps)
fig1 = Plt.Figure()
ax1 = Plt.Axis(
fig1[1, 1],
Expand All @@ -83,24 +88,25 @@ function plot_gaps(F_r_input::FT, ρ_r::FT, λ_min::FT, λ_max::FT, p3::PSP3)
Plt.ylims!(low = 0)
Plt.ylims!(high = 1)

F_rs = range(FT(0.5), stop = FT(1-eps(FT)), length = 1000)
F_rs = range(FT(0), stop = 1-eps(FT), length = 1001) # 1-eps(FT)
for F_r in F_rs

(xs, ysFr, ysρr) = find_gaps_in_solver(F_r, ρ_r, λ_min, λ_max, p3)
(xs, ysFr, ysρr) = find_gaps_in_solver(F_r, ρ_r, λ_min, λ_max, p3, numSteps)
println("xs = ", xs)

Plt.linesegments!(ax1, xs, ysFr, color = :black)

end

# Plt.hlines!(ax1, F_r, xmin = min, xmax = max)

Plt.save("LambdaTesting1.svg", fig1)
Plt.save("LambdaTesting4.svg", fig1)


fig2 = Plt.Figure()
ax2 = Plt.Axis(
fig2[1, 1],
title = Plt.L" Working λ ",
title = Plt.L" Working λ for F_r = 0.4",
xlabel = Plt.L" λ ",
ylabel = Plt.L" ρ_r ",
#xscale = Plt.log10,
Expand All @@ -114,23 +120,216 @@ function plot_gaps(F_r_input::FT, ρ_r::FT, λ_min::FT, λ_max::FT, p3::PSP3)
#limits = ((0.02, 10.0), nothing),
)
Plt.ylims!(low = 0)
Plt.ylims!(high = 1000)s
Plt.ylims!(high = 1000)

ρ_rs = range(FT(100), stop = FT(900), length = 1000)
ρ_rs = range(FT(100), stop = FT(900), length = 500)
for ρ_r in ρ_rs
(xs, ysFr, ysρr) = find_gaps_in_solver(F_r_input, ρ_r, λ_min, λ_max, p3)
(xs, ysFr, ysρr) = find_gaps_in_solver(F_r_input, ρ_r, λ_min, λ_max, p3, numSteps)

Plt.linesegments!(ax2, xs, ysρr, color = :black)
end

Plt.save("LambdaTesting2.svg", fig2)
Plt.save("LambdaTesting3.svg", fig2)
end

function plot_all_gaps(λ_min::FT, λ_max::FT, F_r_min::FT, F_r_max::FT, ρ_r_min::FT, ρ_r_max::FT, p3::PSP3, λSteps, F_rSteps, numPlots) where {FT}
f = Plt.Figure()

λs = range(λ_min, stop = λ_max, length = λSteps)
F_rs = range(F_r_min, stop = F_r_max, length = F_rSteps)
#ρ_rs = (FT(100), FT(200), FT(300), FT(400), FT(500), FT(600), FT(700), FT(800))
#l = length(ρ_rs)
ρ_rs = range(ρ_r_min, stop = ρ_r_max, length = numPlots)

for i in 1:numPlots
ρ_r = ρ_rs[i]

j = Int(ceil(i/2))
k = Int(mod(i - 1, 2) + 1)
Plt.Axis(f[i, 1], xlabel = "λ", ylabel = "F_r", title = string("ρ_r = " , string(ρ_r)), width = 400, height = 300)

for F_r in F_rs
(xs, ysFr,) = find_gaps_in_solver(F_r, ρ_r, λ_min, λ_max, p3, λSteps)

Plt.linesegments!(f[i, 1], xs, ysFr, color = :black)
end
end
Plt.resize_to_layout!(f)
Plt.save("allFrsbig.svg", f)
end

function rel_error(λ, F_r)
ρ_r = FT(400)
N = FT(1e8)
p3 = CMP.ParametersP3(FT)

er = log(λ_diff(F_r, ρ_r, N, λ, p3)/λ)
if isequal(er, NaN)
er = Inf
end
return er
end

F_r = FT(0.8)
function plot_relerrors(λ_min::FT, λ_max::FT, p3::PSP3, λSteps, F_rSteps, numPlots) where{FT}
N = FT(1e8)

#λs = range(λ_min, stop = λ_max, length = λSteps)
#F_rs = range(FT(0), stop = 1-eps(FT), length = F_rSteps)
# ρ_rs = range(FT(10), stop = 990, length = numPlots)
ρ_r = FT(400)

#= xs = Vector{FT}()
ys = Vector{FT}()
zs = Vector{FT}()
for λ in λs
for F_r in F_rs
append!(xs, λ)
append!(ys, F_r)
er = log(λ_diff(F_r, ρ_r, N, λ, p3)/λ)
if isequal(er, NaN)
er = Inf
end
append!(zs, er)
end
end =#

f = Plt.Figure()
Plt.Axis(f[1, 1])

(λs, F_rs, E, min, max) = get_errors(λ_min, λ_max, p3, ρ_r, N, λSteps, F_rSteps)

println(min)
println(max)

Plt.heatmap!(λs, F_rs, E)
Plt.Colorbar(f[1, 2], limits = (min, max), colormap = :viridis, flipaxis = false)

# Plt.Colorbar(f[1, 2], label = "error", limits = FT.(range_val))

Plt.save("ContourAttempt1.svg", f)

end

function get_errors(λ_min::FT, λ_max::FT, p3::PSP3, ρ_r::FT, N::FT, λSteps, F_rSteps) where{FT}
λs = range(λ_min, stop = λ_max, length = λSteps)
F_rs = range(FT(0), stop = 1-eps(FT), length = F_rSteps)
E = zeros(λSteps, F_rSteps)
min = Inf
max = -Inf

for i in 1:λSteps
for j in 1:F_rSteps
λ = λs[i]
F_r = F_rs[j]

er = log(λ_diff(F_r, ρ_r, N, λ, p3)/λ)
#= if isequal(er, NaN)
er = Inf
end =#

E[i, j] = er

if er > max && er < Inf
max = er
end
if er < min && er > -Inf
min = er
end

end
end
return (λs = λs, F_rs = F_rs, E = E, min = min, max = max)
end

function buckets(λ_min::FT, λ_max::FT, F_r_min::FT, F_r_max::FT, p3::PSP3, ρ_r::FT, N::FT, λSteps, F_rSteps) where{FT}
λs = range(λ_min, stop = λ_max, length = λSteps)
F_rs = range(F_r_min, stop = F_r_max, length = F_rSteps)
E = zeros(λSteps, F_rSteps)
min_error = eps(FT)
max_error = sqrt(eps(FT))

for i in 1:λSteps
for j in 1:F_rSteps
λ = λs[i]
F_r = F_rs[j]

(error, converged) = λ_diff(F_r, ρ_r, N, λ, p3)

er = log(error/λ)

if converged
if er <= min_error
E[i, j] = FT(0)
elseif er <= max_error
E[i, j] = FT(1)
elseif er <= FT(1)
E[i, j] = FT(2)
else
E[i, j] = FT(3)
end
else
if er <= min_error
E[i, j] = FT(0.5)
elseif er <= max_error
E[i, j] = FT(1.5)
elseif er <= FT(1)
E[i, j] = FT(2.5)
else
E[i, j] = FT(3.5)
end
end
end
end
return (λs = λs, F_rs = F_rs, E = E)
end

function plot_buckets(λ_min::FT, λ_max::FT, F_r_min::FT, F_r_max::FT, ρ_r_min::FT, ρ_r_max::FT,p3::PSP3, λSteps, F_rSteps, numPlots) where{FT}
N = FT(1e8)
ρ_rs = range(ρ_r_min, stop = ρ_r_max, length = numPlots)

f = Plt.Figure()
for i in 1:numPlots
ρ_r = ρ_rs[i]

Plt.Axis(f[i, 1], xlabel = "λ", ylabel = "F_r", title = string("ρ_r = " , string(ρ_r)), width = 400, height = 300)

(λs, F_rs, E) = buckets(λ_min, λ_max, F_r_min, F_r_max, p3, ρ_r, N, λSteps, F_rSteps)

Plt.heatmap!(λs, F_rs, E)
Plt.Colorbar(f[i, 2], limits = (0, 3.5), colormap = :viridis, flipaxis = false)

# Plt.Colorbar(f[1, 2], label = "error", limits = FT.(range_val))
end

Plt.resize_to_layout!(f)
Plt.save("ContourAttempt2.svg", f)

end

F_r = FT(0.4)
ρ_r = FT(400)
λ_min = FT(15000) # anything above 400 giving errors that λ <= FT(0)
λ_max = FT(100000)
λ_min = FT(1) # anything above 400 giving errors that λ <= FT(0)
λ_max = FT(1e5)
F_r_min = FT(0.0)
F_r_max = FT(1-eps(FT))
ρ_r_min = FT(100)
ρ_r_max = FT(900)

# plot_gaps(F_r, ρ_r, λ_min, λ_max, p3, 250)

#constant_Fr_ρr(FT(0.8), FT(400), FT(0), FT(2000))

#plot_all_gaps(λ_min, λ_max, F_r_min, F_r_max, ρ_r_min, ρ_r_max, p3, 200, 200, 10)

#plot_relerrors(λ_min, λ_max, p3, 200, 200, 1)

F = FT(0.34)
ρ = FT(650)
λ = FT(35000)
N = FT(1e8)

plot_gaps(F_r, ρ_r, λ_min, λ_max, p3)
#diff = λ_diff(F, ρ, N, λ, p3)

#constant_Fr_ρr(FT(0.8), FT(400), FT(0), FT(2000))
#println(λ_diff(FT(0.5), FT(400), FT(1e8), FT(15000), p3))
plot_buckets(λ_min, λ_max, F_r_min, F_r_max, ρ_r_min, ρ_r_max, p3, 250, 250, 9)

0 comments on commit cbbb4a1

Please sign in to comment.