Merge #100

100: Parameter Box integration with InternalClimaParams r=charleskawczynski a=odunbar Co-authored-by: odunbar <[email protected]> Co-authored-by: Charles Kawczynski <[email protected]>
CliMA · Jun 17, 2022 · a7cadd6 · a7cadd6 · charleskawczynski · Jun 17, 2022
2 parents 92184df + 18cab0b
commit a7cadd6
Show file tree

Hide file tree

Showing 18 changed files with 436 additions and 334 deletions.
diff --git a/Project.toml b/Project.toml
@@ -1,17 +1,15 @@
 name = "Thermodynamics"
 uuid = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c"
 authors = ["Climate Modeling Alliance"]
-version = "0.8.0"
+version = "0.9.0"
 
 [deps]
-CLIMAParameters = "6eacf6c3-8458-43b9-ae03-caf5306d3d53"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 RootSolvers = "7181ea78-2dcb-4de3-ab41-2b8ab5a31e74"
 
 [compat]
-CLIMAParameters = "0.6"
 DocStringExtensions = "0.8.1, 0.9"
 KernelAbstractions = "0.7.2"
 RootSolvers = "0.2, 0.3"

diff --git a/docs/src/TemperatureProfiles.md b/docs/src/TemperatureProfiles.md
@@ -21,10 +21,12 @@ to one of the temperature profile constructors.
 import Thermodynamics as TD
 import Plots
 import CLIMAParameters as CP
-import Thermodynamics.InternalClimaParams as ICP
+import Thermodynamics.Parameters as TP
 FT = Float64;
-struct EarthParameterSet <: CP.AbstractEarthParameterSet end;
-const param_set = EarthParameterSet();
+toml_dict = CP.create_toml_dict(FT; dict_type = "alias")
+aliases = string.(fieldnames(TP.ThermodynamicsParameters))
+param_pairs = CP.get_parameter_values!(toml_dict, aliases, "Thermodynamics")
+param_set = TP.ThermodynamicsParameters{FT}(; param_pairs...)
 z = range(FT(0), stop = FT(2.5e4), length = 50);
 
 isothermal = TD.TemperatureProfiles.IsothermalProfile(param_set, FT);
@@ -46,10 +48,12 @@ Plots.savefig("isothermal.svg");
 import Thermodynamics as TD
 import Plots
 import CLIMAParameters as CP
-import Thermodynamics.InternalClimaParams as ICP
+import Thermodynamics.Parameters as TP
 FT = Float64;
-struct EarthParameterSet <: CP.AbstractEarthParameterSet end;
-const param_set = EarthParameterSet();
+toml_dict = CP.create_toml_dict(FT; dict_type = "alias")
+aliases = string.(fieldnames(TP.ThermodynamicsParameters))
+param_pairs = CP.get_parameter_values!(toml_dict, aliases, "Thermodynamics")
+param_set = TP.ThermodynamicsParameters{FT}(; param_pairs...)
 z = range(FT(0), stop = FT(2.5e4), length = 50);
 
 decaying = TD.TemperatureProfiles.DecayingTemperatureProfile{FT}(param_set);
@@ -70,10 +74,12 @@ Plots.savefig("decaying.svg")
 import Thermodynamics as TD
 import Plots
 import CLIMAParameters as CP
-import Thermodynamics.InternalClimaParams as ICP
+import Thermodynamics.Parameters as TP
 FT = Float64;
-struct EarthParameterSet <: CP.AbstractEarthParameterSet end;
-const param_set = EarthParameterSet();
+toml_dict = CP.create_toml_dict(FT; dict_type = "alias")
+aliases = string.(fieldnames(TP.ThermodynamicsParameters))
+param_pairs = CP.get_parameter_values!(toml_dict, aliases, "Thermodynamics")
+param_set = TP.ThermodynamicsParameters{FT}(; param_pairs...)
 z = range(FT(0), stop = FT(2.5e4), length = 50);
 
 dry_adiabatic = TD.TemperatureProfiles.DryAdiabaticProfile{FT}(param_set);

diff --git a/docs/src/TestedProfiles.md b/docs/src/TestedProfiles.md
@@ -8,10 +8,12 @@ Thermodynamics.jl is tested using a set of profiles specified in `src/TestedProf
 import Thermodynamics as TD
 import Plots
 import CLIMAParameters as CP
-import Thermodynamics.InternalClimaParams as ICP
-FT = Float32;
-struct EarthParameterSet <: CP.AbstractEarthParameterSet end;
-const param_set = EarthParameterSet();
+import Thermodynamics.Parameters as TP
+FT = Float64;
+toml_dict = CP.create_toml_dict(FT; dict_type = "alias")
+aliases = string.(fieldnames(TP.ThermodynamicsParameters))
+param_pairs = CP.get_parameter_values!(toml_dict, aliases, "Thermodynamics")
+param_set = TP.ThermodynamicsParameters{FT}(; param_pairs...)
 
 profiles = TD.TestedProfiles.PhaseDryProfiles(param_set, Array{FT});
 (;T, ρ, z) = profiles
@@ -28,10 +30,12 @@ Plots.savefig("tested_profiles_dry.svg");
 import Thermodynamics as TD
 import Plots
 import CLIMAParameters as CP
-import Thermodynamics.InternalClimaParams as ICP
-FT = Float32;
-struct EarthParameterSet <: CP.AbstractEarthParameterSet end;
-const param_set = EarthParameterSet();
+import Thermodynamics.Parameters as TP
+FT = Float64;
+toml_dict = CP.create_toml_dict(FT; dict_type = "alias")
+aliases = string.(fieldnames(TP.ThermodynamicsParameters))
+param_pairs = CP.get_parameter_values!(toml_dict, aliases, "Thermodynamics")
+param_set = TP.ThermodynamicsParameters{FT}(; param_pairs...)
 
 profiles = TD.TestedProfiles.PhaseEquilProfiles(param_set, Array{FT});
 (;T, ρ, q_tot, z) = profiles

diff --git a/docs/src/saturation_adjustment.jl b/docs/src/saturation_adjustment.jl
@@ -4,12 +4,15 @@ import Plots
 
 FT = Float64;
 import Thermodynamics as TD
-const ICP = TD.InternalClimaParams
+const TP = TD.Parameters
 import CLIMAParameters as CP
-struct EarthParameterSet <: CP.AbstractEarthParameterSet end;
-const param_set = EarthParameterSet();
 
 TD.print_warning() = false
+toml_dict = CP.create_toml_dict(FT; dict_type = "alias")
+aliases = string.(fieldnames(TP.ThermodynamicsParameters))
+param_pairs = CP.get_parameter_values!(toml_dict, aliases, "Thermodynamics")
+const param_set = TP.ThermodynamicsParameters{FT}(; param_pairs...)
+
 profiles = TD.TestedProfiles.PhaseEquilProfiles(param_set, Array{FT});
 (; ρ, q_tot) = profiles
 T_true = profiles.T

diff --git a/perf/common.jl b/perf/common.jl
@@ -1,6 +1,7 @@
 using Test
 import Thermodynamics
 const TD = Thermodynamics
+const TP = TD.Parameters
 
 import UnPack
 import BenchmarkTools
@@ -10,11 +11,13 @@ const RS = RootSolvers
 
 import CLIMAParameters
 const CP = CLIMAParameters
+const FT = Float64
+toml_dict = CP.create_toml_dict(FT; dict_type = "alias")
+aliases = string.(fieldnames(TP.ThermodynamicsParameters))
+param_pairs = CP.get_parameter_values!(toml_dict, aliases, "Thermodynamics")
+const param_set = TP.ThermodynamicsParameters{FT}(; param_pairs...)
 
-struct EarthParameterSet <: CP.AbstractEarthParameterSet end
-const param_set = EarthParameterSet()
-
-ArrayType = Array{Float64}
+ArrayType = Array{FT}
 profiles = TD.TestedProfiles.PhaseEquilProfiles(param_set, ArrayType)
 UnPack.@unpack e_int, T, ρ, p, θ_liq_ice, q_tot = profiles
 

diff --git a/perf/common_micro_bm.jl b/perf/common_micro_bm.jl
@@ -3,7 +3,7 @@ import StatsBase
 import PrettyTables
 import OrderedCollections
 const TD = Thermodynamics
-const ICP = TD.InternalClimaParams
+const TP = TD.Parameters
 using JET
 using Test
 
@@ -12,16 +12,18 @@ import BenchmarkTools
 
 import CLIMAParameters
 const CP = CLIMAParameters
-
-struct EarthParameterSet <: CP.AbstractEarthParameterSet end
-const param_set = EarthParameterSet()
+const FT = Float64
+toml_dict = CP.create_toml_dict(FT; dict_type = "alias")
+aliases = string.(fieldnames(TP.ThermodynamicsParameters))
+param_pairs = CP.get_parameter_values!(toml_dict, aliases, "Thermodynamics")
+const param_set = TP.ThermodynamicsParameters{FT}(; param_pairs...)
 
 #####
 ##### Finding indexes in profiles satisfying certain conditions
 #####
 
 function find_freezing_index(profiles)
-    i = findfirst(T -> T === ICP.T_freeze(param_set), profiles.T)
+    i = findfirst(T -> T === TP.T_freeze(param_set), profiles.T)
     isnothing(i) && error("Freezing index not found")
     return i
 end

diff --git a/src/InternalClimaParams.jl b/src/InternalClimaParams.jl
diff --git a/src/Parameters.jl b/src/Parameters.jl
@@ -0,0 +1,68 @@
+module Parameters
+
+"""
+    ThermodynamicsParameters
+
+Parameters for Thermodynamics.jl.
+
+# Example
+```
+import CLIMAParameters as CP
+import Thermodynamics.Parameters as TP
+FT = Float64;
+toml_dict = CP.create_toml_dict(FT; dict_type = "alias");
+aliases = string.(fieldnames(TP.ThermodynamicsParameters));
+param_pairs = CP.get_parameter_values!(toml_dict, aliases, "Thermodynamics");
+param_set = TP.ThermodynamicsParameters{FT}(; param_pairs...);
+```
+"""
+Base.@kwdef struct ThermodynamicsParameters{FT}
+    T_0::FT
+    MSLP::FT
+    cp_v::FT
+    cp_l::FT
+    cp_i::FT
+    LH_v0::FT
+    LH_s0::FT
+    press_triple::FT
+    T_triple::FT
+    T_freeze::FT
+    T_min::FT
+    T_max::FT
+    entropy_reference_temperature::FT
+    entropy_dry_air::FT
+    entropy_water_vapor::FT
+    kappa_d::FT
+    gas_constant::FT
+    molmass_dryair::FT
+    molmass_water::FT
+    T_surf_ref::FT
+    T_min_ref::FT
+    grav::FT
+    T_icenuc::FT
+    pow_icenuc::FT
+end
+
+const ATP = ThermodynamicsParameters
+
+Base.broadcastable(ps::ATP) = Ref(ps)
+
+# wrappers
+for fn in fieldnames(ATP)
+    @eval $(fn)(ps::ATP) = ps.$(fn)
+end
+
+# Derived parameters
+R_d(ps::ATP) = ps.gas_constant / ps.molmass_dryair
+R_v(ps::ATP) = ps.gas_constant / ps.molmass_water
+molmass_ratio(ps::ATP) = ps.molmass_dryair / ps.molmass_water
+LH_f0(ps::ATP) = ps.LH_s0 - ps.LH_v0
+e_int_v0(ps::ATP) = ps.LH_v0 - R_v(ps) * ps.T_0
+e_int_i0(ps::ATP) = LH_f0(ps)
+cp_d(ps::ATP) = R_d(ps) / ps.kappa_d
+cv_d(ps::ATP) = cp_d(ps) - R_d(ps)
+cv_v(ps::ATP) = ps.cp_v - R_v(ps)
+cv_l(ps::ATP) = ps.cp_l
+cv_i(ps::ATP) = ps.cp_i
+
+end
diff --git a/src/TemperatureProfiles.jl b/src/TemperatureProfiles.jl
@@ -6,12 +6,9 @@ const DSE = DocStringExtensions
 export TemperatureProfile,
     IsothermalProfile, DecayingTemperatureProfile, DryAdiabaticProfile
 
-import CLIMAParameters
-const CP = CLIMAParameters
-const APS = CP.AbstractParameterSet
-
-import ..InternalClimaParams
-const ICP = InternalClimaParams
+import ..Parameters
+const TP = Parameters
+const APS = TP.ThermodynamicsParameters
 
 """
     TemperatureProfile
@@ -37,7 +34,7 @@ IsothermalProfile(param_set::APS, T_virt::FT) where {FT} =
     DecayingTemperatureProfile{FT}(param_set, T_virt, T_virt)
 
 function IsothermalProfile(param_set::APS, ::Type{FT}) where {FT}
-    T_virt = FT(ICP.T_surf_ref(param_set))
+    T_virt = FT(TP.T_surf_ref(param_set))
     return DecayingTemperatureProfile{FT}(param_set, T_virt, T_virt)
 end
 
@@ -58,8 +55,8 @@ struct DryAdiabaticProfile{FT} <: TemperatureProfile{FT}
     T_min_ref::FT
     function DryAdiabaticProfile{FT}(
         param_set::APS,
-        T_surface::FT = FT(ICP.T_surf_ref(param_set)),
-        _T_min_ref::FT = FT(ICP.T_min_ref(param_set)),
+        T_surface::FT = FT(TP.T_surf_ref(param_set)),
+        _T_min_ref::FT = FT(TP.T_min_ref(param_set)),
     ) where {FT}
         return new{FT}(T_surface, _T_min_ref)
     end
@@ -77,10 +74,10 @@ to be greater than or equal to `profile.T_min_ref`.
 """
 function (profile::DryAdiabaticProfile)(param_set::APS, z::FT) where {FT}
 
-    R_d::FT = ICP.R_d(param_set)
-    cp_d::FT = ICP.cp_d(param_set)
-    grav::FT = ICP.grav(param_set)
-    MSLP::FT = ICP.MSLP(param_set)
+    R_d::FT = TP.R_d(param_set)
+    cp_d::FT = TP.cp_d(param_set)
+    grav::FT = TP.grav(param_set)
+    MSLP::FT = TP.MSLP(param_set)
 
     # Temperature
     Γ = grav / cp_d
@@ -120,20 +117,19 @@ struct DecayingTemperatureProfile{FT} <: TemperatureProfile{FT}
     H_t::FT
     function DecayingTemperatureProfile{FT}(
         param_set::APS,
-        _T_virt_surf::FT = FT(ICP.T_surf_ref(param_set)),
-        _T_min_ref::FT = FT(ICP.T_min_ref(param_set)),
-        H_t::FT = FT(ICP.R_d(param_set)) * _T_virt_surf /
-                  FT(ICP.grav(param_set)),
+        _T_virt_surf::FT = FT(TP.T_surf_ref(param_set)),
+        _T_min_ref::FT = FT(TP.T_min_ref(param_set)),
+        H_t::FT = FT(TP.R_d(param_set)) * _T_virt_surf / FT(TP.grav(param_set)),
     ) where {FT}
         return new{FT}(_T_virt_surf, _T_min_ref, H_t)
     end
 end
 
 
 function (profile::DecayingTemperatureProfile)(param_set::APS, z::FT) where {FT}
-    R_d::FT = ICP.R_d(param_set)
-    grav::FT = ICP.grav(param_set)
-    MSLP::FT = ICP.MSLP(param_set)
+    R_d::FT = TP.R_d(param_set)
+    grav::FT = TP.grav(param_set)
+    MSLP::FT = TP.MSLP(param_set)
 
     # Scale height for surface temperature
     H_sfc = R_d * profile.T_virt_surf / grav