Skip to content

Commit

Permalink
wip
Browse files Browse the repository at this point in the history
  • Loading branch information
@trontrytel
trontrytel committed Mar 1, 2024
1 parent d5c26d4 commit 79a0b14
Showing 1 changed file with 158 additions and 183 deletions.
341 changes: 158 additions & 183 deletions ml_test/PreprocessAerosolData.jl
View file
Original file line number Diff line number Diff line change
@@ -1,13 +1,21 @@
import DataFrames as DF
using DataFramesMeta
#module PreprocessAerosolData

import CLIMAParameters
import Thermodynamics as TD
import CloudMicrophysics as CM

import CloudMicrophysics.AerosolActivation as AA
import CloudMicrophysics.AerosolModel as AM
import CloudMicrophysics.Parameters as CMP
import CloudMicrophysics.Parameters as CMP
import CloudMicrophysics.Common as CO

import ClimaParams as CP
import Thermodynamics as TD
import DataFrames as DF
using DataFramesMeta

const FT = Float64
#const APS = CMP.ParametersType

function get_num_modes(df::DataFrame)
i = 1
while true
Expand All @@ -28,170 +36,151 @@ function get_num_modes(data_row::NamedTuple)
end
end

function convert_to_ARG_params(data_row::NamedTuple)
#TODO
FT = Float64
tps = TD.Parameters.ThermodynamicsParameters(FT)

num_modes = get_num_modes(data_row)
@assert num_modes > 0
mode_Ns = []
mode_means = []
mode_stdevs = []
mode_kappas = []
w = data_row.velocity
T = data_row.initial_temperature
p = data_row.initial_pressure
for i in 1:num_modes
push!(mode_Ns, data_row[Symbol("mode_$(i)_N")])
push!(mode_means, data_row[Symbol("mode_$(i)_mean")])
push!(mode_stdevs, data_row[Symbol("mode_$(i)_stdev")])
push!(mode_kappas, data_row[Symbol("mode_$(i)_kappa")])
end
ad = AM.AerosolDistribution(
Tuple(
AM.Mode_κ(
mode_means[i],
mode_stdevs[i],
mode_Ns[i],
FT(1),
FT(1),
FT(0),
mode_kappas[i],
1,
) for i in 1:num_modes
),
)
pv0 = TD.saturation_vapor_pressure(tps, T, TD.Liquid())
vapor_mix_ratio =
pv0 / TD.Parameters.molmass_ratio(tps) / (p - pv0)
q_vap = vapor_mix_ratio / (vapor_mix_ratio + 1)
q = TD.PhasePartition(q_vap, FT(0), FT(0))
return (; ad, T, p, w, q, mode_Ns)
end

function convert_to_ARG_intermediates(
data_row::NamedTuple,
)
num_modes = get_num_modes(data_row)
@assert num_modes > 0

# TODO
FT = Float64
tps = TD.Parameters.ThermodynamicsParameters(FT)
aip = CMP.AirProperties(FT)
ap = CMP.AerosolActivationParameters(FT)

(; ad, T, p, w, q) = convert_to_ARG_params(data_row, tps)

ϵ::FT = 1 / TD.Parameters.molmass_ratio(tps)
R_m::FT = TD.gas_constant_air(tps, q)
cp_m::FT = TD.cp_m(tps, q)

L::FT = TD.latent_heat_vapor(tps, T)
p_vs::FT = TD.saturation_vapor_pressure(tps, T, TD.Liquid())
G::FT = CO.G_func(aip, tps, T, TD.Liquid()) / ap.ρ_w

α::FT = L * ap.g * ϵ / R_m / cp_m / T^2 - ap.g / R_m / T
γ::FT = R_m * T / ϵ / p_vs + ϵ * L^2 / cp_m / T / p

A::FT = AA.coeff_of_curvature(ap, T)
ζ::FT = 2 * A / 3 * sqrt(α * w / G)

Sm = AA.critical_supersaturation(ap, ad, T)
η = [
(α * w / G)^FT(3 / 2) / (FT(2 * pi) * ap.ρ_w * γ * ad.Modes[i].N) for i in 1:num_modes
]

per_mode_intermediates = [
(;
Symbol("mode_$(i)_stdev") => ad.Modes[i].stdev,
Symbol("mode_$(i)_η") => η[i],
Symbol("mode_$(i)_Sm") => Sm[i],
) for i in 1:num_modes
]
return merge(reduce(merge, per_mode_intermediates), (; ζ))
end

function get_ARG_S_max(data_row::NamedTuple)
(; ad, T, p, w, q) = convert_to_ARG_params(data_row)
# TODO
FT = Float64
tps = TD.Parameters.ThermodynamicsParameters(FT)
aip = CMP.AirProperties(FT)
ap = CMP.AerosolActivationParameters(FT)

max_supersaturation =
AA.max_supersaturation(ap, ad, aip, tps, T, p, w, q)
return max_supersaturation
end

function get_ARG_S_max(X::DataFrame)
return get_ARG_S_max.(NamedTuple.(eachrow(X)))
end

function get_ARG_S_crit(data_row::NamedTuple)
#TODO
FT = Float64
tps = TD.Parameters.ThermodynamicsParameters(FT)
ap = CMP.AerosolActivationParameters(FT)
(; ad, T) = convert_to_ARG_params(data_row, tps)
return AA.critical_supersaturation(ap, ad, T)
end

function get_ARG_S_crit(X::DataFrame)
return get_ARG_S_crit.(NamedTuple.(eachrow(X)))
end

function get_ARG_act_N(data_row::NamedTuple, S_max = nothing)

#TODO
FT = Float64
tps = TD.Parameters.ThermodynamicsParameters(FT)
(; ad, T, p, w, q) = convert_to_ARG_params(data_row, tps)

if S_max === nothing
return collect(
AA.N_activated_per_mode(ap, ad, aip, tps, T, p, w, q)
)
else
return collect(
AA.N_activated_per_mode(ap, ad, aip, tps, T, p, w, q, S_max)
)
end
end

function get_ARG_act_N(data_row::NamedTuple, S_max = nothing)
return get_ARG_act_N(data_row, S_max)
end

function get_ARG_act_N(X::DataFrame, S_max = nothing)
return transpose(
hcat(get_ARG_act_N.(NamedTuple.(eachrow(X)), S_max)...),
)
end

function get_ARG_act_frac(data_row::NamedTuple, S_max = nothing)
(; mode_Ns) = convert_to_ARG_params(data_row)
return get_ARG_act_N(data_row, S_max) ./ mode_Ns
end

function get_ARG_act_frac(X::DataFrame, S_max = nothing)
return transpose(
hcat(get_ARG_act_frac.(NamedTuple.(eachrow(X)), S_max)...),
)
end

function preprocess_aerosol_data(X::DataFrame, add_ARG_act_frac::Bool)
#function convert_to_ARG_params(data_row::NamedTuple, tps)
#
# num_modes = get_num_modes(data_row)
# @assert num_modes > 0
# mode_Ns = []
# mode_means = []
# mode_stdevs = []
# mode_kappas = []
# w = data_row.velocity
# T = data_row.initial_temperature
# p = data_row.initial_pressure
# for i in 1:num_modes
# push!(mode_Ns, data_row[Symbol("mode_$(i)_N")])
# push!(mode_means, data_row[Symbol("mode_$(i)_mean")])
# push!(mode_stdevs, data_row[Symbol("mode_$(i)_stdev")])
# push!(mode_kappas, data_row[Symbol("mode_$(i)_kappa")])
# end
# ad = AM.AerosolDistribution(
# Tuple(
# AM.Mode_κ(
# mode_means[i],
# mode_stdevs[i],
# mode_Ns[i],
# FT(1),
# FT(1),
# FT(0),
# mode_kappas[i],
# 1,
# ) for i in 1:num_modes
# ),
# )
# pv0 = TD.saturation_vapor_pressure(tps, T, TD.Liquid())
# vapor_mix_ratio = pv0 / TD.Parameters.molmass_ratio(tps) / (p - pv0)
# q_vap = vapor_mix_ratio / (vapor_mix_ratio + 1)
# q = TD.PhasePartition(q_vap, FT(0), FT(0))
# return (; ad, T, p, w, q, mode_Ns)
#end
#
#function convert_to_ARG_params(data_row::NamedTuple)
# return convert_to_ARG_params(data_row, default_param_set)
#end
#
#function get_ARG_S_max(data_row::NamedTuple, param_set::APS)
# (; ad, T, p, w, q) = convert_to_ARG_params(data_row, param_set)
# max_supersaturation =
# AA.max_supersaturation(param_set, CT.ARG2000Type(), ad, T, p, w, q)
# return max_supersaturation
#end
#
#function get_ARG_S_max(data_row::NamedTuple)
# return get_ARG_S_max(data_row, default_param_set)
#end
#
#function get_ARG_S_max(X::DataFrame, param_set::APS)
# return get_ARG_S_max.(NamedTuple.(eachrow(X)), param_set)
#end
#
#function get_ARG_S_max(X::DataFrame)
# return get_ARG_S_max(X, default_param_set)
#end
#
#function get_ARG_S_crit(data_row::NamedTuple, param_set::APS)
# (; ad, T) = convert_to_ARG_params(data_row, param_set)
# return AA.critical_supersaturation(param_set, ad, T)
#end
#
#function get_ARG_S_crit(data_row::NamedTuple)
# return get_ARG_S_crit(data_row, default_param_set)
#end
#
#function get_ARG_S_crit(X::DataFrame, param_set::APS)
# return get_ARG_S_crit.(NamedTuple.(eachrow(X)), param_set)
#end
#
#function get_ARG_S_crit(X::DataFrame)
# return get_ARG_S_crit(X, default_param_set)
#end
#
#function get_ARG_act_N(data_row::NamedTuple, param_set::APS, S_max = nothing)
# (; ad, T, p, w, q) = convert_to_ARG_params(data_row, param_set)
# if S_max === nothing
# return collect(
# AA.N_activated_per_mode(
# param_set,
# CT.ARG2000Type(),
# ad,
# T,
# p,
# w,
# q,
# ),
# )
# else
# critical_supersaturation = AA.critical_supersaturation(param_set, ad, T)
# return collect(
# AA.N_activated_per_mode(
# param_set,
# CT.ARG2000Type(),
# ad,
# T,
# p,
# w,
# q,
# S_max,
# critical_supersaturation,
# ),
# )
# end
#end
#
#function get_ARG_act_N(data_row::NamedTuple, S_max = nothing)
# return get_ARG_act_N(data_row, default_param_set, S_max)
#end
#
#function get_ARG_act_N(X::DataFrame, param_set::APS, S_max = nothing)
# return transpose(
# hcat(get_ARG_act_N.(NamedTuple.(eachrow(X)), param_set, S_max)...),
# )
#end
#
#function get_ARG_act_N(X::DataFrame, S_max = nothing)
# return get_ARG_act_N(X, default_param_set, S_max)
#end
#
#function get_ARG_act_frac(data_row::NamedTuple, param_set::APS, S_max = nothing)
# (; mode_Ns) = convert_to_ARG_params(data_row, param_set)
# return get_ARG_act_N(data_row, param_set, S_max) ./ mode_Ns
#end
#
#function get_ARG_act_frac(data_row::NamedTuple, S_max = nothing)
# return get_ARG_act_frac(data_row, default_param_set, S_max)
#end
#
#function get_ARG_act_frac(X::DataFrame, param_set::APS, S_max = nothing)
# return transpose(
# hcat(get_ARG_act_frac.(NamedTuple.(eachrow(X)), param_set, S_max)...),
# )
#end
#
#function get_ARG_act_frac(X::DataFrame, S_max = nothing)
# return get_ARG_act_frac(X, default_param_set, S_max)
#end

function preprocess_aerosol_data(X::DataFrame)
num_modes = get_num_modes(X)
if add_ARG_act_frac
X = DF.transform(
X,
AsTable(All()) =>
ByRow(x -> get_ARG_act_frac(x)) =>
[Symbol("mode_$(i)_ARG_act_frac") for i in 1:num_modes],
)
end
for i in 1:num_modes
X = DF.transform(
X,
Expand All @@ -207,26 +196,12 @@ function preprocess_aerosol_data(X::DataFrame, add_ARG_act_frac::Bool)
return X
end

function preprocess_aerosol_data_standard(X::DataFrame)
return preprocess_aerosol_data(X, false)
end

function preprocess_aerosol_data_with_ARG_act_frac(X::DataFrame)
return preprocess_aerosol_data(X, true)
end

function preprocess_aerosol_data_with_ARG_intermediates(X::DataFrame)
return DF.DataFrame(
convert_to_ARG_intermediates.(
NamedTuple.(eachrow(X))
)
)
end

function target_transform(act_frac)
return @. atanh(2.0 * 0.99 * (act_frac - 0.5))
end

function inverse_target_transform(transformed_act_frac)
return @. (1.0 / (2.0 * 0.99)) * tanh(transformed_act_frac) + 0.5
end

#end # module

0 comments on commit 79a0b14

Please sign in to comment.