calibrating AIDA HET experiments

.buildkite/pipeline.yml

@@ -40,3 +40,12 @@ steps:
       CLIMACOMMS_DEVICE: "CUDA"
     agents:
       slurm_gpus: 1
+
+  - label: ":cyclone: ClimaCore + GPU unit tests"
+    key: "clima_core_gpu_unittests"
+    command:
+      - "julia --project=test --color=yes test/gpu_clima_core_test.jl"
+    env:
+      CLIMACOMMS_DEVICE: "CUDA"
+    agents:
+      slurm_gpus: 1

papers/ice_nucleation_2024/AIDA_calibrations.jl

@@ -33,11 +33,11 @@ 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"]
-plot_names = ["IN0517", "IN0518", "IN0519"]
+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"]
 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
+end_time_list = [290, 290, 170, 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
 
@@ -49,6 +49,7 @@ R_d = TD.Parameters.R_d(tps)
 
 
 for (exp_index, data_file_name) in enumerate(data_file_names)
+    @info(data_file_name)
     #! format: off
     ### Unpacking experiment-specific variables.
     plot_name = plot_names[exp_index]
@@ -58,6 +59,8 @@ 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" ? "ABDINM" : "ABHOM"
+
     ## 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))
@@ -72,10 +75,16 @@ for (exp_index, data_file_name) in enumerate(data_file_names)
     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]
-    # S_l_profile = e_profile ./ (TD.saturation_vapor_pressure.(tps, T_profile, TD.Liquid()))
+    S_l_profile = e_profile ./ (TD.saturation_vapor_pressure.(tps, T_profile, TD.Liquid()))
 
-    params = AIDA_IN05_params(FT, w, t_max, t_profile, T_profile, P_profile)
-    IC = AIDA_IN05_IC(FT, data_file_name)
+    params =
+        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 == "ABHOM" ?
+        AIDA_IN05_IC(FT, data_file_name) :
+        AIDA_IN07_IC(FT, data_file_name)
 
     Nₜ = IC[7] + IC[8] + IC[9]
     frozen_frac = AIDA_ICNC[start_time_index:end_time_index] ./ Nₜ
@@ -84,8 +93,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]
@@ -95,16 +104,16 @@ 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.
     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, 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.lines!(ax1, t_profile .+ moving_average_start_index .- 100, 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)
 
@@ -131,30 +140,30 @@ for (exp_index, data_file_name) in enumerate(data_file_names)
 
     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, 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, S_l_profile, label = "chamber", color = :blue)
+    MK.lines!(ax_parcel_1, t_profile, S_l_profile, label = "chamber", color = :blue)
 
     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_2, parcel_default.t, parcel_default[5,:], 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, 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_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_4, parcel_default.t, parcel_default[6, :], color = :darkorange2)
 
     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_5, parcel_default.t, parcel_default[8, :], color = :darkorange2)
 
     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, 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))
@@ -223,11 +232,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/AIDA_data_artifact/README.md

@@ -12,16 +12,19 @@ After logging in, data is found under "Data Curation of Homeless Data" -->
 bring the user to the EUROCHAMP2020 interface where AIDA chamber data can be found.
 
 ## Data details
-The point of contact and author of the raw homogeneous ice
-nucleation data (IN05_17, IN05_18, and IN05_19) is Ottmar Mohler.
-This data includes a time series of pressure, mean gas temperature, water partial
-pressure, total water concentration, ice number density, droplet number density,
-and droplet diameter. These temperature-dependent rate measurements were made
-during adiabatic expansion experiments in the AIDA chamber.
+The point of contact and author of the raw ice nucleation data
+from the AIDA chamber is Ottmar Mohler. The homogeneous ice nucleation
+data (IN05_17, IN05_18, IN05_19) includes a time series of pressure,
+mean gas temperature, water partial pressure, total water concentration,
+ice number density from WELAS, droplet number density, and droplet diameter.
+The heterogeneous ice nucleation data (IN07_01, IN07_19) includes the same
+first six variables, ice number density from FTIR, and ice water content.
+These measurements were made during adiabatic expansion experiments in the AIDA chamber.
 Experimental details can be found in Benz et al (2005) here:
-https://www.sciencedirect.com/science/article/pii/S1010603005004181. 
-
-More AIDA chamebr data is to be added from various heterogeneous experiments.
+https://www.sciencedirect.com/science/article/pii/S1010603005004181.
+More details on the heterogeneous ice nucleation experiments can be found
+in Mohler et al (2008) here:
+https://iopscience.iop.org/article/10.1088/1748-9326/3/2/025007/pdf.
 
 ## Licensing
 License: Creative Commons Attribution 4.0 International (CC by 4.0)

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"]
 
@@ -178,13 +218,17 @@ function create_prior(FT, IN_mode, ; perfect_model = false)
         end
     elseif perfect_model == false
         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(40), FT(20), FT(0), Inf]
+                c_stats = [FT(0.5), FT(20), -Inf, Inf]
+            elseif aerosol_type == CMP.Illite(FT)
+                m_stats = [FT(30), 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]
@@ -210,11 +254,14 @@ function create_prior(FT, IN_mode, ; perfect_model = false)
 end
 
 function calibrate_J_parameters_EKI(FT, IN_mode, params, IC, y_truth, end_sim, Γ,; perfect_model = false)
+    @info("Starting EKI calibration")
     # Random number generator
     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
 
@@ -233,7 +280,7 @@ function calibrate_J_parameters_EKI(FT, IN_mode, params, IC, y_truth, end_sim, 
     # Carry out the EKI calibration
     # ϕ_n_values[iteration] stores ensembles of calibrated coeffs in that iteration
     global ϕ_n_values = []
-    final_iter = N_iterations
+    global final_iter = N_iterations
     for n in 1:N_iterations
         ϕ_n = EKP.get_ϕ_final(prior, EKI_obj)
         G_ens = hcat(
@@ -269,7 +316,10 @@ 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)
+    @info("Starting UKI calibration")
+    (; aerosol) = params
+
+    prior = create_prior(FT, IN_mode, perfect_model = perfect_model, aerosol_type = aerosol)
     N_iterations = 25
     α_reg = 1.0
     update_freq = 1
@@ -289,8 +339,8 @@ function calibrate_J_parameters_UKI(FT, IN_mode, params, IC, y_truth, end_sim, 
     )
     UKI_obj = EKP.EnsembleKalmanProcess(truth, process; verbose = true)
 
-    err = []
-    final_iter =[N_iterations]
+    global err = []
+    global final_iter =[N_iterations]
     for n in 1:N_iterations
         # Return transformed parameters in physical/constrained space
         ϕ_n = EKP.get_ϕ_final(prior, UKI_obj)
@@ -304,14 +354,14 @@ function calibrate_J_parameters_UKI(FT, IN_mode, params, IC, y_truth, end_sim, 
         # Update ensemble
         terminate = EKP.EnsembleKalmanProcesses.update_ensemble!(UKI_obj, G_ens)
         push!(err, EKP.get_error(UKI_obj)[end])
-        println(
-            "Iteration: " *
-            string(n) *
-            ", Error: " *
-            string(err[n]) *
-            " norm(Cov):" *
-            string(Distributions.norm(EKP.get_process(UKI_obj).uu_cov[n]))
-        )
+        # println(
+        #     "Iteration: " *
+        #     string(n) *
+        #     ", Error: " *
+        #     string(err[n]) *
+        #     " norm(Cov):" *
+        #     string(Distributions.norm(EKP.get_process(UKI_obj).uu_cov[n]))
+        # )
         if !isnothing(terminate)
             final_iter[1] = n - 1
             break