Skip to content

Commit

Permalink
cleaning before the great squash
Browse files Browse the repository at this point in the history
  • Loading branch information
amylu00 committed Sep 17, 2024
1 parent 857326f commit 20f85ea
Show file tree
Hide file tree
Showing 3 changed files with 8 additions and 25 deletions.
12 changes: 5 additions & 7 deletions papers/ice_nucleation_2024/AIDA_calibrations.jl
Original file line number Diff line number Diff line change
Expand Up @@ -54,7 +54,7 @@ for (exp_index, data_file_name) in enumerate(data_file_names)

### Check for data file in AIDA_data folder.
chamber_data = unpack_data(data_file_name)
ICNC = chamber_data[start_time:end_time, 6] .* 1e6 # Nᵢ [m^-3] ; freezing onset to when it starts decreasing [195:390]
ICNC = chamber_data[start_time:end_time, 6] .* 1e6 # Nᵢ [m^-3]

frozen_frac = ICNC ./ (IC[7] + IC[8] + IC[9])
frozen_frac_moving_mean = moving_average(frozen_frac, moving_average_n)
Expand All @@ -63,7 +63,7 @@ for (exp_index, data_file_name) in enumerate(data_file_names)

### Plotting AIDA data.
AIDA_data_fig = MK.Figure(size = (800, 600), fontsize = 24)
ax1 = MK.Axis(AIDA_data_fig[1, 1], ylabel = "Frozen Fraction [-]", xlabel = "time [s]", title = "AIDA data $plot_name")
ax1 = MK.Axis(AIDA_data_fig[1, 1], ylabel = "ICNC [m^-3]", xlabel = "time [s]", title = "AIDA data $plot_name")
MK.lines!(
ax1,
chamber_data[start_time:end_time, 1],
Expand Down Expand Up @@ -98,9 +98,8 @@ for (exp_index, data_file_name) in enumerate(data_file_names)
## Calibrated coefficients.
# Did they converge?
calibrated_coeffs_fig = MK.Figure(size = (600, 600), fontsize = 24)
ax3 = MK.Axis(calibrated_coeffs_fig[1, 1], ylabel = "m coefficient [-]")#, title = "$plot_name")
ax3 = MK.Axis(calibrated_coeffs_fig[1, 1], ylabel = "m coefficient [-]", title = "$plot_name")
ax4 = MK.Axis(calibrated_coeffs_fig[2, 1], ylabel = "c coefficient [-]", xlabel = "iteration #")
#MK.ylims!(ax3, 250, 260)
MK.lines!(ax3, collect(1:n_iterations), calibrated_ensemble_means[1], label = "ensemble mean", color = :orange, linewidth = 2.5)
MK.lines!(ax4, collect(1:n_iterations), calibrated_ensemble_means[2], label = "ensemble mean", color = :orange, linewidth = 2.5)
MK.save("$plot_name"*"_calibrated_coeffs_fig.svg", calibrated_coeffs_fig)
Expand Down Expand Up @@ -173,19 +172,18 @@ for (exp_index, data_file_name) in enumerate(data_file_names)
# Does calibrated parcel look like observations?
ICNC_comparison_fig = MK.Figure(size = (700, 600), fontsize = 24)
ax_compare = MK.Axis(ICNC_comparison_fig[1, 1], ylabel = "Frozen Fraction [-]", xlabel = "time [s]", title = "$plot_name")
#MK.ylims!(ax_compare, 0, 0.14)
# MK.lines!(
# ax_compare,
# parcel_default.t .+ (moving_average_n / 2),
# parcel_default[9, :]./ (IC[7] + IC[8] + IC[9]), # IC[8] if starting with droplets (no aerosol activation)
# parcel_default[9, :]./ (IC[7] + IC[8] + IC[9]),
# label = "CM.jl Parcel (Pre-Calib)",
# linewidth = 2.5,
# color =:orangered,
# )
MK.lines!(
ax_compare,
parcel.t .+ (moving_average_n / 2),
parcel[9, :]./ (IC[7] + IC[8] + IC[9]), # IC[8] if starting with droplets (no aerosol activation)
parcel[9, :]./ (IC[7] + IC[8] + IC[9]),
label = "CM.jl Parcel (Calibrated)",
linewidth = 2.5,
color =:orange,
Expand Down
13 changes: 1 addition & 12 deletions papers/ice_nucleation_2024/calibration.jl
Original file line number Diff line number Diff line change
Expand Up @@ -47,10 +47,6 @@ function run_model(p, coefficients, IN_mode, FT, IC)
ice_size_distribution = ice_size_distribution,
)

# solve ODE
local sol = run_parcel(IC, FT(0), t_max, params)
return sol[9, :] ./ (IC[7] + IC[8] + IC[9]) # frozen fraction

elseif IN_mode == "ABIFM"
# overwriting
override_file = Dict(
Expand All @@ -74,10 +70,6 @@ function run_model(p, coefficients, IN_mode, FT, IC)
ice_size_distribution = ice_size_distribution,
)

# solve ODE
local sol = run_parcel(IC, FT(0), t_max, params)
return sol[9, :] ./ (IC[7] + IC[8] + IC[9]) # frozen fraction

elseif IN_mode == "ABHOM"
# overwriting
override_file = Dict(
Expand All @@ -90,7 +82,6 @@ function run_model(p, coefficients, IN_mode, FT, IC)
overwrite = CMP.IceNucleationParameters(ip_calibrated)

# run parcel with new coefficients
# TODO - need better way to call σ_g. Rename σ_g to σ_g_liq or something.
local params = parcel_params{FT}(
const_dt = const_dt,
w = w,
Expand All @@ -105,12 +96,10 @@ function run_model(p, coefficients, IN_mode, FT, IC)
ips = overwrite,
r_nuc = r_nuc,
)


end
# solve ODE
local sol = run_parcel(IC, FT(0), t_max, params)
return sol[9, :] ./ (IC[7] + IC[8] + IC[9]) # frozen fraction
end
end

function run_calibrated_model(FT, IN_mode, coefficients, p, IC)
Expand Down
8 changes: 2 additions & 6 deletions papers/ice_nucleation_2024/calibration_setup.jl
Original file line number Diff line number Diff line change
Expand Up @@ -75,7 +75,7 @@ function perf_model_IC(FT, IN_mode)
ϵₘ = R_d / R_v

if IN_mode == "ABDINM"
Nₐ = FT(2e8) # according to fig 2 panel b deposition nucleation experiment
Nₐ = FT(2e8)
Nₗ = FT(0)
Nᵢ = FT(0)
r₀ = FT(1e-7)
Expand All @@ -89,9 +89,6 @@ function perf_model_IC(FT, IN_mode)
eₛ = TD.saturation_vapor_pressure(tps, T₀, TD.Liquid())
e = eᵥ(qᵥ, p₀, Rₐ, R_v)
Sₗ = FT(e / eₛ)
Sᵢ = S_i(tps, T₀, Sₗ)
# @info("", Sᵢ)
# error("breaking")
elseif IN_mode == "ABIFM"
Nₐ = FT(0)
Nₗ = FT(2000)
Expand All @@ -115,7 +112,6 @@ function perf_model_IC(FT, IN_mode)
p₀ = FT(9712.183)
T₀ = FT(190)
qₗ = FT(Nₗ * 4 / 3 * FT(π) * r₀^3 * ρₗ / FT(1.2)) # 1.2 should be ρₐ
C_l = FT(qₗ / ((1 - qₗ) * ϵₘ + qₗ)) # concentration/mol fraction of liquid
C_v = FT(5 * 1e-6)
qᵥ = ϵₘ / (ϵₘ - 1 + 1 / C_v)
qᵢ = FT(0)
Expand Down Expand Up @@ -146,7 +142,7 @@ function perf_model_pseudo_data(FT, IN_mode, params, IC)
noise_dist = Distributions.MvNormal(zeros(dim_output), Γ)

y_truth = zeros(length(G_truth), n_samples) # where noisy ICNC will be stored
y_truth = (G_truth .+ rand(noise_dist))
y_truth = G_truth .+ rand(noise_dist)
return [y_truth, Γ, coeff_true]
end

Expand Down

0 comments on commit 20f85ea

Please sign in to comment.