going insane

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")

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

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)

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