Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add ML aerosol activation emulators #172

Closed
wants to merge 5 commits into from
Closed
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
9 changes: 9 additions & 0 deletions Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -5,11 +5,20 @@ version = "0.13.0"

[deps]
CLIMAParameters = "6eacf6c3-8458-43b9-ae03-caf5306d3d53"
DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
DataFramesMeta = "1313f7d8-7da2-5740-9ea0-a2ca25f37964"
Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
Flux = "587475ba-b771-5e3f-ad9e-33799f191a9c"
ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
GaussianProcesses = "891a1506-143c-57d2-908e-e1f8e92e6de9"
KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
MLJ = "add582a8-e3ab-11e8-2d5e-e98b27df1bc7"
MLJFlux = "094fc8d1-fd35-5302-93ea-dabda2abf845"
MLJModels = "d491faf4-2d78-11e9-2867-c94bc002c0b7"
SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
Thermodynamics = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c"

[compat]
Expand Down
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
18 changes: 18 additions & 0 deletions docs/src/ARGdata.jl
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,8 @@ Fig1_y_obs = [
0.34946619217081853,
0.17580071174377226,
]
Fig1_x_PySDM = [100.0, 1325.0, 2550.0, 3775.0, 5000.0]
Fig1_y_PySDM = [0.76, 0.5, 0.38, 0.3, 0.26]
Fig1_x_param = [
0,
12.891344383057117,
Expand Down Expand Up @@ -129,6 +131,8 @@ Fig2a_y_obs = [
0.5986159169550174,
0.4975778546712803,
]
Fig2a_x_PySDM = [100.0, 1325.0, 2550.0, 3775.0, 5000.0]
Fig2a_y_PySDM = [0.8, 0.68, 0.62, 0.6, 0.56]
Fig2a_x_param = [
12.567324955116646,
44.88330341113101,
Expand Down Expand Up @@ -191,6 +195,8 @@ Fig2b_y_obs = [
0.1949367088607595,
0.1329113924050633,
]
Fig2b_x_PySDM = [100.0, 1325.0, 2550.0, 3775.0, 5000.0]
Fig2b_y_PySDM = [0.38, 0.26, 0.22, 0.2, 0.18]
Fig2b_x_param = [
78.80910683012257,
113.83537653239921,
Expand Down Expand Up @@ -233,6 +239,8 @@ Fig3a_x_obs = [
1.0034965034965033,
]
Fig3a_y_obs = [0.7586666666666667, 0.7453333333333334, 0.732, 0.72, 0.712]
Fig3a_x_PySDM = [0.1, 0.325, 0.55, 0.775, 1.0]
Fig3a_y_PySDM = [0.8, 0.78, 0.76, 0.76, 0.76]
Fig3a_x_param = [
0.10069930069930072,
0.14405594405594407,
Expand Down Expand Up @@ -269,6 +277,8 @@ Fig3b_y_obs = [
0.6734982332155477,
0.7074204946996467,
]
Fig3b_x_PySDM = [0.1, 0.325, 0.55, 0.775, 1.0]
Fig3b_y_PySDM = [0.38, 0.58, 0.66, 0.72, 0.76]
Fig3b_x_param = [
0.10270270270270272,
0.12000000000000002,
Expand Down Expand Up @@ -321,6 +331,8 @@ Fig4a_y_obs = [
0.4013333333333333,
0.12266666666666667,
]
Fig4a_x_PySDM = [0.01, 0.02659148, 0.07071068, 0.18803015, 0.5]
Fig4a_y_PySDM = [0.82, 0.78, 0.74, 0.58, 0.18]
Fig4a_x_param = [
0.010398558176237407,
0.011922822170487846,
Expand Down Expand Up @@ -379,6 +391,8 @@ Fig4b_y_obs = [
0.9805194805194806,
0.9844155844155846,
]
Fig4b_x_PySDM = [0.01, 0.02659148, 0.07071068, 0.18803015, 0.5]
Fig4b_y_PySDM = [0.08, 0.46, 0.86, 0.98, 1.0]
Fig4b_x_param = [
0.010822678662544178,
0.012760436250146357,
Expand Down Expand Up @@ -441,6 +455,8 @@ Fig5a_y_obs = [
0.9223427331887204,
0.975704989154013,
]
Fig5a_x_PySDM = [0.01, 0.04728708, 0.2236068, 1.05737126, 5.0]
Fig5a_y_PySDM = [0.12, 0.42, 0.68, 0.88, 0.98]
Fig5a_x_param = [
0.012498045610787637,
0.0206913808111479,
Expand Down Expand Up @@ -493,6 +509,8 @@ Fig5b_y_obs = [
0.6230508474576271,
0.7993220338983051,
]
Fig5b_x_PySDM = [0.01, 0.04728708, 0.2236068, 1.05737126, 5.0]
Fig5b_y_PySDM = [0.02, 0.1, 0.26, 0.54, 0.84]
Fig5b_x_param = [
0.011331300304886677,
0.013602279391211175,
Expand Down
159 changes: 138 additions & 21 deletions docs/src/ARGplots.jl
Original file line number Diff line number Diff line change
@@ -1,15 +1,14 @@
import Plots
import Plots as PL

import CloudMicrophysics
import CLIMAParameters
import Thermodynamics

const PL = Plots
const AM = CloudMicrophysics.AerosolModel
const AA = CloudMicrophysics.AerosolActivation
const CP = CLIMAParameters
const CMP = CloudMicrophysics.Parameters
const TD = Thermodynamics
import CloudMicrophysics:
AerosolModel as AM,
AerosolActivation as AA,
Parameters as CMP,
CommonTypes as CMT
import CLIMAParameters as CP
import Thermodynamics as TD
import DataFrames as DF

include(joinpath(pkgdir(CloudMicrophysics), "test", "create_parameters.jl"))
FT = Float64
Expand Down Expand Up @@ -47,6 +46,16 @@ M_insol = 0.044 # molar mass of insol
ρ_insol = 1770.0 # density of insol
κ_insol = 0.0 # hygroscopicity of insol

# Schemes
ARG_scheme = CMT.ARG2000Type()
ML_scheme = AA.MLEmulatedAerosolActivation(
joinpath(
pkgdir(CloudMicrophysics),
"aerosol_activation_emulators",
"2modal_nn_machine_naive.jls",
),
)

function mass2vol(mass_mixing_ratios)
if length(mass_mixing_ratios) == 2
densities = (ρ_sulfate, ρ_insol)
Expand All @@ -61,7 +70,7 @@ end

# Abdul-Razzak and Ghan 2000
# https://doi.org/10.1029/1999JD901161
function make_ARG_figX(X)
function make_ARG_figX(X, scheme = ARG_scheme)
p1 = PL.plot()
p2 = PL.plot()

Expand Down Expand Up @@ -176,19 +185,39 @@ function make_ARG_figX(X)
end
AD = AM.AerosolDistribution((paper_mode_1, paper_mode_2))
act_frac1[it] =
AA.N_activated_per_mode(param_set, AD, T, p, w, q)[1] / N_1
AA.N_activated_per_mode(
param_set,
scheme,
AD,
T,
p,
w,
q,
)[1] / N_1
act_frac2[it] =
AA.N_activated_per_mode(param_set, AD, T, p, w, q)[2] / N2i
AA.N_activated_per_mode(
param_set,
scheme,
AD,
T,
p,
w,
q,
)[2] / N2i
global it += 1
end

x1_obs = Fig1_x_obs
y1_obs = Fig1_y_obs
x1_PySDM = Fig1_x_PySDM
y1_PySDM = Fig1_y_PySDM
x1_param = Fig1_x_param
y1_param = Fig1_y_param

x2_obs = Fig1_x_obs
y2_obs = Fig1_y_obs
x2_PySDM = Fig1_x_PySDM
y2_PySDM = Fig1_y_PySDM
x2_param = Fig1_x_param
y2_param = Fig1_y_param

Expand Down Expand Up @@ -230,19 +259,39 @@ function make_ARG_figX(X)
end
AD = AM.AerosolDistribution((paper_mode_1, paper_mode_2))
act_frac1[it] =
AA.N_activated_per_mode(param_set, AD, T, p, w, q)[1] / N_1
AA.N_activated_per_mode(
param_set,
scheme,
AD,
T,
p,
w,
q,
)[1] / N_1
act_frac2[it] =
AA.N_activated_per_mode(param_set, AD, T, p, w, q)[2] / N2i
AA.N_activated_per_mode(
param_set,
scheme,
AD,
T,
p,
w,
q,
)[2] / N2i
global it += 1
end

x1_obs = Fig2a_x_obs
y1_obs = Fig2a_y_obs
x1_PySDM = Fig2a_x_PySDM
y1_PySDM = Fig2a_y_PySDM
x1_param = Fig2a_x_param
y1_param = Fig2a_y_param

x2_obs = Fig2b_x_obs
y2_obs = Fig2b_y_obs
x2_PySDM = Fig2b_x_PySDM
y2_PySDM = Fig2b_y_PySDM
x2_param = Fig2b_x_param
y2_param = Fig2b_y_param

Expand Down Expand Up @@ -286,19 +335,39 @@ function make_ARG_figX(X)
end
AD = AM.AerosolDistribution((paper_mode_1, paper_mode_2))
act_frac1[it] =
AA.N_activated_per_mode(param_set, AD, T, p, w, q)[1] / N_1
AA.N_activated_per_mode(
param_set,
scheme,
AD,
T,
p,
w,
q,
)[1] / N_1
act_frac2[it] =
AA.N_activated_per_mode(param_set, AD, T, p, w, q)[2] / N_2
AA.N_activated_per_mode(
param_set,
scheme,
AD,
T,
p,
w,
q,
)[2] / N_2
global it += 1
end

x1_obs = Fig3a_x_obs
y1_obs = Fig3a_y_obs
x1_PySDM = Fig3a_x_PySDM
y1_PySDM = Fig3a_y_PySDM
x1_param = Fig3a_x_param
y1_param = Fig3a_y_param

x2_obs = Fig3b_x_obs
y2_obs = Fig3b_y_obs
x2_PySDM = Fig3b_x_PySDM
y2_PySDM = Fig3b_y_PySDM
x2_param = Fig3b_x_param
y2_param = Fig3b_y_param

Expand Down Expand Up @@ -339,19 +408,39 @@ function make_ARG_figX(X)
end
AD = AM.AerosolDistribution((paper_mode_1, paper_mode_2))
act_frac1[it] =
AA.N_activated_per_mode(param_set, AD, T, p, w, q)[1] / N_1
AA.N_activated_per_mode(
param_set,
scheme,
AD,
T,
p,
w,
q,
)[1] / N_1
act_frac2[it] =
AA.N_activated_per_mode(param_set, AD, T, p, w, q)[2] / N_2
AA.N_activated_per_mode(
param_set,
scheme,
AD,
T,
p,
w,
q,
)[2] / N_2
global it += 1
end

x1_obs = Fig4a_x_obs
y1_obs = Fig4a_y_obs
x1_PySDM = Fig4a_x_PySDM
y1_PySDM = Fig4a_y_PySDM
x1_param = Fig4a_x_param
y1_param = Fig4a_y_param

x2_obs = Fig4b_x_obs
y2_obs = Fig4b_y_obs
x2_PySDM = Fig4b_x_PySDM
y2_PySDM = Fig4b_y_PySDM
x2_param = Fig4b_x_param
y2_param = Fig4b_y_param

Expand Down Expand Up @@ -394,19 +483,39 @@ function make_ARG_figX(X)

for wi in w
act_frac1[it] =
AA.N_activated_per_mode(param_set, AD, T, p, wi, q)[1] / N_1
AA.N_activated_per_mode(
param_set,
scheme,
AD,
T,
p,
wi,
q,
)[1] / N_1
act_frac2[it] =
AA.N_activated_per_mode(param_set, AD, T, p, wi, q)[2] / N_2
AA.N_activated_per_mode(
param_set,
scheme,
AD,
T,
p,
wi,
q,
)[2] / N_2
global it += 1
end

x1_obs = Fig5a_x_obs
y1_obs = Fig5a_y_obs
x1_PySDM = Fig5a_x_PySDM
y1_PySDM = Fig5a_y_PySDM
x1_param = Fig5a_x_param
y1_param = Fig5a_y_param

x2_obs = Fig5b_x_obs
y2_obs = Fig5b_y_obs
x2_PySDM = Fig5b_x_PySDM
y2_PySDM = Fig5b_y_PySDM
x2_param = Fig5b_x_param
y2_param = Fig5b_y_param

Expand Down Expand Up @@ -449,6 +558,14 @@ function make_ARG_figX(X)
)
PL.scatter!(p2, x2_obs, y2_obs, markercolor = :black)
PL.plot!(p2, x2_param, y2_param, linecolor = :black)
PL.scatter!(
p1,
x1_PySDM,
y1_PySDM,
markercolor = :green,
label = "PySDM observations",
)
PL.scatter!(p2, x2_PySDM, y2_PySDM, markercolor = :green)
end
end

Expand Down
Loading