Add min area limiter

CliMA · Dec 1, 2023 · 9b33c30 · 9b33c30
1 parent dab93c7
commit 9b33c30
Show file tree

Hide file tree

Showing 5 changed files with 35 additions and 5 deletions.
diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml
@@ -1,6 +1,6 @@
 env:
   JULIA_VERSION: "1.9.3"
-  OPENMPI_VERSION: "4.0.4"
+  OPENMPI_VERSION: "4.1.5"
   CUDA_VERSION: "11.2"
   CLIMACOMMS_DEVICE: "CPU"
   JULIA_CPU_TARGET: 'broadwell;skylake'

diff --git a/driver/generate_namelist.jl b/driver/generate_namelist.jl
@@ -91,7 +91,7 @@ function default_namelist(
     namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["surface_area_bc"] = "Fixed" #{"Fixed", "Prognostic", "Closure"}
     namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["surface_area_bc_params"] = [1e-3, 1e-3, 1e-3]
     namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["max_area"] = 0.9
-    namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["min_area"] = 1e-5
+    namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["min_area"] = 1e-10
 
     # mixing_length
     namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["tke_ed_coeff"] = 0.14
@@ -113,6 +113,9 @@ function default_namelist(
     namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["updraft_mixing_frac"] = 0.25
     namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["area_limiter_scale"] = 4.0
     namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["area_limiter_power"] = 10.0
+    namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["limit_min_area"] = false
+    namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["min_area_limiter_scale"] = 10.0
+    namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["min_area_limiter_power"] = 2000.0
     namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["entrainment_scale"] = 0.0004
     namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["sorting_power"] = 2.0
     namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["min_upd_velocity"] = 0.001
@@ -159,7 +162,7 @@ function default_namelist(
     namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["entr_dim_scale"] = "buoy_vel" # {"buoy_vel", "inv_scale_height", "inv_z", "none"}
     namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["detr_dim_scale"] = "buoy_vel" # {"buoy_vel", "inv_scale_height", "inv_z", "none"}
     namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["entr_pi_subset"] = ntuple(i -> i, 6) # or, e.g., (1, 3, 6)
-    namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["pi_norm_consts"] = [478.298, 1.0, 1.0, 1.0, 1.0, 1.0] # normalization constants for Pi groups
+    namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["pi_norm_consts"] = [10.0, 5.0, 1.0, 1.0, 1.0, 1.0] # normalization constants for Pi groups
     namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["stochastic_entrainment"] = "deterministic"  # {"deterministic", "noisy_relaxation_process", "lognormal_scaling", "prognostic_noisy_relaxation_process"}
 
     namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["pressure_closure_buoy"] = "normalmode"

diff --git a/src/closures/entr_detr.jl b/src/closures/entr_detr.jl
@@ -171,10 +171,20 @@ function εδ_dyn(εδ_model, εδ_vars, entr_dim_scale, detr_dim_scale, ε_nond
     δ_dim_scale = entrainment_inv_length_scale(εδ_model, εδ_vars, detr_dim_scale)
 
     area_limiter = max_area_limiter(εδ_model, εδ_vars.max_area, εδ_vars.a_up)
+    min_limiter = min_area_limiter(εδ_model, εδ_vars.a_up)
+
+    ε_dim_scale = min(ε_dim_scale, 1.0)
+    δ_dim_scale = min(δ_dim_scale, 1.0)
 
     # fractional dynamical entrainment / detrainment [1 / m]
     ε_dyn = ε_dim_scale * ε_nondim
-    δ_dyn = δ_dim_scale * (δ_nondim + area_limiter)
+    δ_dyn = δ_dim_scale * δ_nondim
+
+    if εδ_params(εδ_model).limit_min_area
+        ε_dyn += max(δ_dyn, ε_dyn) * min_limiter
+    end
+
+    δ_dyn += max(ε_dyn, δ_dyn) * area_limiter
 
     return ε_dyn, δ_dyn
 end

diff --git a/src/closures/nondimensional_exchange_functions.jl b/src/closures/nondimensional_exchange_functions.jl
@@ -7,6 +7,13 @@ function max_area_limiter(εδ_model, max_area, a_up)
     return A * exp(-k * (max_area - a_up))
 end
 
+function min_area_limiter(εδ_model, a_up)
+    FT = eltype(a_up)
+    A = εδ_params(εδ_model).min_area_limiter_scale
+    k = εδ_params(εδ_model).min_area_limiter_power
+    return A * exp(-k * a_up)
+end
+
 function non_dimensional_groups(εδ_model, εδ_model_vars)
     FT = eltype(εδ_model_vars.tke_en)
     Δw = get_Δw(εδ_model, εδ_model_vars.w_up, εδ_model_vars.w_en)
@@ -19,6 +26,10 @@ function non_dimensional_groups(εδ_model, εδ_model_vars)
     Π₅ = εδ_model_vars.zc_i / εδ_model_vars.H_up / Π_norm[5]
     Π₆ = εδ_model_vars.zc_i / εδ_model_vars.ref_H / Π_norm[6]
 
+    # Π₁, Π₂ are unbounded, so clip values that blow up
+    Π₁ = clamp(Π₁, -1.0, 1.0)
+    Π₂ = clamp(Π₂, -1.0, 1.0)
+
     Π_groups = (Π₁, Π₂, Π₃, Π₄, Π₅, Π₆)
 
     return map(i -> Π_groups[i], εδ_model_vars.entr_Π_subset)

diff --git a/src/types.jl b/src/types.jl
@@ -51,6 +51,9 @@ Base.@kwdef struct εδModelParams{FT, AFT}
     c_λ::FT # scaling factor for TKE in entrainment scale calculations
     γ_lim::FT
     β_lim::FT
+    limit_min_area::Bool
+    min_area_limiter_scale::FT
+    min_area_limiter_power::FT
     c_γ::FT # scaling factor for turbulent entrainment rate
     c_δ::FT # factor multiplier for moist term in entrainment/detrainment
     Π_norm::AFT
@@ -617,12 +620,15 @@ function EDMFModel(::Type{FT}, namelist, precip_model, rain_formation_model) whe
     c_λ = parse_namelist(namelist, "turbulence", "EDMF_PrognosticTKE", "entrainment_smin_tke_coeff")
     γ_lim = parse_namelist(namelist, "turbulence", "EDMF_PrognosticTKE", "area_limiter_scale")
     β_lim = parse_namelist(namelist, "turbulence", "EDMF_PrognosticTKE", "area_limiter_power")
+    limit_min_area = parse_namelist(namelist, "turbulence", "EDMF_PrognosticTKE", "limit_min_area")
+    min_area_limiter_scale = parse_namelist(namelist, "turbulence", "EDMF_PrognosticTKE", "min_area_limiter_scale")
+    min_area_limiter_power = parse_namelist(namelist, "turbulence", "EDMF_PrognosticTKE", "min_area_limiter_power")
     c_γ = parse_namelist(namelist, "turbulence", "EDMF_PrognosticTKE", "turbulent_entrainment_factor")
     c_δ = parse_namelist(namelist, "turbulence", "EDMF_PrognosticTKE", "detrainment_factor")
     Π_norm = parse_namelist(namelist, "turbulence", "EDMF_PrognosticTKE", "pi_norm_consts")
     Π_norm = SA.SVector{length(Π_norm), FT}(Π_norm)
 
-    εδ_params = εδModelParams{FT, typeof(Π_norm)}(; w_min, c_ε, μ_0, β, χ, c_λ, γ_lim, β_lim, c_γ, c_δ, Π_norm)
+    εδ_params = εδModelParams{FT, typeof(Π_norm)}(; w_min, c_ε, μ_0, β, χ, c_λ, γ_lim, β_lim, limit_min_area, min_area_limiter_scale, min_area_limiter_power, c_γ, c_δ, Π_norm)
 
     entr_pi_subset = Tuple(parse_namelist(namelist, "turbulence", "EDMF_PrognosticTKE", "entr_pi_subset"))