From 7081d1a57c9bca42daf083eaec9fa97770cb1786 Mon Sep 17 00:00:00 2001 From: amylu00 Date: Wed, 4 Dec 2024 16:38:04 -0800 Subject: [PATCH] refactor to prep for AIDA HET calibs --- papers/Project.toml | 1 + .../ice_nucleation_2024/AIDA_calibrations.jl | 219 ++++++++---------- parcel/ParcelCommon.jl | 2 + 3 files changed, 102 insertions(+), 120 deletions(-) diff --git a/papers/Project.toml b/papers/Project.toml index 3438a730f2..a603753994 100644 --- a/papers/Project.toml +++ b/papers/Project.toml @@ -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" diff --git a/papers/ice_nucleation_2024/AIDA_calibrations.jl b/papers/ice_nucleation_2024/AIDA_calibrations.jl index 4adbceb2a0..88bd1babd5 100644 --- a/papers/ice_nucleation_2024/AIDA_calibrations.jl +++ b/papers/ice_nucleation_2024/AIDA_calibrations.jl @@ -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) @@ -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, Γ) @@ -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) @@ -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) @@ -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, @@ -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) @@ -265,7 +244,7 @@ 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, @@ -273,50 +252,50 @@ for (exp_index, data_file_name) in enumerate(data_file_names) ) 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 diff --git a/parcel/ParcelCommon.jl b/parcel/ParcelCommon.jl index 3c70661f74..756575f01e 100644 --- a/parcel/ParcelCommon.jl +++ b/parcel/ParcelCommon.jl @@ -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)