Fix perfect model ice nucleation calibs

papers/ice_nucleation_2024/AIDA_calibrations.jl

@@ -175,15 +175,15 @@ for (exp_index, data_file_name) in enumerate(data_file_names)
     # MK.lines!(
     #     ax_compare,
     #     parcel_default.t .+ (moving_average_n / 2),
-    #     parcel_default[9, :]./  (IC[7] + IC[8] + IC[9]),       # IC[8] if starting with droplets (no aerosol activation)
+    #     parcel_default[9, :]./  (IC[7] + IC[8] + IC[9]),
     #     label = "CM.jl Parcel (Pre-Calib)",
     #     linewidth = 2.5,
     #     color =:orangered,
     # )
     MK.lines!(
         ax_compare,
         parcel.t .+ (moving_average_n / 2),
-        parcel[9, :]./ (IC[7] + IC[8] + IC[9]),               # IC[8] if starting with droplets (no aerosol activation)
+        parcel[9, :]./ (IC[7] + IC[8] + IC[9]),
         label = "CM.jl Parcel (Calibrated)",
         linewidth = 2.5,
         color =:orange,

papers/ice_nucleation_2024/calibration.jl

@@ -47,10 +47,6 @@ function run_model(p, coefficients, IN_mode, FT, IC)
             ice_size_distribution = ice_size_distribution,
         )
 
-        # solve ODE
-        local sol = run_parcel(IC, FT(0), t_max, params)
-        return sol[9, :] .* 1e10 # ICNC magnified
-
     elseif IN_mode == "ABIFM"
         # overwriting
         override_file = Dict(
@@ -74,10 +70,6 @@ function run_model(p, coefficients, IN_mode, FT, IC)
             ice_size_distribution = ice_size_distribution,
         )
 
-        # solve ODE
-        local sol = run_parcel(IC, FT(0), t_max, params)
-        return sol[9, :] # ICNC
-
     elseif IN_mode == "ABHOM"
         # overwriting
         override_file = Dict(
@@ -104,22 +96,76 @@ function run_model(p, coefficients, IN_mode, FT, IC)
             ips = overwrite,
             r_nuc = r_nuc,
         )
-
-        # solve ODE
-        local sol = run_parcel(IC, FT(0), t_max, params)
-        return sol[9, :] ./  (IC[7] + IC[8] + IC[9])  # frozen fraction
     end
+
+    # solve ODE
+    local sol = run_parcel(IC, FT(0), t_max, params)
+    return sol[9, :] ./ (IC[7] + IC[8] + IC[9])  # frozen fraction
 end
 
 function run_calibrated_model(FT, IN_mode, coefficients, p, IC)
     # grabbing parameters
     m_calibrated, c_calibrated = coefficients
-    (; const_dt, w, t_max, aerosol_act, aerosol) = p
-    (; deposition_growth, condensation_growth) = p
-    (; homogeneous) = p
+    (; const_dt, w, t_max, aerosol_act, aerosol, r_nuc) = p
     (; liq_size_distribution, ice_size_distribution, aero_σ_g) = p
+    (; deposition_growth, condensation_growth) = p
+
+    if IN_mode == "ABDINM"
+        (; dep_nucleation) = p
+        # 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)
+        overwrite = CMP.Kaolinite(kaolinite_calibrated)
+
+        # run parcel with new coefficients
+        local params = parcel_params{FT}(
+            const_dt = const_dt,
+            w = w,
+            aerosol_act = aerosol_act,
+            aerosol = overwrite,
+            aero_σ_g = aero_σ_g,
+            deposition = dep_nucleation,
+            condensation_growth = condensation_growth,
+            deposition_growth = deposition_growth,
+            liq_size_distribution = liq_size_distribution,
+            ice_size_distribution = ice_size_distribution,
+            ips = overwrite,
+        )
+
+    elseif IN_mode == "ABIFM"
+        (; heterogeneous) = p
+        # overwriting
+        override_file = Dict(
+            "KnopfAlpert2013_J_ABIFM_m_Kaolinite" =>
+                Dict("value" => m_calibrated, "type" => "float"),
+            "KnopfAlpert2013_J_ABIFM_c_Kaolinite" =>
+                Dict("value" => c_calibrated, "type" => "float"),
+        )
+        kaolinite_calibrated = CP.create_toml_dict(FT; override_file)
+        overwrite = CMP.Kaolinite(kaolinite_calibrated)
+
+        # run parcel with new coefficients
+        local params = parcel_params{FT}(
+            const_dt = const_dt,
+            w = w,
+            aerosol_act = aerosol_act,
+            aerosol = overwrite,
+            aero_σ_g = aero_σ_g,
+            heterogeneous = heterogeneous,
+            condensation_growth = condensation_growth,
+            deposition_growth = deposition_growth,
+            liq_size_distribution = liq_size_distribution,
+            ice_size_distribution = ice_size_distribution,
+            ips = overwrite,
+        )
 
-    if IN_mode == "ABHOM"
+    elseif IN_mode == "ABHOM"
+        (; homogeneous) = p
         # overwriting
         override_file = Dict(
             "Linear_J_hom_coeff2" =>
@@ -143,12 +189,12 @@ function run_calibrated_model(FT, IN_mode, coefficients, p, IC)
             liq_size_distribution = liq_size_distribution,
             ice_size_distribution = ice_size_distribution,
             ips = overwrite,
+            r_nuc = r_nuc,
         )
-
-        # solve ODE
-        local sol = run_parcel(IC, FT(0), t_max, params)
-        return sol
     end
+    # solve ODE
+    local sol = run_parcel(IC, FT(0), t_max, params)
+    return sol
 end
 
 function create_prior(FT, IN_mode, ; perfect_model = false)
@@ -159,14 +205,14 @@ function create_prior(FT, IN_mode, ; perfect_model = false)
     # for perfect model calibration
     if perfect_model == true
         if IN_mode == "ABDINM"
-            m_stats = [FT(20), FT(1), FT(0), Inf]
-            c_stats = [FT(-1), FT(1), -Inf, Inf]
+            m_stats = [FT(20), FT(8), FT(0), Inf]
+            c_stats = [FT(-1), FT(2), -Inf, Inf]
         elseif IN_mode == "ABIFM"
-            m_stats = [FT(50), FT(1), FT(0), Inf]
-            c_stats = [FT(-7), FT(1), -Inf, Inf]
+            m_stats = [FT(50), FT(7), FT(0), Inf]
+            c_stats = [FT(-7), FT(4), -Inf, Inf]
         elseif IN_mode == "ABHOM"
-            m_stats = [FT(260), FT(1), FT(0), Inf]
-            c_stats = [FT(-70), FT(1), -Inf, Inf]
+            m_stats = [FT(251), FT(6), FT(0), Inf]
+            c_stats = [FT(-66), FT(3), -Inf, Inf]
         end
     elseif perfect_model == false
         if IN_mode == "ABDINM"
@@ -207,8 +253,8 @@ function calibrate_J_parameters(FT, IN_mode, params, IC, y_truth, Γ,; perfect_m
     rng = Random.seed!(Random.GLOBAL_RNG, rng_seed)
 
     prior = create_prior(FT, IN_mode, perfect_model = perfect_model)
-    N_ensemble = 10      # runs N_ensemble trials per iteration
-    N_iterations = 150    # number of iterations the inverse problem goes through
+    N_ensemble = 25      # runs N_ensemble trials per iteration
+    N_iterations = 100    # number of iterations the inverse problem goes through
 
     # Generate initial ensemble and set up EKI
     initial_ensemble = EKP.construct_initial_ensemble(rng, prior, N_ensemble)

papers/ice_nucleation_2024/calibration_setup.jl

@@ -7,10 +7,10 @@ import Thermodynamics as TD
 function perf_model_params(FT, IN_mode)
     if IN_mode == "ABDINM"
         const_dt = FT(1)
-        w = FT(0.4)
-        t_max = FT(100)
+        w = FT(5)
+        t_max = FT(300)
         aerosol_act = "None"
-        aerosol = "None"
+        aerosol = CMP.Kaolinite(FT)
         aero_σ_g = FT(0)
         r_nuc = FT(1e-7)
         dep_nucleation = "ABDINM"
@@ -22,10 +22,10 @@ function perf_model_params(FT, IN_mode)
         ice_size_distribution = "Monodisperse"
     elseif IN_mode == "ABIFM"
         const_dt = FT(1)
-        w = FT(0.4)
+        w = FT(1)
         t_max = FT(100)
         aerosol_act = "None"
-        aerosol = "None"
+        aerosol = CMP.Kaolinite(FT)
         aero_σ_g = FT(0)
         r_nuc = FT(1e-7)
         dep_nucleation = "None"
@@ -70,14 +70,13 @@ function perf_model_IC(FT, IN_mode)
     ϵₘ = R_d / R_v
 
     if IN_mode == "ABDINM"
-        Nₐ = FT(8e6) # according to fig 2 panel b deposition nucleation experiment
+        Nₐ = FT(2e8)
         Nₗ = FT(0)
         Nᵢ = FT(0)
-        r₀ = FT(0.5e-6)
-        p₀ = FT(987.018 * 1e2)
-        T₀ = FT(212.978)
-        C_v = FT(10.8509 * 1e-6)
-        qᵥ = ϵₘ / (ϵₘ - 1 + 1 / C_v)
+        r₀ = FT(1e-7)
+        p₀ = FT(80000)
+        T₀ = FT(235)
+        qᵥ = FT(8.8e-5)
         qₗ = FT(0)
         qᵢ = FT(0)
         q = TD.PhasePartition.(qᵥ + qₗ + qᵢ, qₗ, qᵢ)
@@ -87,15 +86,13 @@ function perf_model_IC(FT, IN_mode)
         Sₗ = FT(e / eₛ)
     elseif IN_mode == "ABIFM"
         Nₐ = FT(0)
-        Nₗ = FT(8e6)
+        Nₗ = FT(2000)
         Nᵢ = FT(0)
-        r₀ = FT(0.5e-6)
-        p₀ = FT(987.018 * 1e2)
-        T₀ = FT(212.978)
+        r₀ = FT(1e-6)
+        p₀ = FT(800 * 1e2)
+        T₀ = FT(251)
         qₗ = FT(Nₗ * 4 / 3 * FT(π) * r₀^3 * ρₗ / FT(1.2)) # 1.2 should be ρₐ
-        C_l = FT(qₗ / ((1 - qₗ) * ϵₘ + qₗ))  # concentration/mol fraction of liquid
-        C_v = FT(10.8509 * 1e-6 - C_l)     # concentration/mol fraction of vapor
-        qᵥ = ϵₘ / (ϵₘ - 1 + 1 / C_v)
+        qᵥ = FT(8.1e-4)
         qᵢ = FT(0)
         q = TD.PhasePartition.(qᵥ + qₗ + qᵢ, qₗ, qᵢ)
         Rₐ = TD.gas_constant_air(tps, q)
@@ -110,7 +107,6 @@ function perf_model_IC(FT, IN_mode)
         p₀ = FT(9712.183)
         T₀ = FT(190)
         qₗ = FT(Nₗ * 4 / 3 * FT(π) * r₀^3 * ρₗ / FT(1.2)) # 1.2 should be ρₐ
-        C_l = FT(qₗ / ((1 - qₗ) * ϵₘ + qₗ))  # concentration/mol fraction of liquid
         C_v = FT(5 * 1e-6)
         qᵥ = ϵₘ / (ϵₘ - 1 + 1 / C_v)
         qᵢ = FT(0)
@@ -137,7 +133,7 @@ function perf_model_pseudo_data(FT, IN_mode, params, IC)
     G_truth = run_model(params, coeff_true, IN_mode, FT, IC)
     dim_output = length(G_truth)
 
-    Γ = 0.005 * LinearAlgebra.I * (maximum(G_truth) - minimum(G_truth))
+    Γ = 0.03 * LinearAlgebra.I * (maximum(G_truth) - minimum(G_truth))
     noise_dist = Distributions.MvNormal(zeros(dim_output), Γ)
 
     y_truth = zeros(length(G_truth), n_samples) # where noisy ICNC will be stored

papers/ice_nucleation_2024/perfect_calib.jl

@@ -66,19 +66,16 @@ for IN_mode in IN_mode_list
     plot_name = "perfect_calibration_$mode_label.svg"
     MK.save(plot_name, fig)
 
-    # Plotting calibrated parcel's ICNC if ABHOM
-    if IN_mode == "ABHOM"
-        fig2 = MK.Figure(size = (800, 600))
-        ax3 = MK.Axis(fig2[1, 1], ylabel = "ICNC [m^3]", xlabel = "time [s]", title = IN_mode)
-
-        soln = run_calibrated_model(FT, IN_mode, calibrated_parameters, params, IC)
-        soln_dflt = run_calibrated_model(FT, IN_mode, [FT(255.927125), FT(-68.553283)], params, IC)
-
-        MK.lines!(ax3, soln.t, soln[9,:], label = "calibrated")
-        MK.lines!(ax3, soln_dflt.t, soln_dflt[9,:], label = "default", color = :orange)
-        plot_name = "perfect_calibration_ICNC_$mode_label.svg"
-        MK.save(plot_name, fig2)
-    end
+    # Plotting calibrated parcel's ICNC
+    fig2 = MK.Figure(size = (800, 600))
+    ax3 = MK.Axis(fig2[1, 1], ylabel = "Frozen Fraction [-]", xlabel = "time [s]", title = IN_mode)
+
+    soln = run_calibrated_model(FT, IN_mode, calibrated_parameters, params, IC)
+    soln_dflt = run_calibrated_model(FT, IN_mode, coeff_true, params, IC)
 
+    MK.lines!(ax3, soln.t, soln[9,:]./ (IC[7] + IC[8] + IC[9]), label = "calibrated")
+    MK.lines!(ax3, soln_dflt.t, soln_dflt[9,:]./ (IC[7] + IC[8] + IC[9]), label = "default", color = :orange)
+    plot_name = "perfect_calibration_ICNC_$mode_label.svg"
+    MK.save(plot_name, fig2)
 end
 #! format: on

parcel/ParcelModel.jl

@@ -299,7 +299,8 @@ function run_parcel(IC, t_0, t_end, pp)
     else
         throw("Unrecognized deposition growth mode")
     end
-    @info info
+    # To see parcel parameters in output, uncomment
+    # @info info
 
     # Parameters for the ODE solver
     p = (

test/ice_nucleation_calibration.jl

@@ -24,9 +24,26 @@ function test_J_calibration(FT, IN_mode)
         perfect_model = true,
     )
     calibrated_parameters = [output[1], output[2]]
-
-    TT.@test calibrated_parameters[1] ≈ coeff_true[1] rtol = FT(0.2)
-    TT.@test calibrated_parameters[2] ≈ coeff_true[2] rtol = FT(0.2)
+    calibrated_soln =
+        run_calibrated_model(FT, IN_mode, calibrated_parameters, params, IC)
+    true_soln = run_calibrated_model(FT, IN_mode, coeff_true, params, IC)
+
+
+    # test that coeffs are close to "true" values
+    if IN_mode == "ABDINM"
+        TT.@test calibrated_parameters[1] ≈ coeff_true[1] rtol = FT(0.3)
+        TT.@test calibrated_parameters[2] ≈ coeff_true[2] atol = FT(3)
+    elseif IN_mode == "ABIFM"
+        TT.@test calibrated_parameters[1] ≈ coeff_true[1] rtol = FT(0.3)
+        TT.@test calibrated_parameters[2] ≈ coeff_true[2] atol = FT(7)
+    elseif IN_mode == "ABHOM"
+        TT.@test calibrated_parameters[1] ≈ coeff_true[1] rtol = FT(0.3)
+        TT.@test calibrated_parameters[2] ≈ coeff_true[2] atol = FT(20)
+    end
+
+    # test that resulting ICNC are similar
+    TT.@test (calibrated_soln[9, end] ./ (IC[7] + IC[8] + IC[9])) ≈
+             (true_soln[9, end] ./ (IC[7] + IC[8] + IC[9])) rtol = FT(0.1)
 end
 
 @info "Ice Nucleation Calibration Test"

test/performance_tests.jl

@@ -266,7 +266,7 @@ function benchmark_test(FT)
         (liquid, rain, blk1mvel.rain, ce, q_liq, q_rai, ρ_air),
         350,
     )
-    bench_press(CM1.radar_reflectivity, (rain, q_rai, ρ_air), 250)
+    bench_press(CM1.radar_reflectivity, (rain, q_rai, ρ_air), 300)
 
     # 2-moment
     for sb in [sb2006, sb2006_no_limiters]
@@ -293,7 +293,7 @@ function benchmark_test(FT)
         bench_press(
             CM2.rain_terminal_velocity,
             (sb, ch2022.rain, q_rai, ρ_air, N_rai),
-            2000,
+            2200,
         )
         bench_press(
             CM2.radar_reflectivity,