From 20f85ea797f3d7941852da0d08614f861619b8d6 Mon Sep 17 00:00:00 2001 From: amylu00 Date: Tue, 17 Sep 2024 15:03:20 -0700 Subject: [PATCH] cleaning before the great squash --- papers/ice_nucleation_2024/AIDA_calibrations.jl | 12 +++++------- papers/ice_nucleation_2024/calibration.jl | 13 +------------ papers/ice_nucleation_2024/calibration_setup.jl | 8 ++------ 3 files changed, 8 insertions(+), 25 deletions(-) diff --git a/papers/ice_nucleation_2024/AIDA_calibrations.jl b/papers/ice_nucleation_2024/AIDA_calibrations.jl index 3f4cff200..e35f84bc0 100644 --- a/papers/ice_nucleation_2024/AIDA_calibrations.jl +++ b/papers/ice_nucleation_2024/AIDA_calibrations.jl @@ -54,7 +54,7 @@ for (exp_index, data_file_name) in enumerate(data_file_names) ### Check for data file in AIDA_data folder. chamber_data = unpack_data(data_file_name) - ICNC = chamber_data[start_time:end_time, 6] .* 1e6 # Nᵢ [m^-3] ; freezing onset to when it starts decreasing [195:390] + ICNC = chamber_data[start_time:end_time, 6] .* 1e6 # Nᵢ [m^-3] frozen_frac = ICNC ./ (IC[7] + IC[8] + IC[9]) frozen_frac_moving_mean = moving_average(frozen_frac, moving_average_n) @@ -63,7 +63,7 @@ for (exp_index, data_file_name) in enumerate(data_file_names) ### Plotting AIDA data. AIDA_data_fig = MK.Figure(size = (800, 600), fontsize = 24) - ax1 = MK.Axis(AIDA_data_fig[1, 1], ylabel = "Frozen Fraction [-]", xlabel = "time [s]", title = "AIDA data $plot_name") + ax1 = MK.Axis(AIDA_data_fig[1, 1], ylabel = "ICNC [m^-3]", xlabel = "time [s]", title = "AIDA data $plot_name") MK.lines!( ax1, chamber_data[start_time:end_time, 1], @@ -98,9 +98,8 @@ for (exp_index, data_file_name) in enumerate(data_file_names) ## Calibrated coefficients. # Did they converge? calibrated_coeffs_fig = MK.Figure(size = (600, 600), fontsize = 24) - ax3 = MK.Axis(calibrated_coeffs_fig[1, 1], ylabel = "m coefficient [-]")#, title = "$plot_name") + ax3 = MK.Axis(calibrated_coeffs_fig[1, 1], ylabel = "m coefficient [-]", title = "$plot_name") ax4 = MK.Axis(calibrated_coeffs_fig[2, 1], ylabel = "c coefficient [-]", xlabel = "iteration #") - #MK.ylims!(ax3, 250, 260) MK.lines!(ax3, collect(1:n_iterations), calibrated_ensemble_means[1], label = "ensemble mean", color = :orange, linewidth = 2.5) MK.lines!(ax4, collect(1:n_iterations), calibrated_ensemble_means[2], label = "ensemble mean", color = :orange, linewidth = 2.5) MK.save("$plot_name"*"_calibrated_coeffs_fig.svg", calibrated_coeffs_fig) @@ -173,11 +172,10 @@ for (exp_index, data_file_name) in enumerate(data_file_names) # Does calibrated parcel look like observations? ICNC_comparison_fig = MK.Figure(size = (700, 600), fontsize = 24) ax_compare = MK.Axis(ICNC_comparison_fig[1, 1], ylabel = "Frozen Fraction [-]", xlabel = "time [s]", title = "$plot_name") - #MK.ylims!(ax_compare, 0, 0.14) # 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, @@ -185,7 +183,7 @@ for (exp_index, data_file_name) in enumerate(data_file_names) 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, diff --git a/papers/ice_nucleation_2024/calibration.jl b/papers/ice_nucleation_2024/calibration.jl index 21b15baf1..6ffe268c0 100644 --- a/papers/ice_nucleation_2024/calibration.jl +++ b/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, :] ./ (IC[7] + IC[8] + IC[9]) # frozen fraction - 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, :] ./ (IC[7] + IC[8] + IC[9]) # frozen fraction - elseif IN_mode == "ABHOM" # overwriting override_file = Dict( @@ -90,7 +82,6 @@ function run_model(p, coefficients, IN_mode, FT, IC) overwrite = CMP.IceNucleationParameters(ip_calibrated) # run parcel with new coefficients - # TODO - need better way to call σ_g. Rename σ_g to σ_g_liq or something. local params = parcel_params{FT}( const_dt = const_dt, w = w, @@ -105,12 +96,10 @@ function run_model(p, coefficients, IN_mode, FT, IC) ips = overwrite, r_nuc = r_nuc, ) - - + 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 end function run_calibrated_model(FT, IN_mode, coefficients, p, IC) diff --git a/papers/ice_nucleation_2024/calibration_setup.jl b/papers/ice_nucleation_2024/calibration_setup.jl index ac73e5ff4..9f32dd599 100644 --- a/papers/ice_nucleation_2024/calibration_setup.jl +++ b/papers/ice_nucleation_2024/calibration_setup.jl @@ -75,7 +75,7 @@ function perf_model_IC(FT, IN_mode) ϵₘ = R_d / R_v if IN_mode == "ABDINM" - Nₐ = FT(2e8) # according to fig 2 panel b deposition nucleation experiment + Nₐ = FT(2e8) Nₗ = FT(0) Nᵢ = FT(0) r₀ = FT(1e-7) @@ -89,9 +89,6 @@ function perf_model_IC(FT, IN_mode) eₛ = TD.saturation_vapor_pressure(tps, T₀, TD.Liquid()) e = eᵥ(qᵥ, p₀, Rₐ, R_v) Sₗ = FT(e / eₛ) - Sᵢ = S_i(tps, T₀, Sₗ) - # @info("", Sᵢ) - # error("breaking") elseif IN_mode == "ABIFM" Nₐ = FT(0) Nₗ = FT(2000) @@ -115,7 +112,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) @@ -146,7 +142,7 @@ function perf_model_pseudo_data(FT, IN_mode, params, IC) noise_dist = Distributions.MvNormal(zeros(dim_output), Γ) y_truth = zeros(length(G_truth), n_samples) # where noisy ICNC will be stored - y_truth = (G_truth .+ rand(noise_dist)) + y_truth = G_truth .+ rand(noise_dist) return [y_truth, Γ, coeff_true] end