Merge pull request #552 from CliMA/al/p3

Adding P3 to compare w AIDA calibs

Author: amylu00

papers/ice_nucleation_2024/AIDA_calibrations.jl

@@ -32,7 +32,6 @@ data_file_names = [
     ["ACI04_22", "EXP19"],
 ]
 batch_names = ["HOM", "IN0701", "IN0719", "ACI04_22", "EXP19", "EXP45", "DEP", "IMM"]
-end_sim = 25            # Loss func looks at last end_sim timesteps only
 start_time_list = [     # freezing onset
     [Int32(150), Int32(180), Int32(0)],
     [Int32(50)],
@@ -75,12 +74,13 @@ global EKI_calibrated_coeff_dict = Dict()
 global UKI_calibrated_coeff_dict = Dict()
 mutable struct all_calibration_data
     UKI_calibrated_parcel::Vector{ODE.ODESolution}
+    P3_parcel::Vector{ODE.ODESolution}
     Nₜ_list::Vector{Any}
     t_profile_list::Vector{Any}
     frozen_frac_moving_mean_list::Vector{Any}
     frozen_frac_list::Vector{Any}
 end
-global overview_data = all_calibration_data(Vector{ODE.ODESolution}(), [], [], [], [])
+global overview_data = all_calibration_data(Vector{ODE.ODESolution}(), Vector{ODE.ODESolution}(), [], [], [], [])
 
 for (batch_index, batch_name) in enumerate(batch_names)
     @info(batch_name)
@@ -94,23 +94,28 @@ for (batch_index, batch_name) in enumerate(batch_names)
 
     if batch_name == "HOM"
         nuc_mode = "ABHOM"
+        P3_nuc_mode  = "P3_hom"
     elseif batch_name in ["IN0701", "IN0719", "EXP45", "DEP"]
         nuc_mode = "ABDINM"
+        P3_nuc_mode = "P3_dep"
     elseif batch_name in ["ACI04_22", "EXP19", "IMM"]
         nuc_mode = "ABIFM"
+        P3_nuc_mode = "P3_het"
     end
 
     ### Check for and grab data in AIDA_data folder.
     calib_variables = data_to_calib_inputs(
         FT, batch_name, data_file_name_list, nuc_mode,
         start_time, end_time,
-        w, t_max, end_sim,
+        w, t_max,
+        P3_nuc_mode,
     )
     (;
         t_profile, T_profile, P_profile, ICNC_profile, e_profile, S_l_profile,
         params_list, IC_list,
         frozen_frac, frozen_frac_moving_mean, ICNC_moving_avg,
         Nₜ, y_truth,
+        P3_params_list,
     ) = calib_variables
 
     if batch_name in ["HOM", "DEP", "IMM"]
@@ -123,7 +128,7 @@ for (batch_index, batch_name) in enumerate(batch_names)
     Γ = (0.1 / 3)^2 * LinearAlgebra.I   # 0.1 is the uncertainty in observed ICNC
 
     ### Calibration.
-    UKI_output = calibrate_J_parameters_UKI(FT, nuc_mode, params_list, IC_list, y_truth, end_sim, Γ)
+    UKI_output = calibrate_J_parameters_UKI(FT, nuc_mode, params_list, IC_list, y_truth, Γ)
 
     UKI_calibrated_parameters = UKI_output[1] # MEAN of parameters from ensembles in FINAL iteration 
     UKI_all_params = UKI_output[2]            # parameters from EACH ensemble in EACH iteration
@@ -151,11 +156,14 @@ for (batch_index, batch_name) in enumerate(batch_names)
         end
 
         ## Calibrated parcel.
-        UKI_parcel = run_model([params_list[exp_index]], nuc_mode, UKI_calibrated_parameters, FT, [IC_list[exp_index]], end_sim)
+        UKI_parcel = run_model([params_list[exp_index]], nuc_mode, UKI_calibrated_parameters, FT, [IC_list[exp_index]])
+        @info(exp_names[exp_index])
+        P3_parcel = run_model([P3_params_list[exp_index]], P3_nuc_mode, [0, 0], FT, [IC_list[exp_index]], P3 = true)
         if batch_name in ["HOM", "DEP", "IMM"]
             # the order it appends is "in05_17_aida.edf", "in05_18_aida.edf", "TROPIC04",
             # "in07_01_aida.edf", "in07_19_aida.edf", "EXP45", "ACI04_22", "EXP19".
             push!(overview_data.UKI_calibrated_parcel, UKI_parcel)
+            push!(overview_data.P3_parcel, P3_parcel)
         end
 
         ## Plots.
@@ -192,12 +200,12 @@ for (batch_index, batch_name) in enumerate(batch_names)
             t_profile[exp_index], frozen_frac_moving_mean[exp_index], frozen_frac[exp_index],
         )
 
-        ## Looking at spread in UKI calibrated parameters
-        UKI_parcel_1 = run_model([params_list[exp_index]], nuc_mode, [UKI_final_iter_spread[1,1], UKI_final_iter_spread[2,1]], FT, [IC_list[exp_index]], end_sim)
-        UKI_parcel_2 = run_model([params_list[exp_index]], nuc_mode, [UKI_final_iter_spread[1,2], UKI_final_iter_spread[2,2]], FT, [IC_list[exp_index]], end_sim)
-        UKI_parcel_3 = run_model([params_list[exp_index]], nuc_mode, [UKI_final_iter_spread[1,3], UKI_final_iter_spread[2,3]], FT, [IC_list[exp_index]], end_sim)
-        UKI_parcel_4 = run_model([params_list[exp_index]], nuc_mode, [UKI_final_iter_spread[1,4], UKI_final_iter_spread[2,4]], FT, [IC_list[exp_index]], end_sim)
-        UKI_parcel_5 = run_model([params_list[exp_index]], nuc_mode, [UKI_final_iter_spread[1,5], UKI_final_iter_spread[2,5]], FT, [IC_list[exp_index]], end_sim)
+        # Looking at spread in UKI calibrated parameters
+        UKI_parcel_1 = run_model([params_list[exp_index]], nuc_mode, [UKI_final_iter_spread[1,1], UKI_final_iter_spread[2,1]], FT, [IC_list[exp_index]])
+        UKI_parcel_2 = run_model([params_list[exp_index]], nuc_mode, [UKI_final_iter_spread[1,2], UKI_final_iter_spread[2,2]], FT, [IC_list[exp_index]])
+        UKI_parcel_3 = run_model([params_list[exp_index]], nuc_mode, [UKI_final_iter_spread[1,3], UKI_final_iter_spread[2,3]], FT, [IC_list[exp_index]])
+        UKI_parcel_4 = run_model([params_list[exp_index]], nuc_mode, [UKI_final_iter_spread[1,4], UKI_final_iter_spread[2,4]], FT, [IC_list[exp_index]])
+        UKI_parcel_5 = run_model([params_list[exp_index]], nuc_mode, [UKI_final_iter_spread[1,5], UKI_final_iter_spread[2,5]], FT, [IC_list[exp_index]])
 
         UKI_spread_fig = plot_UKI_spread(
             exp_names[exp_index],

papers/ice_nucleation_2024/calibration.jl

@@ -16,7 +16,7 @@ import StatsBase as SB
 include(joinpath(pkgdir(CM), "parcel", "Parcel.jl"))
 
 # Define model which wraps around parcel and overwrites calibrated parameters
-function run_model(p_list, IN_mode, coefficients, FT, IC_list, end_sim::Int64; calibration = false)
+function run_model(p_list, IN_mode, coefficients, FT, IC_list; calibration = false, P3 = false)
     # grabbing calibrated m and c
     m_calibrated, c_calibrated = coefficients
 
@@ -33,7 +33,9 @@ function run_model(p_list, IN_mode, coefficients, FT, IC_list, end_sim::Int64; c
         IC = IC_list[exp_index]
 
         # overwriting
-        if aerosol == CMP.Kaolinite(FT)
+        if P3 == true
+            nothing
+        elseif aerosol == CMP.Kaolinite(FT)
             if IN_mode == "ABDINM"
                 override_file = Dict(
                     "China2017_J_deposition_m_Kaolinite" =>
@@ -194,7 +196,6 @@ function run_model(p_list, IN_mode, coefficients, FT, IC_list, end_sim::Int64; c
         local sol = run_parcel(IC, FT(0), t_max, params)
 
         if calibration == true
-            # loss_func_i = sol[9, (end - end_sim):end] ./ (IC[7] + IC[8] + IC[9])
             loss_func_i = sol[9, end] / (IC[7] + IC[8] + IC[9])
             append!(loss_func, loss_func_i)
         elseif calibration == false
@@ -292,7 +293,7 @@ function create_prior(FT, IN_mode, ; perfect_model = false, aerosol_type = nothi
     return prior
 end
 
-function calibrate_J_parameters_EKI(FT, IN_mode, params, IC, y_truth, end_sim, Γ,; perfect_model = false)
+function calibrate_J_parameters_EKI(FT, IN_mode, params, IC, y_truth, Γ,; perfect_model = false)
     @info("Starting EKI calibration")
     # Random number generator
     rng_seed = 24
@@ -326,7 +327,7 @@ function calibrate_J_parameters_EKI(FT, IN_mode, params, IC, y_truth, end_sim, 
         ϕ_n = EKP.get_ϕ_final(prior, EKI_obj)
         G_ens = hcat(
             [
-                run_model(params, IN_mode, ϕ_n[:, i], FT, IC, end_sim, calibration = true) for
+                run_model(params, IN_mode, ϕ_n[:, i], FT, IC, calibration = true) for
                 i in 1:N_ensemble
             ]...,
         )
@@ -362,7 +363,7 @@ function calibrate_J_parameters_EKI(FT, IN_mode, params, IC, y_truth, end_sim, 
     return [calibrated_coeffs, ϕ_n_values, mean_each_iter, error]
 end
 
-function calibrate_J_parameters_UKI(FT, IN_mode, params, IC, y_truth, end_sim, Γ,; perfect_model = false)
+function calibrate_J_parameters_UKI(FT, IN_mode, params, IC, y_truth, Γ,; perfect_model = false)
     @info("Starting UKI calibration")
     (; aerosol) = params[1]
 
@@ -393,7 +394,7 @@ function calibrate_J_parameters_UKI(FT, IN_mode, params, IC, y_truth, end_sim, 
         ϕ_n = EKP.get_ϕ_final(prior, UKI_obj)
         # Evaluate forward map
         G_n = [
-            run_model(params, IN_mode, ϕ_n[:, i], FT, IC, end_sim, calibration = true) for
+            run_model(params, IN_mode, ϕ_n[:, i], FT, IC, calibration = true) for
             i in 1:size(ϕ_n)[2]  #i in 1:N_ensemble
         ]
         # Reformat into `d x N_ens` matrix

papers/ice_nucleation_2024/calibration_setup.jl

@@ -138,7 +138,7 @@ function perf_model_IC(FT, IN_mode)
     return [Sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, FT(0)]
 end
 
-function perf_model_pseudo_data(FT, IN_mode, params, IC, end_sim)
+function perf_model_pseudo_data(FT, IN_mode, params, IC)
     n_samples = 10
 
     if IN_mode == "ABDINM"
@@ -149,7 +149,7 @@ function perf_model_pseudo_data(FT, IN_mode, params, IC, end_sim)
         coeff_true = [FT(255.927125), FT(-68.553283)]
     end
 
-    G_truth = run_model(params, IN_mode, coeff_true, FT, IC, end_sim, calibration = true)
+    G_truth = run_model(params, IN_mode, coeff_true, FT, IC, calibration = true)
     dim_output = length(G_truth)
 
     Γ = (0.1 / 3)^2 * LinearAlgebra.I
@@ -160,15 +160,15 @@ function perf_model_pseudo_data(FT, IN_mode, params, IC, end_sim)
     return [y_truth, Γ, coeff_true]
 end
 
-function AIDA_IN05_params(FT, w, t_max, t_profile, T_profile, P_profile)
+function AIDA_IN05_params(FT, w, t_max, t_profile, T_profile, P_profile; P3 = false)
     IN_mode = "ABHOM"
     const_dt = FT(1)
     prescribed_thermodynamics = true
     aerosol_act = "AeroAct"
     aerosol = CMP.Sulfate(FT)
-    dep_nucleation = "ABDINM"
-    heterogeneous = "ABIFM"
-    homogeneous = "ABHOM"
+    dep_nucleation = P3 ? "P3_dep" : "ABDINM"
+    heterogeneous = P3 ? "P3_het" : "ABIFM"
+    homogeneous = P3 ? "P3_hom" : "ABHOM"
     condensation_growth = "Condensation"
     deposition_growth = "Deposition"
     liq_size_distribution = "Gamma"
@@ -177,8 +177,9 @@ function AIDA_IN05_params(FT, w, t_max, t_profile, T_profile, P_profile)
     r_nuc = FT(1e-7) #FT(3.057e-6)
     A_aer = FT(4 * π * r_nuc^2)
     ips = CMP.IceNucleationParameters(FT)
+    aap = CMP.AerosolActivationParameters(FT)
 
-    params = (; const_dt, w, t_max, ips,
+    params = (; const_dt, w, t_max, ips, aap,
         prescribed_thermodynamics, t_profile, T_profile, P_profile,
         aerosol_act, aerosol, r_nuc, aero_σ_g, A_aer,   # aerosol activation
         condensation_growth, deposition_growth,         # growth
@@ -236,7 +237,7 @@ function AIDA_IN05_IC(FT, data_file)
     return [Sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, FT(0)]
 end
 
-function AIDA_IN07_params(FT, w, t_max, t_profile, T_profile, P_profile, batch_name)
+function AIDA_IN07_params(FT, w, t_max, t_profile, T_profile, P_profile, batch_name; P3 = false)
     IN_mode = "ABDINM"
     const_dt = FT(1)
     prescribed_thermodynamics = true
@@ -248,9 +249,9 @@ function AIDA_IN07_params(FT, w, t_max, t_profile, T_profile, P_profile, batch_n
     elseif batch_name == "DEP"
         aerosol = CMP.Dust(FT)
     end
-    dep_nucleation = "ABDINM"
-    heterogeneous = "ABIFM"
-    homogeneous = "ABHOM"
+    dep_nucleation = P3 ? "P3_dep" : "ABDINM"
+    heterogeneous = P3 ? "P3_het" : "ABIFM"
+    homogeneous = P3 ? "P3_hom" : "ABHOM"
     condensation_growth = "Condensation"
     deposition_growth = "Deposition"
     liq_size_distribution = "Gamma"
@@ -260,8 +261,9 @@ function AIDA_IN07_params(FT, w, t_max, t_profile, T_profile, P_profile, batch_n
     A_aer = FT(4 * π * r_nuc^2)
     # r_nuc = r₀ in IC
     ips = CMP.IceNucleationParameters(FT)
+    aap = CMP.AerosolActivationParameters(FT)
 
-    params = (; const_dt, w, t_max, ips,
+    params = (; const_dt, w, t_max, ips, aap,
         prescribed_thermodynamics, t_profile, T_profile, P_profile,
         aerosol_act, aerosol, r_nuc, aero_σ_g, A_aer,   # aerosol activation
         condensation_growth, deposition_growth,         # growth
@@ -320,15 +322,15 @@ function AIDA_IN07_IC(FT, data_file)
     return [Sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, FT(0)]
 end
 
-function TROPIC04_params(FT, w, t_max, t_profile, T_profile, P_profile)
+function TROPIC04_params(FT, w, t_max, t_profile, T_profile, P_profile; P3 = false)
     IN_mode = "ABHOM"
     const_dt = FT(1)
     prescribed_thermodynamics = true
     aerosol_act = "AeroAct"
     aerosol = CMP.Sulfate(FT)
-    dep_nucleation = "ABDINM"
-    heterogeneous = "ABIFM"
-    homogeneous = "ABHOM"
+    dep_nucleation = "None"
+    heterogeneous = "None"
+    homogeneous = P3 ? "P3_hom" : "ABHOM"
     condensation_growth = "Condensation"
     deposition_growth = "Deposition"
     liq_size_distribution = "Gamma"
@@ -337,8 +339,9 @@ function TROPIC04_params(FT, w, t_max, t_profile, T_profile, P_profile)
     r_nuc = FT(1.15 / 2 * 1e-6)
     A_aer = FT(4 * π * r_nuc^2)
     ips = CMP.IceNucleationParameters(FT)
+    aap = CMP.AerosolActivationParameters(FT)
 
-    params = (; const_dt, w, t_max, ips,
+    params = (; const_dt, w, t_max, ips, aap,
         prescribed_thermodynamics, t_profile, T_profile, P_profile,
         aerosol_act, aerosol, r_nuc, aero_σ_g, A_aer,   # aerosol activation
         condensation_growth, deposition_growth,         # growth
@@ -376,16 +379,16 @@ function TROPIC04_IC(FT)
     return [Sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, FT(0)]
 end
 
-function ACI04_22_params(FT, w, t_max, t_profile, T_profile, P_profile, batch_name)
+function ACI04_22_params(FT, w, t_max, t_profile, T_profile, P_profile, batch_name; P3 = false)
     # Niemand et al (2012)
     IN_mode = "ABIFM"
     const_dt = FT(1)
     prescribed_thermodynamics = true
     aerosol_act = "None"
     aerosol = batch_name == "ACI04_22" ? CMP.MiddleEasternDust(FT) : CMP.Dust(FT)
-    dep_nucleation = "ABDINM"
-    heterogeneous = "ABIFM"
-    homogeneous = "ABHOM"
+    dep_nucleation = P3 ? "P3_dep" : "ABDINM"
+    heterogeneous = P3 ? "P3_het" : "ABIFM"
+    homogeneous = P3 ? "P3_hom" : "ABHOM"
     condensation_growth = "Condensation"
     deposition_growth = "Deposition"
     liq_size_distribution = "Gamma"
@@ -394,8 +397,9 @@ function ACI04_22_params(FT, w, t_max, t_profile, T_profile, P_profile, batch_na
     r_nuc = FT(0.645 / 2 * 1e-6)  # avg of 2 modes
     A_aer = FT(4 * π * r_nuc^2)
     ips = CMP.IceNucleationParameters(FT)
+    aap = CMP.AerosolActivationParameters(FT)
 
-    params = (; const_dt, w, t_max, ips,
+    params = (; const_dt, w, t_max, ips, aap,
         prescribed_thermodynamics, t_profile, T_profile, P_profile,
         aerosol_act, aerosol, r_nuc, aero_σ_g, A_aer,   # aerosol activation
         condensation_growth, deposition_growth,         # growth
@@ -433,16 +437,16 @@ function ACI04_22_IC(FT)
     return [Sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, FT(0)]
 end
 
-function EXP19_params(FT, w, t_max, t_profile, T_profile, P_profile, batch_name)
+function EXP19_params(FT, w, t_max, t_profile, T_profile, P_profile, batch_name; P3 = false)
     # Cotten et al (2007)
     IN_mode = "ABIFM"
     const_dt = FT(1)
     prescribed_thermodynamics = true
     aerosol_act = "None"
     aerosol = batch_name == "EXP19" ? CMP.AsianDust(FT) : CMP.Dust(FT)
-    dep_nucleation = "ABDINM"
-    heterogeneous = "ABIFM"
-    homogeneous = "ABHOM"
+    dep_nucleation = P3 ? "P3_dep" : "ABDINM"
+    heterogeneous = P3 ? "P3_het" : "ABIFM"
+    homogeneous = P3 ? "P3_hom" : "ABHOM"
     condensation_growth = "Condensation"
     deposition_growth = "Deposition"
     liq_size_distribution = "Gamma"
@@ -451,8 +455,9 @@ function EXP19_params(FT, w, t_max, t_profile, T_profile, P_profile, batch_name)
     r_nuc = FT(0.4 / 2 * 1e-6)   # value is mode radius, not mean
     A_aer = FT(4 * π * r_nuc^2)
     ips = CMP.IceNucleationParameters(FT)
+    aap = CMP.AerosolActivationParameters(FT)
 
-    params = (; const_dt, w, t_max, ips,
+    params = (; const_dt, w, t_max, ips, aap,
         prescribed_thermodynamics, t_profile, T_profile, P_profile,
         aerosol_act, aerosol, r_nuc, aero_σ_g, A_aer,   # aerosol activation
         condensation_growth, deposition_growth,         # growth
@@ -490,16 +495,16 @@ function EXP19_IC(FT)
     return [Sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, FT(0)]
 end
 
-function EXP45_params(FT, w, t_max, t_profile, T_profile, P_profile, batch_name)
+function EXP45_params(FT, w, t_max, t_profile, T_profile, P_profile, batch_name; P3 = false)
     # Cotten et al (2007)
     IN_mode = "ABDINM"
     const_dt = FT(1)
     prescribed_thermodynamics = true
     aerosol_act = "None"
     aerosol = batch_name == "EXP45" ? CMP.SaharanDust(FT) : CMP.Dust(FT)
-    dep_nucleation = "ABDINM"
-    heterogeneous = "ABIFM"
-    homogeneous = "ABHOM"
+    dep_nucleation = P3 ? "P3_dep" : "ABDINM"
+    heterogeneous = P3 ? "P3_het" : "ABIFM"
+    homogeneous = P3 ? "P3_hom" : "ABHOM"
     condensation_growth = "Condensation"
     deposition_growth = "Deposition"
     liq_size_distribution = "Gamma"
@@ -508,8 +513,9 @@ function EXP45_params(FT, w, t_max, t_profile, T_profile, P_profile, batch_name)
     r_nuc = FT(0.4 / 2 * 1e-6)   # value is mode radius, not mean
     A_aer = FT(4 * π * r_nuc^2)
     ips = CMP.IceNucleationParameters(FT)
+    aap = CMP.AerosolActivationParameters(FT)
 
-    params = (; const_dt, w, t_max, ips,
+    params = (; const_dt, w, t_max, ips, aap,
         prescribed_thermodynamics, t_profile, T_profile, P_profile,
         aerosol_act, aerosol, r_nuc, aero_σ_g, A_aer,   # aerosol activation
         condensation_growth, deposition_growth,         # growth