Merge #128

128: Use dispatch for entr-detr, + DocStringExtensions r=charleskawczynski a=charleskawczynski This PR - changes the mutable `entr_detr_model` type to a immutable `MyEntrModel` type - Adds a `EntrResults` type so that the results can be more safely unpacked ( not order-dependent ) - Adds doc strings via `DocStringExtensions` Co-authored-by: Charles Kawczynski <[email protected]>
CliMA · Aug 16, 2021 · 1104d43 · 1104d43
2 parents 2fa9779 + 3e87cde
commit 1104d43
Show file tree

Hide file tree

Showing 8 changed files with 105 additions and 89 deletions.
diff --git a/Project.toml b/Project.toml
@@ -10,6 +10,7 @@ CLIMAParameters = "6eacf6c3-8458-43b9-ae03-caf5306d3d53"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 Dierckx = "39dd38d3-220a-591b-8e3c-4c3a8c710a94"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
+DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
 Glob = "c27321d9-0574-5035-807b-f59d2c89b15c"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"

diff --git a/docs/Project.toml b/docs/Project.toml
@@ -1,6 +1,7 @@
 [deps]
 ArgParse = "c7e460c6-2fb9-53a9-8c5b-16f535851c63"
 CLIMAParameters = "6eacf6c3-8458-43b9-ae03-caf5306d3d53"
+DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244"
 Glob = "c27321d9-0574-5035-807b-f59d2c89b15c"

diff --git a/src/TurbulenceConvection.jl b/src/TurbulenceConvection.jl
@@ -3,6 +3,7 @@ module TurbulenceConvection
 using StaticArrays
 using StatsBase
 using LinearAlgebra
+import DocStringExtensions
 
 import Dierckx
 import Statistics

diff --git a/src/Turbulence_PrognosticTKE.jl b/src/Turbulence_PrognosticTKE.jl
@@ -895,27 +895,17 @@ function get_env_covar_from_GMV(
 end
 
 function compute_updraft_closures(self::EDMF_PrognosticTKE, GMV::GridMeanVariables, Case::CasesBase)
-    εδ_model = entr_detr_model()
-    sa = eos_struct()
-    quadrature_order = 3
     grid = get_grid(self)
     param_set = parameter_set(self)
-
-    upd_cloud_diagnostics(self.UpdVar, reference_state(self))
-
-    # TODO: merge this with grid.cinterior
-    kc_surf = kc_surface(grid)
-    kc_toa = kc_top_of_atmos(grid)
-    cinterior = kc_surf:kc_toa
-
     ref_state = reference_state(self)
-    l = zeros(2)
+
+    upd_cloud_diagnostics(self.UpdVar, ref_state)
 
     @inbounds for k in real_center_indicies(grid)
         @inbounds for i in xrange(self.n_updrafts)
             # entrainment
-            w_upd = interpf2c(self.UpdVar.W.values, grid, k, i)
             if self.UpdVar.Area.values[i, k] > 0.0
+                w_upd = interpf2c(self.UpdVar.W.values, grid, k, i)
                 # compute dMdz at half levels
                 gmv_w_k = interpf2c(GMV.W.values, grid, k)
                 gmv_w_km = interpf2c(GMV.W.values, grid, k - 1)
@@ -925,24 +915,24 @@ function compute_updraft_closures(self::EDMF_PrognosticTKE, GMV::GridMeanVariabl
                 dMdz = (M - Mm) * grid.dzi
                 w_min = 0.001
 
-                εδ_model.b_upd = self.UpdVar.B.values[i, k]
-                εδ_model.b_env = self.EnvVar.B.values[k]
-                εδ_model.w_upd = interpf2c(self.UpdVar.W.values, grid, k, i)
-                εδ_model.w_env = interpf2c(self.EnvVar.W.values, grid, k)
-                εδ_model.a_upd = self.UpdVar.Area.values[i, k]
-                εδ_model.a_env = self.EnvVar.Area.values[k]
-                εδ_model.ql_upd = self.UpdVar.QL.values[i, k]
-                εδ_model.ql_env = self.EnvVar.QL.values[k]
-                εδ_model.RH_upd = self.UpdVar.RH.values[i, k]
-                εδ_model.RH_env = self.EnvVar.RH.values[k]
-                εδ_model.M = M
-                εδ_model.dMdz = dMdz
-                εδ_model.tke = self.EnvVar.TKE.values[k]
-                εδ_model.n_up = self.n_updrafts
-                εδ_model.ρ = ref_state.rho0[k]
-                εδ_model.R_up = self.pressure_plume_spacing[i]
-
-                self.entr_sc[i, k], self.detr_sc[i, k], self.frac_turb_entr[i, k], self.horiz_K_eddy[i, k] = entr_detr(
+                εδ_model = MoistureDeficitEntr(;
+                    q_liq_up = self.UpdVar.QL.values[i, k],
+                    q_liq_en = self.EnvVar.QL.values[k],
+                    w_up = interpf2c(self.UpdVar.W.values, grid, k, i),
+                    w_en = interpf2c(self.EnvVar.W.values, grid, k),
+                    b_up = self.UpdVar.B.values[i, k],
+                    b_en = self.EnvVar.B.values[k],
+                    tke = self.EnvVar.TKE.values[k],
+                    dMdz = dMdz,
+                    M = M,
+                    a_up = self.UpdVar.Area.values[i, k],
+                    a_en = self.EnvVar.Area.values[k],
+                    R_up = self.pressure_plume_spacing[i],
+                    RH_up = self.UpdVar.RH.values[i, k],
+                    RH_en = self.EnvVar.RH.values[k],
+                )
+
+                er = entr_detr(
                     param_set,
                     w_min,
                     self.sorting_power,
@@ -955,6 +945,10 @@ function compute_updraft_closures(self::EDMF_PrognosticTKE, GMV::GridMeanVariabl
                     self.turbulent_entrainment_factor,
                     εδ_model,
                 )
+                self.entr_sc[i, k] = er.ε_dyn
+                self.detr_sc[i, k] = er.δ_dyn
+                self.frac_turb_entr[i, k] = er.ε_turb
+                self.horiz_K_eddy[i, k] = er.K_ε
             else
                 self.entr_sc[i, k] = 0.0
                 self.detr_sc[i, k] = 0.0

diff --git a/src/closures/entr_detr.jl b/src/closures/entr_detr.jl
@@ -1,4 +1,5 @@
 #### Entrainment-Detrainment kernels
+
 """
     entr_detr(
         param_set,
@@ -21,35 +22,35 @@ Cohen et al. (JAMES, 2020), given:
  - `β`: sorting power for moist mixing
  - `c_δ`: detrainment factor
  - `c_div`: divergence factor for bubble case (zero otherwise)
- - `c_μ`: logisitc function scale 
+ - `c_μ`: logisitc function scale
  - `c_μ0`: logisitc function timescale
  - `χ_upd`: updraft mixing fraction
  - `c_λ`: tke scale factor
  - `c_εt`: turbulent entrainment factor
- - `εδ_model`: entrainment detrainment model type
+ - `εδ_model`: a [`MoistureDeficitEntr`](@ref)
 """
-function entr_detr(param_set, w_min, β, c_δ, c_div, c_μ, c_μ0, χ_upd, c_λ, c_εt, εδ_model)
+function entr_detr(param_set, w_min, β, c_δ, c_div, c_μ, c_μ0, χ_upd, c_λ, c_εt, εδ_model::MoistureDeficitEntr)
 
-    l = zeros(2)
     c_ε = CPEDMF.c_ε(param_set)
 
-    if (εδ_model.ql_upd + εδ_model.ql_env) == 0.0
+    if (εδ_model.q_liq_up + εδ_model.q_liq_en) == 0.0
         c_δ = 0.0
     end
 
-    Δw = εδ_model.w_upd - εδ_model.w_env
+    Δw = εδ_model.w_up - εδ_model.w_en
     if Δw < 0.0
         Δw -= w_min
     else
         Δw += w_min
     end
 
-    Δb = (εδ_model.b_upd - εδ_model.b_env)
+    Δb = (εδ_model.b_up - εδ_model.b_en)
 
-    D_ε, D_δ, M_δ, M_ε = nondimensional_exchange_functions(c_ε, c_δ, c_μ, c_μ0, β, χ_upd, Δw, Δb, εδ_model)
+    D_ε, D_δ, M_δ, M_ε = nondimensional_exchange_functions(param_set, c_δ, c_μ, c_μ0, β, χ_upd, Δw, Δb, εδ_model)
 
-    l[1] = c_λ * abs(Δb / sqrt(εδ_model.tke + 1e-8))
-    l[2] = abs(Δb / Δw)
+    l_1 = c_λ * abs(Δb / sqrt(εδ_model.tke + 1e-8))
+    l_2 = abs(Δb / Δw)
+    l = (l_1, l_2)
     λ = lamb_smooth_minimum(l, 0.1, 0.0005)
 
     MdMdz = max(εδ_model.dMdz / max(εδ_model.M, 1e-12), 0.0)
@@ -58,15 +59,16 @@ function entr_detr(param_set, w_min, β, c_δ, c_div, c_μ, c_μ0, χ_upd, c_λ,
 
 
     # turbulent entrainment
-    k_ε = εδ_model.a_upd * c_εt * sqrt(max(εδ_model.tke, 0.0)) * εδ_model.R_up
-    if εδ_model.w_upd * εδ_model.a_upd > 0.0
-        ε_turb = (2.0 / εδ_model.R_up^2.0) * k_ε / (εδ_model.w_upd * εδ_model.a_upd)
+    K_ε = εδ_model.a_up * c_εt * sqrt(max(εδ_model.tke, 0.0)) * εδ_model.R_up
+    if εδ_model.w_up * εδ_model.a_up > 0.0
+        ε_turb = (2.0 / εδ_model.R_up^2.0) * K_ε / (εδ_model.w_up * εδ_model.a_up)
     else
         ε_turb = 0.0
     end
 
     ε_dyn = λ / Δw * (c_ε * D_ε + c_δ * M_ε) + MdMdz * c_div
     δ_dyn = λ / Δw * (c_ε * D_δ + c_δ * M_δ) + MdMdz * c_div
 
-    return ε_dyn, δ_dyn, ε_turb, k_ε
+
+    return EntrDetr(ε_dyn, δ_dyn, ε_turb, K_ε)
 end
diff --git a/src/closures/nondimensional_exchange_functions.jl b/src/closures/nondimensional_exchange_functions.jl
@@ -1,6 +1,6 @@
 """
     nondimensional_exchange_functions(
-        c_ε,
+        param_set,
         c_δ,
         c_μ,
         c_μ0,
@@ -9,29 +9,29 @@
         Δw,
         Δb,
         εδ_model,
-    ) 
+    )
 
 Returns the nondimensional entrainment and detrainment
 functions following Cohen et al. (JAMES, 2020), given:
- - `c_ε`,  entrainment factor  
- - `c_δ`,  detrainment factor   
- - `c_μ`,  logistic function factor  
- - `c_μ0`, logistic function timescale (sec)   
+ - `param_set` a parameter set
+ - `c_δ`,  detrainment factor
+ - `c_μ`,  logistic function factor
+ - `c_μ0`, logistic function timescale (sec)
  - `β`,    sorting power
  - `χ_up`, updraft mixing fraction
  - `Δw`, updraft - environment vertical velocity differnce
  - `Δb`, updraft - environment buoynacy differnce
  - `εδ_model`, entrainment detrainment model type
 """
-function nondimensional_exchange_functions(c_ε, c_δ, c_μ, c_μ0, β, χ_upd, Δw, Δb, εδ_model)
+function nondimensional_exchange_functions(param_set, c_δ, c_μ, c_μ0, β, χ_upd, Δw, Δb, εδ_model::MoistureDeficitEntr)
 
     # should be: c_δ = sign(condensate(ts_en) + condensate(ts_up[i])) * entr.c_δ
-    μ_ij = (χ_upd - εδ_model.a_upd / (εδ_model.a_upd + εδ_model.a_env)) * Δb / Δw
+    μ_ij = (χ_upd - εδ_model.a_up / (εδ_model.a_up + εδ_model.a_en)) * Δb / Δw
     D_ε = 1.0 / (1.0 + exp(-c_μ / c_μ0 * μ_ij))
     D_δ = 1.0 / (1.0 + exp(c_μ / c_μ0 * μ_ij))
 
-    M_δ = (max((εδ_model.RH_upd / 100.0)^β - (εδ_model.RH_env / 100.0)^β, 0.0))^(1.0 / β)
-    M_ε = (max((εδ_model.RH_env / 100.0)^β - (εδ_model.RH_upd / 100.0)^β, 0.0))^(1.0 / β)
+    M_δ = (max((εδ_model.RH_up / 100.0)^β - (εδ_model.RH_en / 100.0)^β, 0.0))^(1.0 / β)
+    M_ε = (max((εδ_model.RH_en / 100.0)^β - (εδ_model.RH_up / 100.0)^β, 0.0))^(1.0 / β)
 
     return D_ε, D_δ, M_δ, M_ε
 end;
diff --git a/src/types.jl b/src/types.jl
@@ -19,41 +19,59 @@ Base.@kwdef mutable struct mph_struct
     thl_rain_src::Float64 = 0
     qr_src::Float64 = 0
 end
+
+"""
+    EntrDetr
+
+$(DocStringExtensions.FIELDS)
+"""
+Base.@kwdef struct EntrDetr{FT}
+    "Dynamical entrainment"
+    ε_dyn::FT
+    "Dynamical detrainment"
+    δ_dyn::FT
+    "Turbulent entrainment"
+    ε_turb::FT
+    "Horizontal eddy-diffusivity"
+    K_ε::FT
+end
+
 """
- Entrainment detrainment model type
- - `εδ_model.b_upd`: updraft buoyancy
- - `εδ_model.b_env`: environment vertical velocity
- - `εδ_model.w_upd`: updraft vertical velocity
- - `εδ_model.w_env`: environment area fraction
- - `εδ_model.a_upd`: updraft area fraction
- - `εδ_model.a_env`: environment buoyancy
- - `εδ_model.ql_up`: updraft liquid water
- - `εδ_model.ql_env`: environment liquid water 
- - `εδ_model.RH_upd`: updraft relative humidity
- - `εδ_model.RH_env`: environment relative humidity
- - `εδ_model.M`: updraft momentum
- - `εδ_model.dMdz`: updraft momentum divergence
- - `εδ_model.tke`: env TKE
- - `εδ_model.N_up`: total number of updrafts
- - `εδ_model.ρ`: referance density
+    MoistureDeficitEntr
+
+My entrainment detrainment model
+
+$(DocStringExtensions.FIELDS)
 """
-Base.@kwdef mutable struct entr_detr_model
-    b_upd::Float64 = 0
-    b_env::Float64 = 0
-    w_upd::Float64 = 0
-    w_env::Float64 = 0
-    a_upd::Float64 = 0
-    a_env::Float64 = 0
-    ql_upd::Float64 = 0
-    ql_env::Float64 = 0
-    RH_upd::Float64 = 0
-    RH_env::Float64 = 0
-    M::Float64 = 0
-    dMdz::Float64 = 0
-    tke::Float64 = 0
-    n_up::Float64 = 0
-    ρ::Float64 = 0
-    R_up::Float64 = 0
+Base.@kwdef struct MoistureDeficitEntr{FT}
+    "updraft liquid water"
+    q_liq_up::FT
+    "environment liquid water"
+    q_liq_en::FT
+    "updraft vertical velocity"
+    w_up::FT
+    "environment vertical velocity"
+    w_en::FT
+    "updraft buoyancy"
+    b_up::FT
+    "environment buoyancy"
+    b_en::FT
+    "environment tke"
+    tke::FT
+    "updraft momentum divergence"
+    dMdz::FT
+    "updraft momentum"
+    M::FT
+    "updraft area fraction"
+    a_up::FT
+    "environment area fraction"
+    a_en::FT
+    "pressure plume spacing"
+    R_up::FT
+    "updraft relative humidity"
+    RH_up::FT
+    "environment relative humidity"
+    RH_en::FT
 end
 
 struct RainVariable{T}

diff --git a/src/utility_functions.jl b/src/utility_functions.jl
@@ -22,7 +22,6 @@ function logistic(x, slope, mid)
 end
 
 function lamb_smooth_minimum(l, lower_bound, upper_bound)
-    leng = size(l)
     x_min = minimum(l)
     λ_0 = max(x_min * lower_bound / real(LambertW.lambertw(2.0 / MathConstants.e)), upper_bound)