Skip to content

Commit

Permalink
Merge branch 'al/HOM_calibs' into al/aida_calibrations
Browse files Browse the repository at this point in the history
  • Loading branch information
@amylu00
amylu00 committed Sep 6, 2024
2 parents e39b147 + 93d7329 commit 0e1bace
Show file tree
Hide file tree
Showing 2 changed files with 60 additions and 53 deletions.
71 changes: 39 additions & 32 deletions papers/ice_nucleation_2024/AIDA_calibrations.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,39 +20,45 @@ end
moving_average(data, n) =
[sum(@view data[i:(i + n - 1)]) / n for i in 1:(length(data) - (n - 1))]

data_file_names = ["in05_17_aida.edf"] #, "in05_18_aida.edf", "in05_19_aida.edf"]
plot_names = ["IN0517"] #, "IN0518", "IN0519"]

start_time = 100 + 195 # AIDA time starts at -100 seconds. Add freezing onset as needed.
end_time = 100 + 290 # might change per experiment. time freezing stops.
moving_average_n = 20 # average every n points
moving_average_start = Int32(start_time + (moving_average_n / 2 - 1))
moving_average_end = Int32(end_time - (moving_average_n / 2))
t_max =
(end_time - 100 - (moving_average_n / 2)) - # AIDA time starts at -100 seconds (t = 0 at index 101)
(start_time - 100 + (moving_average_n / 2 - 1)) # duration of simulation

updrafts = [FT(1.5), FT(1.2), FT(3.7)]
data_file_names = ["in05_17_aida.edf", "in05_18_aida.edf", "in05_19_aida.edf"]
plot_names = ["IN0517", "IN0518", "IN0519"]

start_time_list = [195, 180, 80] .+ 100 # AIDA time starts at -100 seconds. Add freezing onset as needed.
end_time_list = [290, 290, 170] .+ 100 # approximate time freezing stops
moving_average_n = 20 # average every n points

updrafts = [FT(1.5), FT(1.4), FT(5)]
tps = TD.Parameters.ThermodynamicsParameters(FT)
R_v = TD.Parameters.R_v(tps)
R_d = TD.Parameters.R_d(tps)
ϵₘ = R_d / R_v


for (index, data_file_name) in enumerate(data_file_names)
for (exp_index, data_file_name) in enumerate(data_file_names)
#! format: off
### Unpacking experiment-specific variables
plot_name = plot_names[index]
w = updrafts[index]
plot_name = plot_names[exp_index]
w = updrafts[exp_index]
start_time = start_time_list[exp_index]
end_time = end_time_list[exp_index]

## Moving average to smooth data
moving_average_start = Int32(start_time + (moving_average_n / 2 - 1))
moving_average_end = Int32(end_time - (moving_average_n / 2))
t_max =
(end_time - 100 - (moving_average_n / 2)) - # AIDA time starts at -100 seconds (t = 0 at index 101)
(start_time - 100 + (moving_average_n / 2 - 1)) # duration of simulation

params = AIDA_IN05_params(FT, w, t_max)
IC = AIDA_IN05_IC(FT, data_file_name)

### 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]

frozen_frac = ICNC ./ (IC[7] + IC[8] + IC[9])
frozen_frac = ICNC ./ (IC[7] + IC[8] + IC[9])
frozen_frac_moving_mean = moving_average(frozen_frac, moving_average_n)
ICNC_moving_avg = moving_average(ICNC, moving_average_n)
Γ = 0.15 * LinearAlgebra.I * (maximum(frozen_frac_moving_mean) - minimum(frozen_frac_moving_mean)) # coeff is an estimated of the noise

### Plotting AIDA data.
Expand All @@ -67,13 +73,14 @@ for (index, data_file_name) in enumerate(data_file_names)
linestyle =:dash,
linewidth = 2,
)
# MK.lines!(
# ax1,
# chamber_data[moving_average_start:moving_average_end, 1],
# frozen_frac_moving_mean, label = "AIDA moving mean",
# linewidth = 2.5,
# color =:blue,
# )
MK.lines!(
ax1,
chamber_data[moving_average_start:moving_average_end, 1],
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)

Expand All @@ -87,7 +94,7 @@ for (index, data_file_name) in enumerate(data_file_names)
parcel = run_calibrated_model(FT, "ABHOM", calibrated_parameters, params, IC)
#parcel_default = run_calibrated_model(FT, "ABHOM", [FT(255.927125), FT(-68.553283)], params, IC)

### Plots
### Plots.
## Calibrated coefficients.
# Did they converge?
calibrated_coeffs_fig = MK.Figure(size = (600, 600), fontsize = 24)
Expand Down Expand Up @@ -124,13 +131,13 @@ for (index, data_file_name) in enumerate(data_file_names)
MK.lines!(ax_parcel_1, parcel.t .+ (moving_average_n / 2), parcel[1, :], label = "calibrated", color = :orange) # label = "liquid"
#MK.lines!(ax_parcel_1, parcel_default.t .+ (moving_average_n / 2), parcel_default[1, :], label = "default", color = :darkorange2)
MK.lines!(ax_parcel_1, parcel.t .+ (moving_average_n / 2), S_i.(tps, parcel[3, :], parcel[1, :]), label = "ice", color = :green)
MK.lines!(
ax_parcel_1,
chamber_data[start_time:end_time, 1] .- (start_time - 100),
vec(chamber_S_l),
label = "chamber",
color = :blue
)
# MK.lines!(
# ax_parcel_1,
# chamber_data[start_time:end_time, 1] .- (start_time - 100),
# vec(chamber_S_l),
# label = "chamber",
# color = :blue
# )
MK.lines!(ax_parcel_2, parcel.t .+ (moving_average_n / 2), parcel[5, :], color = :orange)
#MK.lines!(ax_parcel_2, parcel_default.t .+ (moving_average_n / 2), parcel_default[5,:], color = :darkorange2)
MK.lines!(ax_parcel_3, parcel.t .+ (moving_average_n / 2), parcel[3, :], color = :orange)
Expand Down
42 changes: 21 additions & 21 deletions papers/ice_nucleation_2024/calibration_setup.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -181,9 +181,9 @@ function AIDA_IN05_IC(FT, data_file)

if data_file == "in05_17_aida.edf"
# starting at t = 205 s (to match moving average freezing onset)
Nₐ = FT(360 * 1e6 - 297.136 * 1e6)
Nₗ = FT(297.136 * 1e6)
Nᵢ = FT(1.49 * 1e6)
Nₐ = FT(360 * 1e6) - Nₗ - Nᵢ
r₀ = FT(6.17664e-6) # FT(1e-7)
p₀ = FT(865.179 * 1e2)
T₀ = FT(236.91)
Expand All @@ -198,38 +198,38 @@ function AIDA_IN05_IC(FT, data_file)
eₛ = TD.saturation_vapor_pressure(tps, T₀, TD.Liquid())
Sₗ = FT(e / eₛ)
elseif data_file == "in05_18_aida.edf"
Nₐ = FT(0)
Nₗ = FT(275 * 1e6)
Nᵢ = FT(0)
r₀ = FT(1e-7)
p₀ = FT(996.472 * 1e2)
T₀ = FT(243.072)
Nₗ = FT(209.46 * 1e6)
Nᵢ = FT(1.53 * 1e6)
Nₐ = FT(275 * 1e6) - Nₗ - Nᵢ
r₀ = FT(7.03467 * 1e-6)
p₀ = FT(873.212 * 1e2)
T₀ = FT(237.134)
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(355.898 * 1e-6 - C_l) # concentration/mol fraction of vapor
qᵥ = ϵₘ / (ϵₘ - 1 + 1 / C_v)
qᵢ = FT(0)
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(28.1324)
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())
e = eᵥ(qᵥ, p₀, Rₐ, R_v)
Sₗ = FT(e / eₛ)
elseif data_file == "in05_19_aida.edf"
Nₐ = FT(0)
Nₗ = FT(180 * 1e6)
Nᵢ = FT(0)
r₀ = FT(1e-7)
p₀ = FT(795.081 * 1e2)
T₀ = FT(242.807)
Nᵢ = FT(0.49 * 1e6)
r₀ = FT(6.5e-6) # !!missing in dataset!!
p₀ = FT(722.852 * 1e2)
T₀ = FT(237.521)
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(432.364 * 1e-6 - C_l)
qᵥ = ϵₘ / (ϵₘ - 1 + 1 / C_v)
qᵢ = FT(0)
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(29.5341) # !!missing in dataset!!
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())
e = eᵥ(qᵥ, p₀, Rₐ, R_v)
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
Expand Down

0 comments on commit 0e1bace

Please sign in to comment.