diff --git a/papers/ice_nucleation_2024/AIDA_calibrations.jl b/papers/ice_nucleation_2024/AIDA_calibrations.jl index cdc3d7ea6e..ebeea265c7 100644 --- a/papers/ice_nucleation_2024/AIDA_calibrations.jl +++ b/papers/ice_nucleation_2024/AIDA_calibrations.jl @@ -33,13 +33,13 @@ moving_average(data, n) = [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 = ["in07_01_aida.edf", "in07_19_aida.edf"] +plot_names = ["IN0701", "IN0719"] end_sim = 25 # Loss func looks at last end_sim timesteps only -start_time_list = [195, 180, 80, 50, 35] # freezing onset -end_time_list = [290, 290, 170, 800, 800] # approximate time freezing stops +start_time_list = [50, 35] # freezing onset +end_time_list = [400, 400] # 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 +updrafts = [FT(1.5), FT(1.5)] # updrafts matching AIDA cooling rate # Additional definitions tps = TD.Parameters.ThermodynamicsParameters(FT) @@ -58,8 +58,10 @@ for (exp_index, data_file_name) in enumerate(data_file_names) end_time = end_time_list[exp_index] end_time_index = end_time .+ 100 - nuc_mode = plot_name == "IN0701" || plot_name == "IN0719" ? "HET" : "HOM" - + nuc_mode = plot_name == "IN0701" || plot_name == "IN0719" ? "ABDINM" : "ABHOM" + println("\n \n \n") + @info(exp_index, plot_name, nuc_mode) + sleep(5) ## Moving average to smooth data. moving_average_start_index = Int32(start_time_index + (moving_average_n / 2)) moving_average_end_index = Int32(end_time_index - (moving_average_n / 2)) @@ -77,11 +79,11 @@ for (exp_index, data_file_name) in enumerate(data_file_names) S_l_profile = e_profile ./ (TD.saturation_vapor_pressure.(tps, T_profile, TD.Liquid())) params = - nuc_mode == "HOM" ? + nuc_mode == "ABHOM" ? AIDA_IN05_params(FT, w, t_max, t_profile, T_profile, P_profile) : AIDA_IN07_params(FT, w, t_max, t_profile, T_profile, P_profile, plot_name) IC = - nuc_mode == "HOM" ? + nuc_mode == "ABHOM" ? AIDA_IN05_IC(FT, data_file_name) : AIDA_IN07_IC(FT, data_file_name) @@ -92,8 +94,8 @@ for (exp_index, data_file_name) in enumerate(data_file_names) Γ = 0.1 * LinearAlgebra.I * (maximum(frozen_frac_moving_mean) - minimum(frozen_frac_moving_mean)) # coeff is an estimated of the noise ### Calibration. - EKI_output = calibrate_J_parameters_EKI(FT, "ABHOM", params, IC, frozen_frac_moving_mean, end_sim, Γ) - UKI_output = calibrate_J_parameters_UKI(FT, "ABHOM", params, IC, frozen_frac_moving_mean, end_sim, Γ) + EKI_output = calibrate_J_parameters_EKI(FT, nuc_mode, params, IC, frozen_frac_moving_mean, end_sim, Γ) + UKI_output = calibrate_J_parameters_UKI(FT, nuc_mode, params, IC, frozen_frac_moving_mean, end_sim, Γ) EKI_n_iterations = size(EKI_output[2])[1] EKI_n_ensembles = size(EKI_output[2][1])[2] @@ -103,9 +105,9 @@ for (exp_index, data_file_name) in enumerate(data_file_names) calibrated_ensemble_means = ensemble_means(EKI_output[2], EKI_n_iterations, EKI_n_ensembles) ## Calibrated parcel. - EKI_parcel = run_calibrated_model(FT, "ABHOM", EKI_calibrated_parameters, params, IC) - UKI_parcel = run_calibrated_model(FT, "ABHOM", UKI_calibrated_parameters, params, IC) - parcel_default = run_calibrated_model(FT, "ABHOM", [FT(255.927125), FT(-68.553283)], params, IC) + EKI_parcel = run_calibrated_model(FT, nuc_mode, EKI_calibrated_parameters, params, IC) + UKI_parcel = run_calibrated_model(FT, nuc_mode, UKI_calibrated_parameters, params, IC) + parcel_default = run_calibrated_model(FT, nuc_mode, [FT(255.927125), FT(-68.553283)], params, IC) # TODO - these are HOM params ### Plots. ## Plotting AIDA data. @@ -231,11 +233,11 @@ for (exp_index, data_file_name) in enumerate(data_file_names) ## Looking at spread in UKI calibrated parameters ϕ_UKI = UKI_output[2] - UKI_parcel_1 = run_calibrated_model(FT, "ABHOM", [ϕ_UKI[1,1], ϕ_UKI[2,1]], params, IC) - UKI_parcel_2 = run_calibrated_model(FT, "ABHOM", [ϕ_UKI[1,2], ϕ_UKI[2,2]], params, IC) - UKI_parcel_3 = run_calibrated_model(FT, "ABHOM", [ϕ_UKI[1,3], ϕ_UKI[2,3]], params, IC) - UKI_parcel_4 = run_calibrated_model(FT, "ABHOM", [ϕ_UKI[1,4], ϕ_UKI[2,4]], params, IC) - UKI_parcel_5 = run_calibrated_model(FT, "ABHOM", [ϕ_UKI[1,5], ϕ_UKI[2,5]], params, IC) + UKI_parcel_1 = run_calibrated_model(FT, nuc_mode, [ϕ_UKI[1,1], ϕ_UKI[2,1]], params, IC) + UKI_parcel_2 = run_calibrated_model(FT, nuc_mode, [ϕ_UKI[1,2], ϕ_UKI[2,2]], params, IC) + UKI_parcel_3 = run_calibrated_model(FT, nuc_mode, [ϕ_UKI[1,3], ϕ_UKI[2,3]], params, IC) + UKI_parcel_4 = run_calibrated_model(FT, nuc_mode, [ϕ_UKI[1,4], ϕ_UKI[2,4]], params, IC) + UKI_parcel_5 = run_calibrated_model(FT, nuc_mode, [ϕ_UKI[1,5], ϕ_UKI[2,5]], params, IC) UKI_spread_fig = MK.Figure(size = (700, 600), fontsize = 24) ax_spread = MK.Axis(UKI_spread_fig[1, 1], ylabel = "Frozen Fraction [-]", xlabel = "time [s]", title = "$plot_name") diff --git a/papers/ice_nucleation_2024/calibration.jl b/papers/ice_nucleation_2024/calibration.jl index e7af60f1fd..e3b59ee391 100644 --- a/papers/ice_nucleation_2024/calibration.jl +++ b/papers/ice_nucleation_2024/calibration.jl @@ -29,14 +29,34 @@ function run_model(p, coefficients, IN_mode, FT, IC, end_sim) if IN_mode == "ABDINM" # overwriting - override_file = Dict( - "China2017_J_deposition_m_Kaolinite" => - Dict("value" => m_calibrated, "type" => "float"), - "China2017_J_deposition_c_Kaolinite" => - Dict("value" => c_calibrated, "type" => "float"), - ) - kaolinite_calibrated = CP.create_toml_dict(FT; override_file) - aerosol = CMP.Kaolinite(kaolinite_calibrated) + if aerosol == CMP.Kaolinite(FT) + override_file = Dict( + "China2017_J_deposition_m_Kaolinite" => + Dict("value" => m_calibrated, "type" => "float"), + "China2017_J_deposition_c_Kaolinite" => + Dict("value" => c_calibrated, "type" => "float"), + ) + kaolinite_calibrated = CP.create_toml_dict(FT; override_file) + aerosol = CMP.Kaolinite(kaolinite_calibrated) + elseif aerosol == CMP.ArizonaTestDust(FT) + override_file = Dict( + "J_ABDINM_m_ArizonaTestDust" => + Dict("value" => m_calibrated, "type" => "float"), + "J_ABDINM_c_ArizonaTestDust" => + Dict("value" => c_calibrated, "type" => "float"), + ) + ATD_calibrated = CP.create_toml_dict(FT; override_file) + aerosol = CMP.ArizonaTestDust(ATD_calibrated) + elseif aerosol == CMP.Illite(FT) + override_file = Dict( + "J_ABDINM_m_Illite" => + Dict("value" => m_calibrated, "type" => "float"), + "J_ABDINM_c_Illite" => + Dict("value" => c_calibrated, "type" => "float"), + ) + illite_calibrated = CP.create_toml_dict(FT; override_file) + aerosol = CMP.Illite(illite_calibrated) + end elseif IN_mode == "ABIFM" # overwriting @@ -100,14 +120,34 @@ function run_calibrated_model(FT, IN_mode, coefficients, p, IC) if IN_mode == "ABDINM" # overwriting - override_file = Dict( - "China2017_J_deposition_m_Kaolinite" => - Dict("value" => m_calibrated, "type" => "float"), - "China2017_J_deposition_c_Kaolinite" => - Dict("value" => c_calibrated, "type" => "float"), - ) - kaolinite_calibrated = CP.create_toml_dict(FT; override_file) - aerosol = CMP.Kaolinite(kaolinite_calibrated) + if aerosol == CMP.Kaolinite(FT) + override_file = Dict( + "China2017_J_deposition_m_Kaolinite" => + Dict("value" => m_calibrated, "type" => "float"), + "China2017_J_deposition_c_Kaolinite" => + Dict("value" => c_calibrated, "type" => "float"), + ) + kaolinite_calibrated = CP.create_toml_dict(FT; override_file) + aerosol = CMP.Kaolinite(kaolinite_calibrated) + elseif aerosol == CMP.ArizonaTestDust(FT) + override_file = Dict( + "J_ABDINM_m_ArizonaTestDust" => + Dict("value" => m_calibrated, "type" => "float"), + "J_ABDINM_c_ArizonaTestDust" => + Dict("value" => c_calibrated, "type" => "float"), + ) + ATD_calibrated = CP.create_toml_dict(FT; override_file) + aerosol = CMP.ArizonaTestDust(ATD_calibrated) + elseif aerosol == CMP.Illite(FT) + override_file = Dict( + "J_ABDINM_m_Illite" => + Dict("value" => m_calibrated, "type" => "float"), + "J_ABDINM_c_Illite" => + Dict("value" => c_calibrated, "type" => "float"), + ) + illite_calibrated = CP.create_toml_dict(FT; override_file) + aerosol = CMP.Illite(illite_calibrated) + end elseif IN_mode == "ABIFM" # overwriting @@ -159,7 +199,7 @@ function run_calibrated_model(FT, IN_mode, coefficients, p, IC) return sol end -function create_prior(FT, IN_mode, ; perfect_model = false) +function create_prior(FT, IN_mode, ; perfect_model = false, aerosol_type = nothing) # TODO - add perfect_model flag to plot_ensemble_mean.jl observation_data_names = ["m_coeff", "c_coeff"] @@ -177,14 +217,20 @@ function create_prior(FT, IN_mode, ; perfect_model = false) c_stats = [FT(-66), FT(3), -Inf, Inf] end elseif perfect_model == false + @info("create_prior.jl", aerosol_type) + sleep(2) if IN_mode == "ABDINM" - # m_stats = [FT(20), FT(1), FT(0), Inf] - # c_stats = [FT(-1), FT(1), -Inf, Inf] - println("Calibration for ABDINM with AIDA not yet implemented.") + if aerosol_type == CMP.ArizonaTestDust(FT) + m_stats = [FT(13), FT(20), FT(0), Inf] + c_stats = [FT(0.7), FT(7), -Inf, Inf] + elseif aerosol_type == CMP.Illite(FT) + m_stats = [FT(13), FT(20), FT(0), Inf] + c_stats = [FT(0.7), FT(7), -Inf, Inf] + end elseif IN_mode == "ABIFM" # m_stats = [FT(50), FT(1), FT(0), Inf] # c_stats = [FT(-7), FT(1), -Inf, Inf] - println("Calibration for ABIFM with AIDA not yet implemented.") + error("ABIFM calibrations not yet supported.") elseif IN_mode == "ABHOM" m_stats = [FT(260.927125), FT(25), FT(0), Inf] c_stats = [FT(-68.553283), FT(10), -Inf, Inf] @@ -214,7 +260,10 @@ function calibrate_J_parameters_EKI(FT, IN_mode, params, IC, y_truth, end_sim, rng_seed = 24 rng = Random.seed!(Random.GLOBAL_RNG, rng_seed) - prior = create_prior(FT, IN_mode, perfect_model = perfect_model) + (; aerosol) = params + + + prior = create_prior(FT, IN_mode, perfect_model = perfect_model, aerosol_type = aerosol) N_ensemble = 25 # runs N_ensemble trials per iteration N_iterations = 50 # number of iterations the inverse problem goes through @@ -269,7 +318,9 @@ function calibrate_J_parameters_EKI(FT, IN_mode, params, IC, y_truth, end_sim, end function calibrate_J_parameters_UKI(FT, IN_mode, params, IC, y_truth, end_sim, Γ,; perfect_model = false) - prior = create_prior(FT, IN_mode, perfect_model = perfect_model) + (; aerosol) = params + + prior = create_prior(FT, IN_mode, perfect_model = perfect_model, aerosol_type = aerosol) N_iterations = 25 α_reg = 1.0 update_freq = 1 diff --git a/papers/ice_nucleation_2024/calibration_setup.jl b/papers/ice_nucleation_2024/calibration_setup.jl index 253b8ad7f3..ef63b94ec6 100644 --- a/papers/ice_nucleation_2024/calibration_setup.jl +++ b/papers/ice_nucleation_2024/calibration_setup.jl @@ -253,8 +253,10 @@ function AIDA_IN07_params(FT, w, t_max, t_profile, T_profile, P_profile, plot_na IN_mode = "ABDINM" const_dt = FT(1) prescribed_thermodynamics = true - aerosol_act = "AeroAct" + aerosol_act = "None" aerosol = plot_name == "IN0701" ? CMP.ArizonaTestDust(FT) : CMP.Illite(FT) + @info(plot_name, aerosol) + sleep(2) dep_nucleation = "ABDINM" heterogeneous = "None" homogeneous = "None" @@ -297,7 +299,7 @@ function AIDA_IN07_IC(FT, data_file) qᵢ = FT(Nᵢ * 4 / 3 * FT(π) * r₀^3 * ρₗ / FT(1.2)) m_l = Nₗ * ρₗ * 4 * π / 3 * r₀^3 m_i = Nᵢ * ρᵢ * 4 * π / 3 * r₀^3 - e = FT(0.067935) + e = FT(0.67935) qᵥ = (e / R_v / T₀) / ((p₀ - e) / (R_d * T₀) + e / R_v / T₀ + m_l + m_i) q = TD.PhasePartition.(qᵥ + qₗ + qᵢ, qₗ, qᵢ) Rₐ = TD.gas_constant_air(tps, q) diff --git a/parcel/ParcelModel.jl b/parcel/ParcelModel.jl index 070ca24660..f9983c0a01 100644 --- a/parcel/ParcelModel.jl +++ b/parcel/ParcelModel.jl @@ -141,6 +141,7 @@ function parcel_model(dY, Y, p, t) dqᵥ_dt = -dqₗ_dt_v2l - dqᵢ_dt_v2i dSₗ_dt = + #prescribed_thermodynamics ? - AIDA_rate(t, t_profile, T_profile) / 0.0065 : a1 * w * Sₗ - (a2 + a3) * Sₗ * dqₗ_dt_v2l - (a2 + a4) * Sₗ * dqᵢ_dt_v2i - a5 * Sₗ * dqᵢ_dt_l2i