Add harmonics limiter for entr/detr, set GABLS config

.buildkite/pipeline.yml

@@ -670,7 +670,6 @@ steps:
           julia --color=yes --project=examples examples/hybrid/driver.jl
           --config_file $CONFIG_PATH/edmfx_gabls_box.yml
         artifact_paths: "edmfx_gabls_box/*"
-        soft_fail: true
         agents:
           slurm_mem: 20GB
 

config/model_configs/edmfx_gabls_box.yml

@@ -20,9 +20,10 @@ x_elem: 2
 y_elem: 2
 z_elem: 8
 z_stretch: false
-dt: "10secs"
+dt: "5secs"
 t_end: "9hours"
 dt_save_to_disk: "10mins"
 perturb_initstate: false
 FLOAT_TYPE: "Float64"
-toml: [toml/edmfx_box_gabls.toml]
+toml: [toml/edmfx_box.toml]
+#toml: [toml/edmfx_box_gabls.toml]

examples/hybrid/driver.jl

@@ -29,9 +29,9 @@ import OrderedCollections
 using ClimaCoreTempestRemap
 using ClimaCorePlots, Plots
 using ClimaCoreMakie, CairoMakie
-if config.comms_ctx isa ClimaComms.SingletonCommsContext
-    include(joinpath(pkgdir(CA), "perf", "jet_report_nfailures.jl"))
-end
+#if config.comms_ctx isa ClimaComms.SingletonCommsContext
+#    include(joinpath(pkgdir(CA), "perf", "jet_report_nfailures.jl"))
+#end
 include(joinpath(pkgdir(CA), "post_processing", "contours_and_profiles.jl"))
 include(joinpath(pkgdir(CA), "post_processing", "post_processing_funcs.jl"))
 include(

src/prognostic_equations/edmfx_entr_detr.jl

@@ -2,6 +2,9 @@
 ##### EDMF entrainment detrainment
 #####
 
+# return a harmonic mean of (a, 1-a)
+hm_limiter(a) = 2 * a * (1 - a)
+
 """
    Return entrainment rate [1/s].
 
@@ -107,6 +110,28 @@ function entrainment(
     return entr
 end
 
+function entrainment(
+    params,
+    ᶜz::FT,
+    z_sfc::FT,
+    ᶜp::FT,
+    ᶜρ::FT,
+    buoy_flux_surface::FT,
+    ᶜaʲ::FT,
+    ᶜwʲ::FT,
+    ᶜRHʲ::FT,
+    ᶜbuoyʲ::FT,
+    ᶜw⁰::FT,
+    ᶜRH⁰::FT,
+    ᶜbuoy⁰::FT,
+    dt::FT,
+    ::ConstantCoefficientHarmonicsEntrainment,
+) where {FT}
+    entr_coeff = CAP.entr_coeff(params)
+    entr = min(entr_coeff * abs(ᶜwʲ) / (ᶜz - z_sfc), 1 / dt)
+    return entr * FT(2) * hm_limiter(ᶜaʲ)
+end
+
 function entrainment(
     params,
     ᶜz::FT,
@@ -236,12 +261,34 @@ function detrainment(
     return detr
 end
 
+function detrainment(
+    params,
+    ᶜz::FT,
+    z_sfc::FT,
+    ᶜp::FT,
+    ᶜρ::FT,
+    buoy_flux_surface::FT,
+    ᶜaʲ::FT,
+    ᶜwʲ::FT,
+    ᶜRHʲ::FT,
+    ᶜbuoyʲ::FT,
+    ᶜw⁰::FT,
+    ᶜRH⁰::FT,
+    ᶜbuoy⁰::FT,
+    dt::FT,
+    ::ConstantCoefficientHarmonicsDetrainment,
+) where {FT}
+    detr_coeff = CAP.detr_coeff(params)
+    detr = min(detr_coeff * abs(ᶜwʲ), 1 / dt)
+    return detr * FT(2) * hm_limiter(ᶜaʲ)
+end
+
 edmfx_entr_detr_tendency!(Yₜ, Y, p, t, colidx, turbconv_model) = nothing
 function edmfx_entr_detr_tendency!(Yₜ, Y, p, t, colidx, turbconv_model::EDMFX)
 
     n = n_mass_flux_subdomains(turbconv_model)
     (; ᶜspecificʲs, ᶜh_totʲs, ᶜentrʲs, ᶜdetrʲs) = p
-    (; ᶜu⁰, ᶜspecific⁰, ᶜh_tot⁰) = p
+    (; ᶜspecific⁰, ᶜh_tot⁰, ᶠu₃⁰) = p
 
     for j in 1:n
 
@@ -262,8 +309,8 @@ function edmfx_entr_detr_tendency!(Yₜ, Y, p, t, colidx, turbconv_model::EDMFX)
             )
 
         @. Yₜ.f.sgsʲs.:($$j).u₃[colidx] +=
-            ᶠinterp(ᶜentrʲs.:($$j)[colidx] * C3(ᶜu⁰[colidx])) -
-            ᶠinterp(ᶜentrʲs.:($$j)[colidx]) * Y.f.sgsʲs.:($$j).u₃[colidx]
+            ᶠinterp(ᶜentrʲs.:($$j)[colidx]) *
+            (ᶠu₃⁰[colidx] - Y.f.sgsʲs.:($$j).u₃[colidx])
     end
     return nothing
 end

src/solver/model_getters.jl

@@ -372,6 +372,8 @@ function get_entrainment_model(parsed_args)
         PiGroupsEntrainment()
     elseif entr_model == "ConstantCoefficient"
         ConstantCoefficientEntrainment()
+    elseif entr_model == "ConstantCoefficientHarmonics"
+        ConstantCoefficientHarmonicsEntrainment()
     elseif entr_model == "ConstantTimescale"
         ConstantTimescaleEntrainment()
     else
@@ -387,6 +389,8 @@ function get_detrainment_model(parsed_args)
         PiGroupsDetrainment()
     elseif detr_model == "ConstantCoefficient"
         ConstantCoefficientDetrainment()
+    elseif detr_model == "ConstantCoefficientHarmonics"
+        ConstantCoefficientHarmonicsDetrainment()
     else
         error("Invalid entr_model $(entr_model)")
     end

src/solver/types.jl

@@ -179,13 +179,15 @@ abstract type AbstractEntrainmentModel end
 struct NoEntrainment <: AbstractEntrainmentModel end
 struct PiGroupsEntrainment <: AbstractEntrainmentModel end
 struct ConstantCoefficientEntrainment <: AbstractEntrainmentModel end
+struct ConstantCoefficientHarmonicsEntrainment <: AbstractEntrainmentModel end
 struct ConstantTimescaleEntrainment <: AbstractEntrainmentModel end
 
 abstract type AbstractDetrainmentModel end
 
 struct NoDetrainment <: AbstractDetrainmentModel end
 struct PiGroupsDetrainment <: AbstractDetrainmentModel end
 struct ConstantCoefficientDetrainment <: AbstractDetrainmentModel end
+struct ConstantCoefficientHarmonicsDetrainment <: AbstractDetrainmentModel end
 
 abstract type AbstractQuadratureType end
 struct LogNormalQuad <: AbstractQuadratureType end

toml/edmfx_box.toml

@@ -10,6 +10,6 @@ type = "float"
 
 [EDMF_min_area]
 alias = "min_area"
-value = 1.0e-2
+value = 1.0e-3
 type = "float"
 description = "Minimum area fraction per updraft. Parameter not described in the literature."

toml/edmfx_box_gabls.toml

@@ -3,7 +3,25 @@ alias = "C_E"
 value = 0.044
 type = "float"
 
+[entr_coeff]
+alias = "entr_coeff"
+value = 1.0
+type = "float"
+description = "TODO: Remove this. Constant entrainment coefficient used for testing EDMF"
+
+[detr_coeff]
+alias = "detr_coeff"
+value = 0.001
+type = "float"
+description = "TODO: Remove this. Constant entrainment coefficient used for testing EDMF"
+
 [EDMF_surface_area]
 alias = "surface_area"
-value = 1e-4
+value = 1e-1
+type = "float"
+
+[EDMF_min_area]
+alias = "min_area"
+value = 1.0e-3
 type = "float"
+description = "Minimum area fraction per updraft. Parameter not described in the literature."