Skip to content

feat(P3): implement additional rates + some fixes #681

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
haakon-e wants to merge 4 commits into
base: main
Choose a base branch
from
Open

feat(P3): implement additional rates + some fixes #681

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
a217ef1
remove dt limiting from p3 melt, nuc
haakon-e Feb 21, 2026
39f6fc5
fix(P3): Manually unroll LogExpFunctions.logsumexp
haakon-e Mar 19, 2026
b752add
feat(P3): add more P3 ice tendencies
haakon-e Mar 11, 2026
f7f4661
wip: number adjustment tendency
haakon-e Mar 22, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
26 changes: 12 additions & 14 deletions docs/src/plots/P3ImmersionFreezing.jl
View file
Open in desktop
Comment thread
haakon-e marked this conversation as resolved.
Original file line number Diff line number Diff line change
Expand Up @@ -30,9 +30,6 @@ p = FT(800 * 1e2)
dd = CMP.DesertDust(FT)
il = CMP.Illite(FT)

# model time step (for limiting)
dt = FT(1)

# plot data
RH_range = range(0.8, stop = 1.2, length = 1000)
T1 = FT(273.15 - 15)
Expand All @@ -41,22 +38,23 @@ T2 = FT(273.15 - 35)
#! format: off

# limiters to not nucleate more mass and number than we have in liquid phase
max_dLdt_T1 = [qₗ* p2ρ(T1, RH) / dt for RH in RH_range]
max_dLdt_T2 = [qₗ* p2ρ(T2, RH) / dt for RH in RH_range]
max_dNdt = [Nₗ / dt for RH in RH_range]
# (implicitly dt = 1)
max_dLdt_T1 = @. qₗ * p2ρ(T1, RH_range)
max_dLdt_T2 = @. qₗ * p2ρ(T2, RH_range)
max_dNdt = fill(Nₗ, size(RH_range))

dLdt_dd_T1 = [P3.het_ice_nucleation(dd, tps, qₗ, Nₗ, RH, T1, p2ρ(T1, RH), dt).dLdt for RH in RH_range]
dNdt_dd_T1 = [P3.het_ice_nucleation(dd, tps, qₗ, Nₗ, RH, T1, p2ρ(T1, RH), dt).dNdt for RH in RH_range]
dLdt_dd_T1 = [P3.het_ice_nucleation(dd, tps, qₗ, Nₗ, RH, T1, p2ρ(T1, RH)).dLdt for RH in RH_range]
dNdt_dd_T1 = [P3.het_ice_nucleation(dd, tps, qₗ, Nₗ, RH, T1, p2ρ(T1, RH)).dNdt for RH in RH_range]

dLdt_il_T1 = [P3.het_ice_nucleation(il, tps, qₗ, Nₗ, RH, T1, p2ρ(T1, RH), dt).dLdt for RH in RH_range]
dNdt_il_T1 = [P3.het_ice_nucleation(il, tps, qₗ, Nₗ, RH, T1, p2ρ(T1, RH), dt).dNdt for RH in RH_range]
dLdt_il_T1 = [P3.het_ice_nucleation(il, tps, qₗ, Nₗ, RH, T1, p2ρ(T1, RH)).dLdt for RH in RH_range]
dNdt_il_T1 = [P3.het_ice_nucleation(il, tps, qₗ, Nₗ, RH, T1, p2ρ(T1, RH)).dNdt for RH in RH_range]


dLdt_dd_T2 = [P3.het_ice_nucleation(dd, tps, qₗ, Nₗ, RH, T2, p2ρ(T2, RH), dt).dLdt for RH in RH_range]
dNdt_dd_T2 = [P3.het_ice_nucleation(dd, tps, qₗ, Nₗ, RH, T2, p2ρ(T2, RH), dt).dNdt for RH in RH_range]
dLdt_dd_T2 = [P3.het_ice_nucleation(dd, tps, qₗ, Nₗ, RH, T2, p2ρ(T2, RH)).dLdt for RH in RH_range]
dNdt_dd_T2 = [P3.het_ice_nucleation(dd, tps, qₗ, Nₗ, RH, T2, p2ρ(T2, RH)).dNdt for RH in RH_range]

dLdt_il_T2 = [P3.het_ice_nucleation(il, tps, qₗ, Nₗ, RH, T2, p2ρ(T2, RH), dt).dLdt for RH in RH_range]
dNdt_il_T2 = [P3.het_ice_nucleation(il, tps, qₗ, Nₗ, RH, T2, p2ρ(T2, RH), dt).dNdt for RH in RH_range]
dLdt_il_T2 = [P3.het_ice_nucleation(il, tps, qₗ, Nₗ, RH, T2, p2ρ(T2, RH)).dLdt for RH in RH_range]
dNdt_il_T2 = [P3.het_ice_nucleation(il, tps, qₗ, Nₗ, RH, T2, p2ρ(T2, RH)).dNdt for RH in RH_range]

# plotting
fig = PL.Figure(size = (1500, 500), fontsize=22, linewidth=3)
Expand Down
116 changes: 58 additions & 58 deletions docs/src/plots/P3Melting.jl
View file
Open in desktop
Comment thread
haakon-e marked this conversation as resolved.
Original file line number Diff line number Diff line change
Expand Up @@ -23,73 +23,73 @@ Nᵢ = FT(2e5)
ΔT_range = range(1e-4, stop = 0.025, length = 1000)

# limiters to not melt more mass and number than we have
max_dLdt = [Lᵢ / dt for ΔT in ΔT_range]
max_dNdt = [Nᵢ / dt for ΔT in ΔT_range]
max_dLdt = Lᵢ / dt
max_dNdt = Nᵢ / dt

#! format: off

ρₐ1 = FT(1.2)
Fᵣ1 = FT(0.8)
ρᵣ1 = FT(800)
state = P3.get_state(params; F_rim = Fᵣ1, ρ_rim = ρᵣ1, L_ice = Lᵢ, N_ice = Nᵢ)
ρₐ = FT(1.2)
Fᵣ = FT(0.8)
ρᵣ = FT(800)
label1 = "ρₐ=$ρₐ kg/m³, Fᵣ=$Fᵣ, ρᵣ=$(Int(ρᵣ)) kg/m³"
state = P3.get_state(params; F_rim = Fᵣ, ρ_rim = ρᵣ, L_ice = Lᵢ, N_ice = Nᵢ)
logλ = P3.get_distribution_logλ(state)
dLdt1 = [P3.ice_melt(vel, aps, tps, params.T_freeze .+ ΔT, ρₐ1, dt, state, logλ).dLdt for ΔT in ΔT_range]
dNdt1 = [P3.ice_melt(vel, aps, tps, params.T_freeze .+ ΔT, ρₐ1, dt, state, logλ).dNdt for ΔT in ΔT_range]
melt1 = P3.ice_melt.(vel, aps, tps, params.T_freeze .+ ΔT_range, ρₐ, state, logλ)
dLdt1 = getfield.(melt1, :dLdt)
dNdt1 = getfield.(melt1, :dNdt)

Fᵣ2 = FT(0.2)
state = P3.get_state(params; F_rim = Fᵣ2, ρ_rim = ρᵣ1, L_ice = Lᵢ, N_ice = Nᵢ)
Fᵣ = FT(0.2)
label2 = "ρₐ=$ρₐ kg/m³, Fᵣ=$Fᵣ, ρᵣ=$(Int(ρᵣ)) kg/m³"
state = P3.get_state(params; F_rim = Fᵣ, ρ_rim = ρᵣ, L_ice = Lᵢ, N_ice = Nᵢ)
logλ = P3.get_distribution_logλ(state)
dLdt2 = [P3.ice_melt(vel, aps, tps, params.T_freeze .+ ΔT, ρₐ1, dt, state, logλ).dLdt for ΔT in ΔT_range]
dNdt2 = [P3.ice_melt(vel, aps, tps, params.T_freeze .+ ΔT, ρₐ1, dt, state, logλ).dNdt for ΔT in ΔT_range]
melt2 = P3.ice_melt.(vel, aps, tps, params.T_freeze .+ ΔT_range, ρₐ, state, logλ)
dLdt2 = getfield.(melt2, :dLdt)
dNdt2 = getfield.(melt2, :dNdt)

ρᵣ2 = FT(200)
state = P3.get_state(params; F_rim = Fᵣ2, ρ_rim = ρᵣ2, L_ice = Lᵢ, N_ice = Nᵢ)
ρᵣ = FT(200)
label3 = "ρₐ=$ρₐ kg/m³, Fᵣ=$Fᵣ, ρᵣ=$(Int(ρᵣ)) kg/m³"
state = P3.get_state(params; F_rim = Fᵣ, ρ_rim = ρᵣ, L_ice = Lᵢ, N_ice = Nᵢ)
logλ = P3.get_distribution_logλ(state)
dLdt3 = [P3.ice_melt(vel, aps, tps, params.T_freeze .+ ΔT, ρₐ1, dt, state, logλ).dLdt for ΔT in ΔT_range]
dNdt3 = [P3.ice_melt(vel, aps, tps, params.T_freeze .+ ΔT, ρₐ1, dt, state, logλ).dNdt for ΔT in ΔT_range]
melt3 = P3.ice_melt.(vel, aps, tps, params.T_freeze .+ ΔT_range, ρₐ, state, logλ)
dLdt3 = getfield.(melt3, :dLdt)
dNdt3 = getfield.(melt3, :dNdt)

ρₐ2 = FT(0.5)
state = P3.get_state(params; F_rim = Fᵣ2, ρ_rim = ρᵣ2, L_ice = Lᵢ, N_ice = Nᵢ)
ρₐ = FT(0.5)
label4 = "ρₐ=$ρₐ kg/m³, Fᵣ=$Fᵣ, ρᵣ=$(Int(ρᵣ)) kg/m³"
state = P3.get_state(params; F_rim = Fᵣ, ρ_rim = ρᵣ, L_ice = Lᵢ, N_ice = Nᵢ)
logλ = P3.get_distribution_logλ(state)
dLdt4 = [P3.ice_melt(vel, aps, tps, params.T_freeze .+ ΔT, ρₐ2, dt, state, logλ).dLdt for ΔT in ΔT_range]
dNdt4 = [P3.ice_melt(vel, aps, tps, params.T_freeze .+ ΔT, ρₐ2, dt, state, logλ).dNdt for ΔT in ΔT_range]
melt4 = P3.ice_melt.(vel, aps, tps, params.T_freeze .+ ΔT_range, ρₐ, state, logλ)
dLdt4 = getfield.(melt4, :dLdt)
dNdt4 = getfield.(melt4, :dNdt)

# plotting
fig = Makie.Figure(size = (1500, 500), fontsize=22, linewidth=3)

ax1 = Makie.Axis(fig[1, 1]; yscale = log10)
ax2 = Makie.Axis(fig[1, 2]; yscale = log10)

ax1.xlabel = "T [C]"
ax1.ylabel = "ice mass melting rate [g/m3/s]"
ax2.xlabel = "T [C]"
ax2.ylabel = "ice number melting rate [1/cm3/s]"

l_max_dLdt = Makie.lines!(ax1, ΔT_range, max_dLdt * 1e3, color = :thistle)
l_max_dNdt = Makie.lines!(ax2, ΔT_range, max_dNdt * 1e-6, color = :thistle)

l_dLdt1 = Makie.lines!(ax1, ΔT_range, dLdt1 * 1e3, color = :skyblue)
l_dNdt1 = Makie.lines!(ax2, ΔT_range, dNdt1 * 1e-6, color = :skyblue)

l_dLdt2 = Makie.lines!(ax1, ΔT_range, dLdt2 * 1e3, color = :blue3)
l_dNdt2 = Makie.lines!(ax2, ΔT_range, dNdt2 * 1e-6, color = :blue3)

l_dLdt3 = Makie.lines!(ax1, ΔT_range, dLdt3 * 1e3, color = :orchid)
l_dNdt3 = Makie.lines!(ax2, ΔT_range, dNdt3 * 1e-6, color = :orchid)

l_dLdt4 = Makie.lines!(ax1, ΔT_range, dLdt4 * 1e3, color = :purple)
l_dNdt4 = Makie.lines!(ax2, ΔT_range, dNdt4 * 1e-6, color = :purple)

Makie.Legend(
fig[1, 3],
[l_max_dNdt, l_dNdt1, l_dNdt2, l_dNdt3, l_dNdt4],
[
"limit",
"ρₐ=1.2 kg/m3, Fᵣ=0.8, ρᵣ=800kg/m3",
"ρₐ=1.2 kg/m3, Fᵣ=0.2, ρᵣ=800kg/m3",
"ρₐ=1.2 kg/m3, Fᵣ=0.2, ρᵣ=200kg/m3",
"ρₐ=0.5 kg/m3, Fᵣ=0.2, ρᵣ=200kg/m3",
],
framevisible = false,
)
Makie.save("P3_ice_melt.svg", fig)
Makie.with_theme(Makie.theme_minimal(), fontsize = 22, linewidth = 3) do
fig = Makie.Figure(size = (800, 600))

ax1 = Makie.Axis(fig[1, 1];
xlabel = "T [°C]", ylabel = "ice mass melting rate [g/m³/s]",
title = "P3 ice melting",
)
ax2 = Makie.Axis(fig[1, 1];
xlabel = "T [°C]", ylabel = "ice number melting rate [1/cm³/s]",
yaxisposition = :right, rightspinevisible = true,
)

l_max_dLdt = Makie.hlines!(ax1, max_dLdt * 1e3; color = :gray, linewidth = 1)
l_max_dNdt = Makie.hlines!(ax2, max_dNdt * 1e-6; color = :gray, linewidth = 1, label = "limit")

l_dLdt1 = Makie.lines!(ax1, ΔT_range, dLdt1 * 1e3; color = :skyblue)
l_dNdt1 = Makie.lines!(ax2, ΔT_range, dNdt1 * 1e-6; color = :skyblue, label = label1)

l_dLdt2 = Makie.lines!(ax1, ΔT_range, dLdt2 * 1e3; color = :blue3)
l_dNdt2 = Makie.lines!(ax2, ΔT_range, dNdt2 * 1e-6; color = :blue3, label = label2)

l_dLdt3 = Makie.lines!(ax1, ΔT_range, dLdt3 * 1e3; color = :orchid)
l_dNdt3 = Makie.lines!(ax2, ΔT_range, dNdt3 * 1e-6; color = :orchid, label = label3)

l_dLdt4 = Makie.lines!(ax1, ΔT_range, dLdt4 * 1e3; color = :purple)
l_dNdt4 = Makie.lines!(ax2, ΔT_range, dNdt4 * 1e-6; color = :purple, label = label4)

Makie.axislegend(ax2; position = :rb, framevisible = false)
Makie.save("P3_ice_melt.svg", fig)
end
Loading
Loading