Skip to content

Commit

Permalink
refactor to prep for AIDA HET calibs
Browse files Browse the repository at this point in the history
  • Loading branch information
@amylu00
amylu00 committed Dec 5, 2024
1 parent ee7d3ac commit 7081d1a
Show file tree
Hide file tree
Showing 3 changed files with 102 additions and 120 deletions.
1 change: 1 addition & 0 deletions papers/Project.toml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@ ClimaUtilities = "b3f4f4ca-9299-4f7f-bd9b-81e1242a7513"
CloudMicrophysics = "6a9e3e04-43cd-43ba-94b9-e8782df3c71b"
Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
EnsembleKalmanProcesses = "aa8a2aa5-91d8-4396-bcef-d4f2ec43552d"
Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
Expand Down
219 changes: 99 additions & 120 deletions papers/ice_nucleation_2024/AIDA_calibrations.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,20 +14,34 @@ include(joinpath(pkgdir(CM), "papers", "ice_nucleation_2024", "calibration.jl"))

# Helper functions
function unpack_data(data_file_name)

file_path = CM.ArtifactCalling.AIDA_ice_nucleation(data_file_name)
return DelimitedFiles.readdlm(file_path, skipstart = 125)
end
function grab_data(unpacked_data)

AIDA_t_profile = unpacked_data[:, 1]
AIDA_T_profile = unpacked_data[:, 3]
AIDA_P_profile = unpacked_data[:, 2] * 1e2 # hPa to Pa
AIDA_ICNC = unpacked_data[:, 6] .* 1e6 # Nᵢ [m^-3]
AIDA_e = unpacked_data[:, 4]

data = (; AIDA_t_profile, AIDA_T_profile, AIDA_P_profile, AIDA_ICNC, AIDA_e)
return data
end
moving_average(data, n) =
[sum(@view data[i:(i + n - 1)]) / n for i in 1:(length(data) - (n - 1))]
[sum(@view data[i:(i + n)]) / n for i in 1:(length(data) - n)]

# Defining data names, start/end times, etc.
# data_file_names = ["in05_17_aida.edf", "in05_18_aida.edf", "in05_19_aida.edf", "in07_01_aida.edf", "in07_19_aida.edf"]
# plot_names = ["IN0517", "IN0518", "IN0519", "IN0701", "IN0719"]
data_file_names = ["in05_17_aida.edf", "in05_18_aida.edf", "in05_19_aida.edf"]
plot_names = ["IN0517", "IN0518", "IN0519"]
end_sim = 25 # Loss func looks at last end_sim timesteps only
start_time_list = [195, 180, 80] .+ 100 # AIDA time starts at -100 seconds. Add freezing onset as needed.
end_time_list = [290, 290, 170] .+ 100 # approximate time freezing stops
moving_average_n = 20 # average every n points
updrafts = [FT(1.5), FT(1.4), FT(5)] # updrafts matching AIDA cooling rate
end_sim = 25 # Loss func looks at last end_sim timesteps only
start_time_index_list = [195, 180, 80, 50, 35] .+ 100 # AIDA time starts at -100 seconds. Add freezing onset as needed.
end_time_index_list = [290, 290, 170, 800, 800] .+ 100 # approximate time freezing stops
moving_average_n = 20 # average every n points
updrafts = [FT(1.5), FT(1.4), FT(5), FT(1.5), FT(1.5)] # updrafts matching AIDA cooling rate

# Additional definitions
tps = TD.Parameters.ThermodynamicsParameters(FT)
Expand All @@ -41,42 +55,44 @@ for (exp_index, data_file_name) in enumerate(data_file_names)
### Unpacking experiment-specific variables.
plot_name = plot_names[exp_index]
w = updrafts[exp_index]
start_time = start_time_list[exp_index]
end_time = end_time_list[exp_index]
start_time_index = start_time_index_list[exp_index]
start_time = start_time_index - 100
end_time_index = end_time_index_list[exp_index]
end_time = end_time_index - 100

## Moving average to smooth data.
moving_average_start = Int32(start_time + (moving_average_n / 2 - 1))
moving_average_end = Int32(end_time - (moving_average_n / 2))
t_max =
(end_time - 100 - (moving_average_n / 2)) - # AIDA time starts at -100 seconds (t = 0 at index 101)
(start_time - 100 + (moving_average_n / 2 - 1)) # duration of simulation
moving_average_start_index = Int32(start_time_index + (moving_average_n / 2))
moving_average_end_index = Int32(end_time_index - (moving_average_n / 2))
t_max = moving_average_end_index - moving_average_start_index
# (end_time - (moving_average_n / 2)) - # AIDA time starts at -100 seconds (t = 0 at index 101)
# (start_time + (moving_average_n / 2 - 1)) # duration of simulation

### 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]
t_profile = chamber_data[moving_average_start+1:moving_average_end, 1] .- moving_average_start .+ 100
T_profile = chamber_data[moving_average_start+1:moving_average_end, 3]
P_profile = chamber_data[moving_average_start+1:moving_average_end, 2] * 1e2
AIDA_data = grab_data(chamber_data)
(; AIDA_t_profile, AIDA_T_profile, AIDA_P_profile, AIDA_ICNC, AIDA_e) = AIDA_data
t_profile = AIDA_t_profile[moving_average_start_index:moving_average_end_index] .- moving_average_start_index .+ 100
T_profile = AIDA_T_profile[moving_average_start_index:moving_average_end_index]
P_profile = AIDA_P_profile[moving_average_start_index:moving_average_end_index]
ICNC_profile = AIDA_ICNC[moving_average_start_index:moving_average_end_index]
e_profile = AIDA_e[moving_average_start_index:moving_average_end_index]

params = AIDA_IN05_params(FT, w, t_max, t_profile, T_profile, P_profile)
IC = AIDA_IN05_IC(FT, data_file_name)

frozen_frac = ICNC ./ (IC[7] + IC[8] + IC[9])
Nₜ = IC[7] + IC[8] + IC[9]
frozen_frac = AIDA_ICNC[start_time_index:end_time_index] ./ Nₜ
frozen_frac_moving_mean = moving_average(frozen_frac, moving_average_n)
ICNC_moving_avg = moving_average(ICNC, moving_average_n)
ICNC_moving_avg = moving_average(AIDA_ICNC[start_time_index:end_time_index], moving_average_n)
Γ = 0.1 * LinearAlgebra.I * (maximum(frozen_frac_moving_mean) - minimum(frozen_frac_moving_mean)) # coeff is an estimated of the noise

## Calculating some variables from AIDA data.
chamber_S_l = zeros(FT, size(chamber_data[start_time:end_time, 4], 1), size(chamber_data[start_time:end_time, 4], 2))
for (i, e) in enumerate(chamber_data[start_time:end_time, 4])
temp = chamber_data[start_time:end_time, 3][i]
## Calculating saturation wrt liquid from AIDA data.
chamber_S_l = zeros(FT, size(e_profile, 1), size(e_profile, 2))
for (i, e) in enumerate(e_profile)
temp = T_profile[i]
eₛ = TD.saturation_vapor_pressure(tps, temp, TD.Liquid())
chamber_S_l[i] = e / eₛ
end
chamber_r_l = zeros(FT, size(chamber_data[start_time:end_time, 8], 1), size(chamber_data[start_time:end_time, 8], 2))
for (i, r_l) in enumerate(chamber_data[start_time:end_time, 8] .* 1e-6)
chamber_r_l[i] = r_l < 0 ? FT(0) : r_l
end

### Calibration.
EKI_output = calibrate_J_parameters_EKI(FT, "ABHOM", params, IC, frozen_frac_moving_mean, end_sim, Γ)
Expand All @@ -98,23 +114,8 @@ 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 = "ICNC [m^-3]", xlabel = "time [s]", title = "AIDA data $plot_name")
MK.lines!(
ax1,
chamber_data[start_time:end_time, 1],
ICNC,
label = "Raw AIDA",
color =:blue,
linestyle =:dash,
linewidth = 2,
)
MK.lines!(
ax1,
chamber_data[moving_average_start:moving_average_end, 1],
ICNC_moving_avg,
label = "AIDA moving mean",
linewidth = 2.5,
color =:blue,
)
MK.lines!(ax1, AIDA_t_profile, AIDA_ICNC, label = "Raw AIDA", color =:blue, linestyle =:dash, linewidth = 2)
MK.lines!(ax1, t_profile, ICNC_moving_avg, label = "AIDA moving mean", linewidth = 2.5, color =:blue)
MK.axislegend(ax1, framevisible = true, labelsize = 12, position = :rc)
MK.save("$plot_name"*"_ICNC.svg", AIDA_data_fig)

Expand All @@ -139,65 +140,43 @@ for (exp_index, data_file_name) in enumerate(data_file_names)
ax_parcel_7 = MK.Axis(calibrated_parcel_fig[1, 3], ylabel = "pressure [Pa]", xlabel = "time [s]")
ax_parcel_8 = MK.Axis(calibrated_parcel_fig[2, 3], ylabel = "Nₐ [m^-3]", xlabel = "time [s]")

MK.lines!(ax_parcel_1, EKI_parcel.t .+ (moving_average_n / 2), EKI_parcel[1, :], label = "EKI Calib Liq", color = :orange) # label = "liquid"
MK.lines!(ax_parcel_1, UKI_parcel.t .+ (moving_average_n / 2), UKI_parcel[1, :], label = "UKI Calib Liq", color = :fuchsia) # label = "liquid"
MK.lines!(ax_parcel_1, parcel_default.t .+ (moving_average_n / 2), parcel_default[1, :], label = "default", color = :darkorange2)
MK.lines!(ax_parcel_1, EKI_parcel.t .+ (moving_average_n / 2), S_i.(tps, EKI_parcel[3, :], EKI_parcel[1, :]), label = "EKI Calib Ice", color = :green)
# MK.lines!(
# ax_parcel_1,
# chamber_data[start_time:end_time, 1] .- (start_time - 100),
# vec(chamber_S_l),
# label = "chamber",
# color = :blue
# )

MK.lines!(ax_parcel_2, EKI_parcel.t .+ (moving_average_n / 2), EKI_parcel[5, :], color = :orange)
MK.lines!(ax_parcel_2, UKI_parcel.t .+ (moving_average_n / 2), UKI_parcel[5, :], color = :fuchsia)
MK.lines!(ax_parcel_2, parcel_default.t .+ (moving_average_n / 2), parcel_default[5,:], color = :darkorange2)
MK.lines!(ax_parcel_1, EKI_parcel.t, EKI_parcel[1, :], label = "EKI Calib Liq", color = :orange) # label = "liquid"
MK.lines!(ax_parcel_1, UKI_parcel.t, UKI_parcel[1, :], label = "UKI Calib Liq", color = :fuchsia) # label = "liquid"
MK.lines!(ax_parcel_1, parcel_default.t, parcel_default[1, :], label = "default", color = :darkorange2)
MK.lines!(ax_parcel_1, EKI_parcel.t, S_i.(tps, EKI_parcel[3, :], EKI_parcel[1, :]), label = "EKI Calib Ice", color = :green)
# MK.lines!(ax_parcel_1, t_profile, vec(chamber_S_l), label = "chamber", color = :blue)

MK.lines!(ax_parcel_3, EKI_parcel.t .+ (moving_average_n / 2), EKI_parcel[3, :], color = :orange)
MK.lines!(ax_parcel_3, UKI_parcel.t .+ (moving_average_n / 2), UKI_parcel[3, :], color = :fuchsia)
MK.lines!(ax_parcel_3, parcel_default.t .+ (moving_average_n / 2), parcel_default[3, :], color = :darkorange2)
MK.lines!(
ax_parcel_3,
chamber_data[start_time:end_time, 1] .- (start_time - 100),
chamber_data[start_time:end_time, 3],
color = :blue,
linestyle =:dash,
)
MK.lines!(ax_parcel_2, EKI_parcel.t, EKI_parcel[5, :], color = :orange)
MK.lines!(ax_parcel_2, UKI_parcel.t, UKI_parcel[5, :], color = :fuchsia)
MK.lines!(ax_parcel_2, parcel_default.t, parcel_default[5,:], color = :darkorange2)

MK.lines!(ax_parcel_4, EKI_parcel.t .+ (moving_average_n / 2), EKI_parcel[6, :], color = :orange)
MK.lines!(ax_parcel_4, UKI_parcel.t .+ (moving_average_n / 2), UKI_parcel[6, :], color = :fuchsia)
MK.lines!(ax_parcel_4, parcel_default.t .+ (moving_average_n / 2), parcel_default[6, :], color = :darkorange2)
MK.lines!(ax_parcel_3, EKI_parcel.t, EKI_parcel[3, :], color = :orange)
MK.lines!(ax_parcel_3, UKI_parcel.t, UKI_parcel[3, :], color = :fuchsia)
MK.lines!(ax_parcel_3, parcel_default.t, parcel_default[3, :], color = :darkorange2)
MK.lines!(ax_parcel_3, t_profile, T_profile, color = :blue, linestyle =:dash)

MK.lines!(ax_parcel_5, EKI_parcel.t .+ (moving_average_n / 2), EKI_parcel[8, :], color = :orange)
MK.lines!(ax_parcel_5, UKI_parcel.t .+ (moving_average_n / 2), UKI_parcel[8, :], color = :fuchsia)
MK.lines!(ax_parcel_5, parcel_default.t .+ (moving_average_n / 2), parcel_default[8, :], color = :darkorange2)
MK.lines!(ax_parcel_4, EKI_parcel.t, EKI_parcel[6, :], color = :orange)
MK.lines!(ax_parcel_4, UKI_parcel.t, UKI_parcel[6, :], color = :fuchsia)
MK.lines!(ax_parcel_4, parcel_default.t, parcel_default[6, :], color = :darkorange2)

MK.lines!(ax_parcel_6, EKI_parcel.t .+ (moving_average_n / 2), EKI_parcel[9, :], color = :orange, label = "EKI")
MK.lines!(ax_parcel_6, UKI_parcel.t .+ (moving_average_n / 2), UKI_parcel[9, :], color = :fuchsia, label = "UKI")
MK.lines!(ax_parcel_6, parcel_default.t .+ (moving_average_n / 2), parcel_default[9, :], color = :darkorange2, label = "Pre-Calib")
MK.lines!(ax_parcel_6,
chamber_data[start_time:end_time, 1] .- (start_time - 100),
ICNC,
color = :blue,
label = "AIDA",
)
error = fill(sqrt(Γ[1,1]) * 2, length(chamber_data[start_time:end_time, 1] .- (start_time - 100)))
MK.errorbars!(ax_parcel_6, chamber_data[start_time:end_time, 1] .- (start_time - 100), ICNC ./ (IC[7] + IC[8] + IC[9]), error)
MK.lines!(ax_parcel_5, EKI_parcel.t, EKI_parcel[8, :], color = :orange)
MK.lines!(ax_parcel_5, UKI_parcel.t, UKI_parcel[8, :], color = :fuchsia)
MK.lines!(ax_parcel_5, parcel_default.t, parcel_default[8, :], color = :darkorange2)

MK.lines!(ax_parcel_7, EKI_parcel.t .+ (moving_average_n / 2), EKI_parcel[2, :], color = :orange)
MK.lines!(ax_parcel_7, UKI_parcel.t .+ (moving_average_n / 2), UKI_parcel[2, :], color = :fuchsia)
MK.lines!(
ax_parcel_7,
chamber_data[start_time:end_time, 1] .- (start_time - 100),
chamber_data[start_time:end_time, 2] .* 1e2,
color = :blue,
linestyle =:dash,
)
MK.lines!(ax_parcel_6, EKI_parcel.t, EKI_parcel[9, :], color = :orange, label = "EKI")
MK.lines!(ax_parcel_6, UKI_parcel.t, UKI_parcel[9, :], color = :fuchsia, label = "UKI")
MK.lines!(ax_parcel_6, parcel_default.t, parcel_default[9, :], color = :darkorange2, label = "Pre-Calib")
MK.lines!(ax_parcel_6, t_profile, ICNC_profile, color = :blue, label = "AIDA",)

error = fill(sqrt(Γ[1,1]) * 2, length(AIDA_t_profile[start_time_index:end_time_index] .- start_time))
MK.errorbars!(ax_parcel_6, AIDA_t_profile[start_time_index:end_time_index] .- start_time, AIDA_ICNC[start_time_index:end_time_index] ./ Nₜ, error)

MK.lines!(ax_parcel_7, EKI_parcel.t, EKI_parcel[2, :], color = :orange)
MK.lines!(ax_parcel_7, UKI_parcel.t, UKI_parcel[2, :], color = :fuchsia)
MK.lines!(ax_parcel_7, t_profile, P_profile, color = :blue) #, linestyle =:dash

MK.lines!(ax_parcel_8, EKI_parcel.t .+ (moving_average_n / 2), EKI_parcel[7, :], color = :orange)
MK.lines!(ax_parcel_8, UKI_parcel.t .+ (moving_average_n / 2), UKI_parcel[7, :], color = :fuchsia)
MK.lines!(ax_parcel_8, EKI_parcel.t, EKI_parcel[7, :], color = :orange)
MK.lines!(ax_parcel_8, UKI_parcel.t, UKI_parcel[7, :], color = :fuchsia)

MK.axislegend(ax_parcel_6, framevisible = false, labelsize = 16, position = :rb)
MK.save("$plot_name"*"_calibrated_parcel_fig.svg", calibrated_parcel_fig)
Expand All @@ -208,31 +187,31 @@ for (exp_index, data_file_name) in enumerate(data_file_names)
ax_compare = MK.Axis(ICNC_comparison_fig[1, 1], ylabel = "Frozen Fraction [-]", xlabel = "time [s]", title = "$plot_name")
# MK.lines!(
# ax_compare,
# parcel_default.t .+ (moving_average_n / 2),
# parcel_default[9, :]./ (IC[7] + IC[8] + IC[9]),
# parcel_default.t,
# parcel_default[9, :]./ Nₜ,
# label = "CM.jl Parcel (Pre-Calib)",
# linewidth = 2.5,
# color =:orangered,
# )
MK.lines!(
ax_compare,
EKI_parcel.t .+ (moving_average_n / 2),
EKI_parcel[9, :]./ (IC[7] + IC[8] + IC[9]),
EKI_parcel.t,
EKI_parcel[9, :]./ Nₜ,
label = "CM.jl Parcel (EKI Calibrated)",
linewidth = 2.5,
color =:orange,
)
MK.lines!(
ax_compare,
UKI_parcel.t .+ (moving_average_n / 2),
UKI_parcel[9, :]./ (IC[7] + IC[8] + IC[9]),
UKI_parcel.t,
UKI_parcel[9, :]./ Nₜ,
label = "CM.jl Parcel (UKI Calibrated)",
linewidth = 2.5,
color =:fuchsia,
)
MK.lines!(
ax_compare,
chamber_data[start_time:end_time, 1] .- (start_time - 100),
AIDA_t_profile[start_time_index:end_time_index] .- start_time,
frozen_frac,
label = "Raw AIDA",
linewidth = 2,
Expand All @@ -241,14 +220,14 @@ for (exp_index, data_file_name) in enumerate(data_file_names)
)
MK.lines!(
ax_compare,
chamber_data[moving_average_start:moving_average_end, 1] .- (start_time - 100),
t_profile,
frozen_frac_moving_mean,
label = "AIDA Moving Avg",
linewidth = 2.5,
color =:blue
)
error = fill(sqrt(Γ[1,1]) * 2, length(frozen_frac_moving_mean))
MK.errorbars!(ax_compare, chamber_data[moving_average_start:moving_average_end, 1] .- (start_time - 100), frozen_frac_moving_mean, error)
MK.errorbars!(ax_compare, t_profile, frozen_frac_moving_mean, error)

MK.axislegend(ax_compare, framevisible = false, labelsize = 20, position = :rb)
MK.save("$plot_name"*"_ICNC_comparison_fig.svg", ICNC_comparison_fig)
Expand All @@ -265,58 +244,58 @@ for (exp_index, data_file_name) in enumerate(data_file_names)
ax_spread = MK.Axis(UKI_spread_fig[1, 1], ylabel = "Frozen Fraction [-]", xlabel = "time [s]", title = "$plot_name")
MK.lines!(
ax_spread,
chamber_data[moving_average_start:moving_average_end, 1] .- (start_time - 100),
t_profile,
frozen_frac_moving_mean,
label = "AIDA Moving Avg",
linewidth = 2.5,
color =:blue
)
MK.lines!(
ax_spread,
UKI_parcel_1.t .+ (moving_average_n / 2),
UKI_parcel_1[9, :]./ (IC[7] + IC[8] + IC[9]),
UKI_parcel_1.t,
UKI_parcel_1[9, :]./ Nₜ,
linewidth = 2.5,
color =:grey80,
)
MK.lines!(
ax_spread,
UKI_parcel_2.t .+ (moving_average_n / 2),
UKI_parcel_2[9, :]./ (IC[7] + IC[8] + IC[9]),
UKI_parcel_2.t,
UKI_parcel_2[9, :]./ Nₜ,
linewidth = 2.5,
color =:grey60,
)
MK.lines!(
ax_spread,
UKI_parcel_3.t .+ (moving_average_n / 2),
UKI_parcel_3[9, :]./ (IC[7] + IC[8] + IC[9]),
UKI_parcel_3.t,
UKI_parcel_3[9, :]./ Nₜ,
linewidth = 2.5,
color =:grey40,
)
MK.lines!(
ax_spread,
UKI_parcel_4.t .+ (moving_average_n / 2),
UKI_parcel_4[9, :]./ (IC[7] + IC[8] + IC[9]),
UKI_parcel_4.t,
UKI_parcel_4[9, :]./ Nₜ,
linewidth = 2.5,
color =:grey25,
)
MK.lines!(
ax_spread,
UKI_parcel_5.t .+ (moving_average_n / 2),
UKI_parcel_5[9, :]./ (IC[7] + IC[8] + IC[9]),
UKI_parcel_5.t,
UKI_parcel_5[9, :]./ Nₜ,
linewidth = 2.5,
color =:grey12,
)
MK.lines!(
ax_spread,
UKI_parcel.t .+ (moving_average_n / 2),
UKI_parcel[9, :]./ (IC[7] + IC[8] + IC[9]),
UKI_parcel.t,
UKI_parcel[9, :]./ Nₜ,
label = "CM.jl Parcel (UKI Calibrated)",
linewidth = 2.5,
color =:fuchsia,
linestyle = :dash,
)
error = fill(sqrt(Γ[1,1]) * 2, length(frozen_frac_moving_mean))
MK.errorbars!(ax_spread, chamber_data[moving_average_start:moving_average_end, 1] .- (start_time - 100), frozen_frac_moving_mean, error)
MK.errorbars!(ax_spread, t_profile, frozen_frac_moving_mean, error)
MK.save("$plot_name"*"_UKI_spread_fig.svg", UKI_spread_fig)

#! format: on
Expand Down
2 changes: 2 additions & 0 deletions parcel/ParcelCommon.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,8 @@ S_i(tps, T, S_liq) = ξ(tps, T) * S_liq

# Interpolating for cooling/expansion rate
function AIDA_rate(model_t, data_t, data)
data_t = data_t[2:end]
data = data[2:end]
index = Int32(model_t) + 1
data_at_t = Intp.linear_interpolation(data_t, data)

Expand Down

0 comments on commit 7081d1a

Please sign in to comment.