From c72834c0ce0b66c50211b11dc540e77ba490f2c2 Mon Sep 17 00:00:00 2001 From: amylu00 Date: Mon, 22 Apr 2024 16:44:22 -0700 Subject: [PATCH] Generalize ice nucleation calibration funcs --- NEWS.md | 1 + .../{perfect_model_J.jl => calibration.jl} | 183 +++--------------- .../ice_nucleation_2024/calibration_setup.jl | 128 ++++++++++++ .../ice_nucleation_2024/plot_ensemble_mean.jl | 118 ++++++----- test/ice_nucleation_calibration.jl | 27 ++- test/performance_tests.jl | 2 +- 6 files changed, 237 insertions(+), 222 deletions(-) rename papers/ice_nucleation_2024/{perfect_model_J.jl => calibration.jl} (53%) create mode 100644 papers/ice_nucleation_2024/calibration_setup.jl diff --git a/NEWS.md b/NEWS.md index 085836086..109fd7e4a 100644 --- a/NEWS.md +++ b/NEWS.md @@ -4,6 +4,7 @@ CloudMicrophysics.jl Release Notes main ------ +- Generalize calibration functions in ice_nucleation_2024 ([#380](https://github.com/CliMA/CloudMicrophysics.jl/pull/380)) v0.19.0 ------ diff --git a/papers/ice_nucleation_2024/perfect_model_J.jl b/papers/ice_nucleation_2024/calibration.jl similarity index 53% rename from papers/ice_nucleation_2024/perfect_model_J.jl rename to papers/ice_nucleation_2024/calibration.jl index ff1f1b838..7aa249166 100644 --- a/papers/ice_nucleation_2024/perfect_model_J.jl +++ b/papers/ice_nucleation_2024/calibration.jl @@ -10,8 +10,10 @@ import CloudMicrophysics as CM import CloudMicrophysics.Parameters as CMP import Thermodynamics as TD +#! format: off # definition of the ODE problem for parcel model include(joinpath(pkgdir(CM), "parcel", "Parcel.jl")) +include(joinpath(pkgdir(CM), "papers", "ice_nucleation_2024", "calibration_setup.jl")) # Define model which wraps around parcel and overwrites calibrated parameters function run_model(p, coefficients, IN_mode, FT, IC) @@ -106,20 +108,34 @@ function run_model(p, coefficients, IN_mode, FT, IC) end # Creating noisy pseudo-observations -function create_prior(FT, IN_mode) +function create_prior(FT, IN_mode, ; perfect_model = false) + # TODO - add perfect_model flag to plot_ensemble_mean.jl observation_data_names = ["m_coeff", "c_coeff"] - # Define prior distributions for the coefficients # stats = [mean, std dev, lower bound, upper bound] - if IN_mode == "ABDINM" - m_stats = [FT(20), FT(1), FT(0), Inf] - c_stats = [FT(-1), FT(1), -Inf, Inf] - elseif IN_mode == "ABIFM" - m_stats = [FT(50), FT(1), FT(0), Inf] - c_stats = [FT(-7), FT(1), -Inf, Inf] - elseif IN_mode == "ABHOM" - m_stats = [FT(260), FT(1), FT(0), Inf] - c_stats = [FT(-70), FT(1), -Inf, Inf] + # 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] + elseif IN_mode == "ABIFM" + m_stats = [FT(50), FT(1), FT(0), Inf] + c_stats = [FT(-7), FT(1), -Inf, Inf] + elseif IN_mode == "ABHOM" + m_stats = [FT(260), FT(1), FT(0), Inf] + c_stats = [FT(-70), FT(1), -Inf, Inf] + end + elseif perfect_model == false + if IN_mode == "ABDINM" # TODO - dependent on dust type + m_stats = [FT(20), FT(1), FT(0), Inf] + c_stats = [FT(-1), FT(1), -Inf, Inf] + elseif IN_mode == "ABIFM" # TODO - dependent on dust type + m_stats = [FT(50), FT(1), FT(0), Inf] + c_stats = [FT(-7), FT(1), -Inf, Inf] + elseif IN_mode == "ABHOM" + m_stats = [FT(255.927125), FT(1), FT(0), Inf] + c_stats = [FT(-68.553283), FT(1), -Inf, Inf] + end end m_prior = EKP.constrained_gaussian( @@ -140,151 +156,16 @@ function create_prior(FT, IN_mode) return prior end -function calibrate_J_parameters(FT, IN_mode) +function calibrate_J_parameters(FT, IN_mode, params, IC, y_truth, Γ,; perfect_model = false) # Random number generator rng_seed = 24 rng = Random.seed!(Random.GLOBAL_RNG, rng_seed) - # Get free parameters - tps = TD.Parameters.ThermodynamicsParameters(FT) - wps = CMP.WaterProperties(FT) - - # Define other parameters and initial conditions - ρₗ = wps.ρw - R_d = TD.Parameters.R_d(tps) - R_v = TD.Parameters.R_v(tps) - ϵₘ = R_d / R_v - - const_dt = FT(1) - - # Create prior - prior = create_prior(FT, IN_mode) - - if IN_mode == "ABDINM" - Nₐ = FT(8e6) # according to fig 2 panel b deposition nucleation experiment - 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) - qₗ = FT(0) - qᵢ = FT(0) - - q = TD.PhasePartition.(qᵥ + qₗ + qᵢ, qₗ, qᵢ) - Rₐ = TD.gas_constant_air(tps, q) - eₛ = TD.saturation_vapor_pressure(tps, T₀, TD.Liquid()) - e = eᵥ(qᵥ, p₀, Rₐ, R_v) - Sₗ = FT(e / eₛ) - w = FT(0.4) - IC = [Sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, FT(0)] - - dep_nucleation = "ABDINM" - deposition_growth = "Deposition" - size_distribution = "Monodisperse" - - params = (; - const_dt, - w, - dep_nucleation, - deposition_growth, - size_distribution, - ) - - coeff_true = [FT(27.551), FT(-2.2209)] - - elseif IN_mode == "ABIFM" - Nₐ = FT(0) - Nₗ = FT(8e6) - Nᵢ = FT(0) - r₀ = FT(0.5e-6) # droplets maybe bigger? - p₀ = FT(987.018 * 1e2) - T₀ = FT(212.978) - - 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(0) - - q = TD.PhasePartition.(qᵥ + qₗ + qᵢ, qₗ, qᵢ) - Rₐ = TD.gas_constant_air(tps, q) - eₛ = TD.saturation_vapor_pressure(tps, T₀, TD.Liquid()) - e = eᵥ(qᵥ, p₀, Rₐ, R_v) - Sₗ = FT(e / eₛ) - w = FT(0.4) - IC = [Sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, FT(0)] - Sᵢ = ξ(tps, T₀) * Sₗ - - heterogeneous = "ABIFM" - condensation_growth = "Condensation" - deposition_growth = "Deposition" - size_distribution = "Monodisperse" - - params = (; - const_dt, - w, - heterogeneous, - condensation_growth, - deposition_growth, - size_distribution, - ) - - coeff_true = [FT(54.58834), FT(-10.54758)] - - elseif IN_mode == "ABHOM" - Nₐ = FT(0) - # Nₗ = FT(360 * 1e6) - Nₗ = FT(300 * 1e6) - Nᵢ = FT(0) - #r₀ = FT(2e-7) # droplets maybe bigger? - r₀ = FT(25 * 1e-9) - #p₀ = FT(987.345 * 1e2) - p₀ = FT(9712.183) - # T₀ = FT(243.134) - 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(357.096 * 1e-6 - C_l) # concentration/mol fraction of vapor - C_v = FT(5 * 1e-6) - qᵥ = ϵₘ / (ϵₘ - 1 + 1 / C_v) - qᵢ = FT(0) - - q = TD.PhasePartition.(qᵥ + qₗ + qᵢ, qₗ, qᵢ) - Rₐ = TD.gas_constant_air(tps, q) - eₛ = TD.saturation_vapor_pressure(tps, T₀, TD.Liquid()) - e = eᵥ(qᵥ, p₀, Rₐ, R_v) - Sₗ = FT(e / eₛ) - w = FT(1) - IC = [Sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, FT(0)] - Sᵢ = ξ(tps, T₀) * Sₗ - - homogeneous = "ABHOM" - deposition_growth = "Deposition" - size_distribution = "Monodisperse" - - params = - (; const_dt, w, homogeneous, deposition_growth, size_distribution) - - coeff_true = [FT(255.927125), FT(-68.553283)] - end - - n_samples = 10 - G_truth = run_model(params, coeff_true, IN_mode, FT, IC) # ICNC from running parcel w default values - y_truth = zeros(length(G_truth), n_samples) # where noisy ICNC will be stored - - dim_output = length(G_truth) - Γ = 0.005 * LinearAlgebra.I * (maximum(G_truth) - minimum(G_truth)) - noise_dist = Distributions.MvNormal(zeros(dim_output), Γ) - y_truth = G_truth .+ rand(noise_dist) - - # Generate initial ensemble and set up EKI - # TODO - make compatible with Float32. Only works with Float64 for now. + prior = create_prior(FT, IN_mode, perfect_model = perfect_model) N_ensemble = 10 # runs N_ensemble trials per iteration N_iterations = 15 # 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) EKI_obj = EKP.EnsembleKalmanProcess( initial_ensemble, @@ -320,7 +201,7 @@ function calibrate_J_parameters(FT, IN_mode) digits = 6, ) - return [m_coeff_ekp, c_coeff_ekp, coeff_true[1], coeff_true[2], ϕ_n_values] + return [m_coeff_ekp, c_coeff_ekp, ϕ_n_values] end function ensemble_means(ϕ_n_values, N_iterations, N_ensemble) diff --git a/papers/ice_nucleation_2024/calibration_setup.jl b/papers/ice_nucleation_2024/calibration_setup.jl new file mode 100644 index 000000000..7e7aa6d22 --- /dev/null +++ b/papers/ice_nucleation_2024/calibration_setup.jl @@ -0,0 +1,128 @@ +import ClimaParams as CP +import CloudMicrophysics as CM +import CloudMicrophysics.Parameters as CMP +import Thermodynamics as TD + +#! format: off +function perf_model_params(FT, IN_mode) + if IN_mode == "ABDINM" + const_dt = FT(1) + w = FT(0.4) + dep_nucleation = "ABDINM" + heterogeneous = "None" + homogeneous = "None" + condensation_growth = "None" + deposition_growth = "Deposition" + size_distribution = "Monodisperse" + elseif IN_mode == "ABIFM" + const_dt = FT(1) + w = FT(0.4) + dep_nucleation = "None" + heterogeneous = "ABIFM" + homogeneous = "None" + condensation_growth = "Condensation" + deposition_growth = "Deposition" + size_distribution = "Monodisperse" + elseif IN_mode == "ABHOM" + const_dt = FT(1) + w = FT(1) + dep_nucleation = "None" + heterogeneous = "None" + homogeneous = "ABHOM" + condensation_growth = "None" + deposition_growth = "Deposition" + size_distribution = "Monodisperse" + end + params = (; const_dt, w, + condensation_growth, deposition_growth, # growth + size_distribution, # size distribution + dep_nucleation, heterogeneous, homogeneous, # nucleation + ) + return params +end + +function perf_model_IC(FT, IN_mode) + tps = TD.Parameters.ThermodynamicsParameters(FT) + wps = CMP.WaterProperties(FT) + + ρₗ = wps.ρw + R_d = TD.Parameters.R_d(tps) + R_v = TD.Parameters.R_v(tps) + ϵₘ = R_d / R_v + + if IN_mode == "ABDINM" + Nₐ = FT(8e6) # according to fig 2 panel b deposition nucleation experiment + 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) + qₗ = FT(0) + qᵢ = FT(0) + q = TD.PhasePartition.(qᵥ + qₗ + qᵢ, qₗ, qᵢ) + Rₐ = TD.gas_constant_air(tps, q) + eₛ = TD.saturation_vapor_pressure(tps, T₀, TD.Liquid()) + e = eᵥ(qᵥ, p₀, Rₐ, R_v) + Sₗ = FT(e / eₛ) + elseif IN_mode == "ABIFM" + Nₐ = FT(0) + Nₗ = FT(8e6) + Nᵢ = FT(0) + r₀ = FT(0.5e-6) + p₀ = FT(987.018 * 1e2) + T₀ = FT(212.978) + 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(0) + q = TD.PhasePartition.(qᵥ + qₗ + qᵢ, qₗ, qᵢ) + Rₐ = TD.gas_constant_air(tps, q) + eₛ = TD.saturation_vapor_pressure(tps, T₀, TD.Liquid()) + e = eᵥ(qᵥ, p₀, Rₐ, R_v) + Sₗ = FT(e / eₛ) + elseif IN_mode == "ABHOM" + Nₐ = FT(0) + Nₗ = FT(300 * 1e6) + Nᵢ = FT(0) + r₀ = FT(25 * 1e-9) + 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) + q = TD.PhasePartition.(qᵥ + qₗ + qᵢ, qₗ, qᵢ) + Rₐ = TD.gas_constant_air(tps, q) + eₛ = TD.saturation_vapor_pressure(tps, T₀, TD.Liquid()) + e = eᵥ(qᵥ, p₀, Rₐ, R_v) + Sₗ = FT(e / eₛ) + end + return [Sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, FT(0)] +end + +function perf_model_pseudo_data(FT, IN_mode, params, IC) + n_samples = 10 + + if IN_mode == "ABDINM" + coeff_true = [FT(27.551), FT(-2.2209)] + elseif IN_mode == "ABIFM" + coeff_true = [FT(54.58834), FT(-10.54758)] + elseif IN_mode == "ABHOM" + coeff_true = [FT(255.927125), FT(-68.553283)] + end + + 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)) + 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) + return [y_truth, Γ, coeff_true] +end +#! format: on diff --git a/papers/ice_nucleation_2024/plot_ensemble_mean.jl b/papers/ice_nucleation_2024/plot_ensemble_mean.jl index 643c9c56a..d11e01789 100644 --- a/papers/ice_nucleation_2024/plot_ensemble_mean.jl +++ b/papers/ice_nucleation_2024/plot_ensemble_mean.jl @@ -3,76 +3,68 @@ import CloudMicrophysics as CM # Float32 does not work because inconsistent types using EKP FT = Float64 -include( - joinpath(pkgdir(CM), "papers", "ice_nucleation_2024", "perfect_model_J.jl"), -) +include(joinpath(pkgdir(CM), "papers", "ice_nucleation_2024", "calibration.jl")) -ABDINM_output = calibrate_J_parameters(FT, "ABDINM") -ABIFM_output = calibrate_J_parameters(FT, "ABIFM") -ABHOM_output = calibrate_J_parameters(FT, "ABHOM") -iterations = collect(1:size(ABHOM_output[5])[1]) +IN_mode_list = ["ABDINM", "ABIFM", "ABHOM"] -ABDINM_calibrated_parameters = [ABDINM_output[1], ABDINM_output[2]] -ABDINM_coeff_true = [ABDINM_output[3], ABDINM_output[4]] -ABIFM_calibrated_parameters = [ABIFM_output[1], ABIFM_output[2]] -ABIFM_coeff_true = [ABIFM_output[3], ABIFM_output[4]] -ABHOM_calibrated_parameters = [ABHOM_output[1], ABHOM_output[2]] -ABHOM_coeff_true = [ABHOM_output[3], ABHOM_output[4]] - -# ensemble_means(ϕ_n_values, N_iterations, ensembles) -ABDINM_ensemble_means = ensemble_means( - ABDINM_output[5], - size(ABDINM_output[5])[1], - size(ABDINM_output[5][1])[2], -) -ABIFM_ensemble_means = ensemble_means( - ABIFM_output[5], - size(ABIFM_output[5])[1], - size(ABIFM_output[5][1])[2], -) -ABHOM_ensemble_means = ensemble_means( - ABHOM_output[5], - size(ABHOM_output[5])[1], - size(ABHOM_output[5][1])[2], -) - -# Plotting calibrated parameters per iteration #! format: off -ABDINM_fig = MK.Figure(size = (800, 600)) -ax1 = MK.Axis(ABDINM_fig[1, 1], ylabel = "m coefficient [-]", xlabel = "iteration number", title = "ABDINM") -ax2 = MK.Axis(ABDINM_fig[2, 1], ylabel = "c coefficient [-]", xlabel = "iteration number") +for IN_mode in IN_mode_list + params = perf_model_params(FT, IN_mode) + IC = perf_model_IC(FT, IN_mode) -MK.lines!(ax1, iterations, ABDINM_ensemble_means[1], label = "ensemble mean", color = :orange) -MK.lines!(ax1, iterations, zeros(length(ABDINM_ensemble_means[1])) .+ ABDINM_output[3], label = "default value") -MK.lines!(ax2, iterations, ABDINM_ensemble_means[2], label = "ensemble mean", color = :orange) -MK.lines!(ax2, iterations, zeros(length(ABDINM_ensemble_means[1])) .+ ABDINM_output[4], label = "default value") + pseudo_data = perf_model_pseudo_data(FT, IN_mode, params, IC) + Γ = pseudo_data[2] + y_truth = pseudo_data[1] + coeff_true = pseudo_data[3] -ABIFM_fig = MK.Figure(size = (800, 600)) -ax3 = MK.Axis(ABIFM_fig[1, 1], ylabel = "m coefficient [-]", xlabel = "iteration number", title = "ABIFM") -ax4 = MK.Axis(ABIFM_fig[2, 1], ylabel = "c coefficient [-]", xlabel = "iteration number") + output = calibrate_J_parameters( + FT, + IN_mode, + params, + IC, + y_truth, + Γ, + perfect_model = true, + ) -MK.lines!(ax3, iterations, ABIFM_ensemble_means[1], label = "ensemble mean", color = :orange) -MK.lines!(ax3, iterations, zeros(length(ABIFM_ensemble_means[1])) .+ ABIFM_output[3], label = "default value") -MK.lines!(ax4, iterations, ABIFM_ensemble_means[2], label = "ensemble mean", color = :orange) -MK.lines!(ax4, iterations, zeros(length(ABIFM_ensemble_means[1])) .+ ABIFM_output[4], label = "default value") + iterations = collect(1:size(output[3])[1]) + calibrated_parameters = [output[1], output[2]] + m_mean = [] + m_mean = ensemble_means( + output[3], + size(output[3])[1], + size(output[3][1])[2], + )[1] + c_mean = [] + c_mean = ensemble_means( + output[3], + size(output[3])[1], + size(output[3][1])[2], + )[2] -ABHOM_fig = MK.Figure(size = (800, 600)) -ax5 = MK.Axis(ABHOM_fig[1, 1], ylabel = "m coefficient [-]", xlabel = "iteration number", title = "ABHOM") -ax6 = MK.Axis(ABHOM_fig[2, 1], ylabel = "c coefficient [-]", xlabel = "iteration number") + # Plotting calibrated parameters per iteration + fig = MK.Figure(size = (800, 600)) + ax1 = MK.Axis(fig[1, 1], ylabel = "m coefficient [-]", xlabel = "iteration number", title = IN_mode) + ax2 = MK.Axis(fig[2, 1], ylabel = "c coefficient [-]", xlabel = "iteration number") + + MK.lines!(ax1, iterations, m_mean, label = "ensemble mean", color = :orange) + MK.lines!(ax1, iterations, zeros(length(m_mean)) .+ coeff_true[1], label = "default value") + MK.lines!(ax2, iterations, c_mean, label = "ensemble mean", color = :orange) + MK.lines!(ax2, iterations, zeros(length(c_mean)) .+ coeff_true[2], label = "default value") + + MK.axislegend(ax1, framevisible = true, labelsize = 12, position = :rc) + MK.axislegend(ax2, framevisible = true, labelsize = 12) -MK.lines!(ax5, iterations, ABHOM_ensemble_means[1], label = "ensemble mean", color = :orange) -MK.lines!(ax5, iterations, zeros(length(ABHOM_ensemble_means[1])) .+ ABHOM_output[3], label = "default value") -MK.lines!(ax6, iterations, ABHOM_ensemble_means[2], label = "ensemble mean", color = :orange) -MK.lines!(ax6, iterations, zeros(length(ABHOM_ensemble_means[1])) .+ ABHOM_output[4], label = "default value") -#! format: on + if IN_mode == "ABDINM" + mode_label = "dep" + elseif IN_mode == "ABIFM" + mode_label = "imm" + elseif IN_mode == "ABHOM" + mode_label = "hom" + end -MK.axislegend(ax1, framevisible = true, labelsize = 12, position = :rc) -MK.axislegend(ax2, framevisible = true, labelsize = 12) -MK.axislegend(ax3, framevisible = true, labelsize = 12, position = :rc) -MK.axislegend(ax4, framevisible = true, labelsize = 12) -MK.axislegend(ax5, framevisible = true, labelsize = 12, position = :rc) -MK.axislegend(ax6, framevisible = true, labelsize = 12) + plot_name = "perfect_calibration_$mode_label.svg" + MK.save(plot_name, fig) -MK.save("perfect_calibration_dep.svg", ABDINM_fig) -MK.save("perfect_calibration_imm.svg", ABIFM_fig) -MK.save("perfect_calibration_hom.svg", ABHOM_fig) +end +#! format: on diff --git a/test/ice_nucleation_calibration.jl b/test/ice_nucleation_calibration.jl index dd06dda83..4dfbe8210 100644 --- a/test/ice_nucleation_calibration.jl +++ b/test/ice_nucleation_calibration.jl @@ -3,20 +3,33 @@ import Test as TT # definition of the ODE problem for parcel model include(joinpath(pkgdir(CM), "parcel", "Parcel.jl")) -include( - joinpath(pkgdir(CM), "papers", "ice_nucleation_2024", "perfect_model_J.jl"), -) +include(joinpath(pkgdir(CM), "papers", "ice_nucleation_2024", "calibration.jl")) function test_J_calibration(FT, IN_mode) - output = calibrate_J_parameters(FT, IN_mode) + params = perf_model_params(FT, IN_mode) + IC = perf_model_IC(FT, IN_mode) + + pseudo_data = perf_model_pseudo_data(FT, IN_mode, params, IC) + Γ = pseudo_data[2] + y_truth = pseudo_data[1] + coeff_true = pseudo_data[3] + + output = calibrate_J_parameters( + FT, + IN_mode, + params, + IC, + y_truth, + Γ, + perfect_model = true, + ) calibrated_parameters = [output[1], output[2]] - coeff_true = [output[3], output[4]] TT.@test calibrated_parameters[1] ≈ coeff_true[1] rtol = FT(0.15) - TT.@test calibrated_parameters[2] ≈ coeff_true[2] rtol = FT(0.15) + TT.@test calibrated_parameters[2] ≈ coeff_true[2] rtol = FT(0.17) end -@info "Ice Nucleation Test" +@info "Ice Nucleation Calibration Test" TT.@testset "Perfect model calibration on ABDINM" begin test_J_calibration(Float64, "ABDINM") diff --git a/test/performance_tests.jl b/test/performance_tests.jl index c1c22400d..b50b18436 100644 --- a/test/performance_tests.jl +++ b/test/performance_tests.jl @@ -187,7 +187,7 @@ function benchmark_test(FT) bench_press(CMN.τ_relax, (liquid,), 10) # 0-moment - bench_press(CM0.remove_precipitation, (p0m, q), 10) + bench_press(CM0.remove_precipitation, (p0m, q), 12) # 1-moment bench_press(