Skip to content

Commit

Permalink
Merge pull request #380 from CliMA/al/aida_calib
Browse files Browse the repository at this point in the history
Generalize ice nucleation calibration funcs
  • Loading branch information
amylu00 authored Apr 26, 2024
2 parents a11f7db + c72834c commit f726c33
Show file tree
Hide file tree
Showing 6 changed files with 237 additions and 222 deletions.
1 change: 1 addition & 0 deletions NEWS.md
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,7 @@ CloudMicrophysics.jl Release Notes
main
------
<!--- # Add changes since the most recent release here --->
- Generalize calibration functions in ice_nucleation_2024 ([#380](https://github.com/CliMA/CloudMicrophysics.jl/pull/380))

v0.19.0
------
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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(
Expand All @@ -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,
Expand Down Expand Up @@ -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)
Expand Down
128 changes: 128 additions & 0 deletions papers/ice_nucleation_2024/calibration_setup.jl
Original file line number Diff line number Diff line change
@@ -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
Loading

0 comments on commit f726c33

Please sign in to comment.