diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml index 883489975..bb423c9f5 100644 --- a/.buildkite/pipeline.yml +++ b/.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 diff --git a/papers/ice_nucleation_2024/AIDA_calibrations.jl b/papers/ice_nucleation_2024/AIDA_calibrations.jl index 16700a0e9..2d294fa46 100644 --- a/papers/ice_nucleation_2024/AIDA_calibrations.jl +++ b/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") diff --git a/papers/ice_nucleation_2024/AIDA_data_artifact/README.md b/papers/ice_nucleation_2024/AIDA_data_artifact/README.md index 0066266af..a3833cd2f 100644 --- a/papers/ice_nucleation_2024/AIDA_data_artifact/README.md +++ b/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) diff --git a/papers/ice_nucleation_2024/calibration.jl b/papers/ice_nucleation_2024/calibration.jl index e7af60f1f..2c2dfcb95 100644 --- a/papers/ice_nucleation_2024/calibration.jl +++ b/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 diff --git a/papers/ice_nucleation_2024/calibration_setup.jl b/papers/ice_nucleation_2024/calibration_setup.jl index 83b2f3683..14d753289 100644 --- a/papers/ice_nucleation_2024/calibration_setup.jl +++ b/papers/ice_nucleation_2024/calibration_setup.jl @@ -195,7 +195,6 @@ function AIDA_IN05_IC(FT, data_file) ϵₘ = R_d / R_v if data_file == "in05_17_aida.edf" - # starting at t = 205 s (to match moving average freezing onset) Nₗ = FT(297.136 * 1e6) Nᵢ = FT(1.49 * 1e6) Nₐ = FT(360 * 1e6) - Nₗ - Nᵢ @@ -247,8 +246,84 @@ function AIDA_IN05_IC(FT, data_file) eₛ = TD.saturation_vapor_pressure(tps, T₀, TD.Liquid()) Sₗ = FT(e / eₛ) end - return [Sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, FT(0)] #remove the last 2 elements, its J & r_l + return [Sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, FT(0)] end +function AIDA_IN07_params(FT, w, t_max, t_profile, T_profile, P_profile, plot_name) + IN_mode = "ABDINM" + const_dt = FT(1) + prescribed_thermodynamics = true + aerosol_act = "None" + aerosol = plot_name == "IN0701" ? CMP.ArizonaTestDust(FT) : CMP.Illite(FT) + dep_nucleation = "ABDINM" + heterogeneous = "None" + homogeneous = "None" + condensation_growth = "Condensation" + deposition_growth = "Deposition" + liq_size_distribution = "Gamma" + ice_size_distribution = "Gamma" + aero_σ_g = FT(2.3) + r_nuc = FT(0.1742857 * 1e-6) + # r_nuc = r₀ in IC + ips = CMP.IceNucleationParameters(FT) + + params = (; const_dt, w, t_max,ips, + prescribed_thermodynamics, t_profile, T_profile, P_profile, + aerosol_act, aerosol, r_nuc, aero_σ_g, # aerosol activation + condensation_growth, deposition_growth, # growth + liq_size_distribution, ice_size_distribution, # size distribution + dep_nucleation, heterogeneous, homogeneous, # ice nucleation + ) + return params +end + +function AIDA_IN07_IC(FT, data_file) + tps = TD.Parameters.ThermodynamicsParameters(FT) + wps = CMP.WaterProperties(FT) + + ρₗ = wps.ρw + ρᵢ = wps.ρi + R_d = TD.Parameters.R_d(tps) + R_v = TD.Parameters.R_v(tps) + ϵₘ = R_d / R_v + if data_file == "in07_01_aida.edf" + Nₗ = FT(0) # TODO - idk what this should be + Nᵢ = FT(1.485 * 1e6) + Nₐ = FT(146 * 1e6) - Nₗ - Nᵢ + r₀ = FT(0.19801 * 1e-6) + # r₀ is weighted avg of two modes as listed in Table 2 of Mohler et al (2008) + p₀ = FT(976.752 * 1e2) + T₀ = FT(208.671) + qₗ = FT(Nₗ * 4 / 3 * FT(π) * r₀^3 * ρₗ / FT(1.2)) # 1.2 should be ρₐ + 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.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) + eₛ = TD.saturation_vapor_pressure(tps, T₀, TD.Liquid()) + Sₗ = FT(e / eₛ) + elseif data_file == "in07_19_aida.edf" + Nₗ = FT(0) # TODO - idk what this should be + Nᵢ = FT(0.21 * 1e6) + Nₐ = FT(95 * 1e6) - Nₗ - Nᵢ + r₀ = FT(0.1742857 * 1e-6) + # r₀ is weighted avg of two modes as listed in Table 2 of Mohler et al (2008) + p₀ = FT(977.948 * 1e2) + T₀ = FT(208.672) + qₗ = FT(Nₗ * 4 / 3 * FT(π) * r₀^3 * ρₗ / FT(1.2)) # 1.2 should be ρₐ + 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.624557) + 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) + eₛ = TD.saturation_vapor_pressure(tps, T₀, TD.Liquid()) + Sₗ = FT(e / eₛ) + end + return [Sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, FT(0)] +end #! format: on diff --git a/parcel/ParcelModel.jl b/parcel/ParcelModel.jl index 070ca2466..88e6d663e 100644 --- a/parcel/ParcelModel.jl +++ b/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 diff --git a/src/IceNucleation.jl b/src/IceNucleation.jl index dff2d5023..801cb29b5 100644 --- a/src/IceNucleation.jl +++ b/src/IceNucleation.jl @@ -83,7 +83,13 @@ see DOI: 10.1002/2016JD025817 Returns zero for unsupported aerosol types. """ function deposition_J( - dust::Union{CMP.Ferrihydrite, CMP.Feldspar, CMP.Kaolinite}, + dust::Union{ + CMP.Ferrihydrite, + CMP.Feldspar, + CMP.Kaolinite, + CMP.Illite, + CMP.ArizonaTestDust + }, Δa_w::FT, ) where {FT} @@ -92,6 +98,7 @@ function deposition_J( return max(FT(0), FT(10)^logJ * FT(1e4)) # converts cm^-2 s^-1 to m^-2 s^-1 end function deposition_J(dust::CMP.AerosolType, Δa_w::FT) where {FT} + println("Aerosol type not supported for ABDINM.") return FT(0) end diff --git a/src/parameters/AerosolATD.jl b/src/parameters/AerosolATD.jl index 08526e933..f37ce2adc 100644 --- a/src/parameters/AerosolATD.jl +++ b/src/parameters/AerosolATD.jl @@ -18,6 +18,14 @@ Base.@kwdef struct ArizonaTestDust{FT} <: AerosolType{FT} a_warm::FT "a for T < T_thr [-]" a_cold::FT + "m coefficient for deposition nucleation J [-]" + deposition_m::FT + "c coefficient for deposition nucleation J [-]" + deposition_c::FT + "m coefficient for immersion freezing J [-]" + ABIFM_m::FT + "c coefficient for immersion freezing J [-]" + ABIFM_c::FT end ArizonaTestDust(::Type{FT}) where {FT <: AbstractFloat} = @@ -29,6 +37,10 @@ function ArizonaTestDust(td::CP.AbstractTOMLDict) :Mohler2006_S0_cold_ArizonaTestDust => :S₀_cold, :Mohler2006_a_warm_ArizonaTestDust => :a_warm, :Mohler2006_a_cold_ArizonaTestDust => :a_cold, + :J_ABDINM_m_ArizonaTestDust => :deposition_m, + :J_ABDINM_c_ArizonaTestDust => :deposition_c, + :J_ABIFM_m_ArizonaTestDust => :ABIFM_m, + :J_ABIFM_c_ArizonaTestDust => :ABIFM_c, ) parameters = CP.get_parameter_values(td, name_map, "CloudMicrophysics") FT = CP.float_type(td) diff --git a/src/parameters/AerosolIllite.jl b/src/parameters/AerosolIllite.jl index 2730bc1e2..e6a3a343c 100644 --- a/src/parameters/AerosolIllite.jl +++ b/src/parameters/AerosolIllite.jl @@ -10,6 +10,10 @@ DOI: 10.1039/C3FD00035D $(DocStringExtensions.FIELDS) """ Base.@kwdef struct Illite{FT} <: AerosolType{FT} + "m coefficient for deposition nucleation J [-]" + deposition_m::FT + "c coefficient for deposition nucleation J [-]" + deposition_c::FT "m coefficient for immersion freezing J [-]" ABIFM_m::FT "c coefficient for immersion freezing J [-]" @@ -20,6 +24,8 @@ Illite(::Type{FT}) where {FT <: AbstractFloat} = Illite(CP.create_toml_dict(FT)) function Illite(td::CP.AbstractTOMLDict) name_map = (; + :J_ABDINM_m_Illite => :deposition_m, + :J_ABDINM_c_Illite => :deposition_c, :KnopfAlpert2013_J_ABIFM_m_Illite => :ABIFM_m, :KnopfAlpert2013_J_ABIFM_c_Illite => :ABIFM_c, ) diff --git a/test/Project.toml b/test/Project.toml index c5c244f47..11b954321 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -4,6 +4,7 @@ BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" ClimaComms = "3a4d1b5c-c61d-41fd-a00a-5873ba7a1b0d" +ClimaCore = "d414da3d-4745-48bb-8d80-42e94e092884" ClimaParams = "5c42b081-d73a-476f-9059-fd94b934656c" CloudMicrophysics = "6a9e3e04-43cd-43ba-94b9-e8782df3c71b" DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" diff --git a/test/gpu_clima_core_test.jl b/test/gpu_clima_core_test.jl new file mode 100644 index 000000000..fd6475c60 --- /dev/null +++ b/test/gpu_clima_core_test.jl @@ -0,0 +1,175 @@ +#ENV["CLIMACOMMS_DEVICE"] = "CUDA" + +import ClimaComms +ClimaComms.@import_required_backends +import ClimaCore as CC + +import ClimaParams as CP +import CloudMicrophysics.Parameters as CMP +import CloudMicrophysics.MicrophysicsNonEq as CMN + +""" + A helper function to create a ClimaCore 1d column space +""" +function make_column(::Type{FT}) where {FT} + + context = ClimaComms.SingletonCommsContext(ClimaComms.CUDADevice()) + + vert_domain = CC.Domains.IntervalDomain( + CC.Geometry.ZPoint{FT}(FT(0)), + CC.Geometry.ZPoint{FT}(FT(1000)); + boundary_names = (:bottom, :top), + ) + vert_mesh = CC.Meshes.IntervalMesh(vert_domain; nelems = 1000) + vert_topology = CC.Topologies.IntervalTopology(context, vert_mesh) + vert_center_space = CC.Spaces.CenterFiniteDifferenceSpace(vert_topology) + return vert_center_space +end + +""" + A helper function to create a ClimaCore extruded sphere space +""" +function make_extruded_sphere(::Type{FT}) where {FT} + + context = ClimaComms.SingletonCommsContext(ClimaComms.CUDADevice()) + + # Define vertical + # domain + vert_domain = CC.Domains.IntervalDomain( + CC.Geometry.ZPoint{FT}(FT(0)), + CC.Geometry.ZPoint{FT}(FT(1000)); + boundary_names = (:bottom, :top), + ) + # mesh + vert_mesh = CC.Meshes.IntervalMesh(vert_domain; nelems = 1000) + # topology + vert_topology = CC.Topologies.IntervalTopology(context, vert_mesh) + # grid + vert_grid = CC.Grids.FiniteDifferenceGrid(vert_topology) + #vert_center_space = CC.Spaces.CenterFiniteDifferenceSpace(vert_topology) + + # Define horizontal: + # domain + horz_domain = CC.Domains.SphereDomain(FT(30)) + # mesh + horz_mesh = CC.Meshes.EquiangularCubedSphere(horz_domain, 4) + # topology + horz_topology = CC.Topologies.Topology2D( + context, + horz_mesh, + CC.Topologies.spacefillingcurve(horz_mesh), + ) + # space + horz_space = CC.Spaces.SpectralElementSpace2D( + horz_topology, + CC.Quadratures.GLL{3 + 1}(); + enable_bubble = true, + ) + # grid + horz_grid = CC.Spaces.grid(horz_space) + + # Define surface + z_surface = zeros(horz_space) + hypsography = CC.Hypsography.Flat() + + # Define grid + deep = false + grid = CC.Grids.ExtrudedFiniteDifferenceGrid( + horz_grid, + vert_grid, + hypsography; + deep, + ) + + # Define 3D space + center_extruded_space = CC.Spaces.CenterExtrudedFiniteDifferenceSpace(grid) + return center_extruded_space +end + +""" + Try to reproduce the setup of how terminal velocity is used in Atmos +""" +function set_sedimentation_precomputed_quantities(Y, p, t) + + (; w) = p + (; params) = p + + @. w = CMN.terminal_velocity( + params.liquid, + params.Ch2022.rain, + Y.ρ, + max(0, Y.ρq / Y.ρ), + ) + return nothing +end + +function main_1d(::Type{FT}) where {FT} + + Ch2022 = CMP.Chen2022VelType(FT) + liquid = CMP.CloudLiquid(FT) + ice = CMP.CloudIce(FT) + rain = CMP.Rain(FT) + snow = CMP.Snow(FT) + + params = (; liquid, ice, Ch2022) + + space_1d_ρq = make_column(FT) + space_1d_ρ = make_column(FT) + space_1d_w = make_column(FT) + + ρq = CC.Fields.ones(space_1d_ρq) .* FT(1e-3) + ρ = CC.Fields.ones(space_1d_ρ) + w = CC.Fields.zeros(space_1d_w) + + Y = (; ρq, ρ) + p = (; w, params) + + t = 1 + + set_sedimentation_precomputed_quantities(Y, p, t) + return nothing +end + +function main_3d(::Type{FT}) where {FT} + + Ch2022 = CMP.Chen2022VelType(FT) + liquid = CMP.CloudLiquid(FT) + ice = CMP.CloudIce(FT) + rain = CMP.Rain(FT) + snow = CMP.Snow(FT) + + params = (; liquid, ice, Ch2022) + + space_3d_ρq = make_extruded_sphere(FT) + space_3d_ρ = make_extruded_sphere(FT) + space_3d_w = make_extruded_sphere(FT) + + ρq = CC.Fields.ones(space_3d_ρq) .* FT(1e-3) + ρ = CC.Fields.ones(space_3d_ρ) + w = CC.Fields.zeros(space_3d_w) + + Y = (; ρq, ρ) + p = (; w, params) + + t = 1 + + set_sedimentation_precomputed_quantities(Y, p, t) + return nothing +end + +using Test +@testset "GPU inference failure 1D Float64" begin + main_1d(Float64) +end + +@testset "GPU inference failure 3D Float64" begin + main_3d(Float64) +end + +@testset "GPU inference failure 1D Float32" begin + main_1d(Float32) +end + +@testset "GPU inference failure 3D Float32" begin + main_3d(Float32) +end