Merge pull request #2882 from CliMA/zs/pigroup

update pi groups in entrainment and detrainment
CliMA · Apr 8, 2024 · 77b10a6 · 77b10a6
2 parents d8dd580 + 7e2f82b
commit 77b10a6
Show file tree

Hide file tree

Showing 7 changed files with 146 additions and 47 deletions.
diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml
@@ -558,6 +558,14 @@ steps:
         agents:
           slurm_mem: 20GB
 
+      - label: ":genie: Prognostic EDMFX Bomex in a column (explicit)"
+        command: >
+          julia --color=yes --project=examples examples/hybrid/driver.jl
+          --config_file $CONFIG_PATH/prognostic_edmfx_bomex_explicit_column.yml
+        artifact_paths: "prognostic_edmfx_bomex_explicit_column/output_active/*"
+        agents:
+          slurm_mem: 20GB
+
       - label: ":genie: Prognostic EDMFX Bomex with precribed TKE in a column"
         command: >
           julia --color=yes --project=examples examples/hybrid/driver.jl

diff --git a/config/model_configs/prognostic_edmfx_bomex_explicit_column.yml b/config/model_configs/prognostic_edmfx_bomex_explicit_column.yml
@@ -0,0 +1,39 @@
+job_id: "prognostic_edmfx_bomex_explicit_column"
+initial_condition: "Bomex"
+subsidence: "Bomex"
+edmf_coriolis: "Bomex"
+ls_adv: "Bomex"
+surface_setup: "Bomex"
+turbconv: "prognostic_edmfx"
+implicit_diffusion: false
+implicit_sgs_advection: false
+approximate_linear_solve_iters: 2
+max_newton_iters_ode: 3
+edmfx_upwinding: first_order
+edmfx_entr_model: "PiGroups"
+edmfx_detr_model: "PiGroups"
+edmfx_sgs_mass_flux: true
+edmfx_sgs_diffusive_flux: true
+edmfx_nh_pressure: true
+edmfx_velocity_relaxation: true
+prognostic_tke: true
+cloud_model: "grid_scale"
+moist: "equil"
+config: "column"
+z_max: 3e3
+z_elem: 60
+z_stretch: false
+perturb_initstate: false
+dt: "1secs"
+t_end: "6hours"
+dt_save_to_disk: "10mins"
+toml: [toml/prognostic_edmfx_bomex_explicit.toml]
+netcdf_output_at_levels: true
+netcdf_interpolation_num_points: [2, 2, 60]
+diagnostics:
+  - short_name: [ts, ta, thetaa, ha, pfull, rhoa, ua, va, wa, hur, hus, cl, clw, cli, hussfc, evspsbl]
+    period: 10mins
+  - short_name: [arup, waup, taup, thetaaup, haup, husup, hurup, clwup, cliup, waen, taen, thetaaen, haen, husen, huren, clwen, clien, tke]
+    period: 10mins
+  - short_name: [entr, detr, lmix, bgrad, strain, edt, evu]
+    period: 10mins
diff --git a/post_processing/ci_plots.jl b/post_processing/ci_plots.jl
@@ -945,6 +945,7 @@ EDMFBoxPlots = Union{
     Val{:prognostic_edmfx_bomex_fixtke_column},
     Val{:prognostic_edmfx_bomex_column},
     Val{:prognostic_edmfx_bomex_stretched_column},
+    Val{:prognostic_edmfx_bomex_explicit_column},
     Val{:prognostic_edmfx_dycoms_rf01_column},
     Val{:prognostic_edmfx_trmm_column_0M},
     Val{:prognostic_edmfx_simpleplume_column},

diff --git a/src/cache/diagnostic_edmf_precomputed_quantities.jl b/src/cache/diagnostic_edmf_precomputed_quantities.jl
@@ -228,7 +228,6 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_do_integral!(
     (; ᶜΦ) = p.core
     (; ᶜp, ᶠu³, ᶜts, ᶜh_tot, ᶜK) = p.precomputed
     (; q_tot) = p.precomputed.ᶜspecific
-    (; buoyancy_flux) = p.precomputed.sfc_conditions
     (;
         ᶜρaʲs,
         ᶠu³ʲs,
@@ -241,7 +240,7 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_do_integral!(
         ᶜdetrʲs,
         ᶠnh_pressure³ʲs,
     ) = p.precomputed
-    (; ᶠu³⁰, ᶜK⁰) = p.precomputed
+    (; ᶠu³⁰, ᶜK⁰, ᶜtke⁰) = p.precomputed
 
     thermo_params = CAP.thermodynamics_params(params)
     microphys_params = CAP.microphysics_params(params)
@@ -256,7 +255,6 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_do_integral!(
 
     z_sfc_halflevel =
         Fields.field_values(Fields.level(Fields.coordinate_field(Y.f).z, half))
-    buoyancy_flux_sfc_halflevel = Fields.field_values(buoyancy_flux)
 
     # integral
     for i in 2:Spaces.nlevels(axes(Y.c))
@@ -342,13 +340,13 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_do_integral!(
                 CAP.R_d(params) * CAP.T_surf_ref(params) / CAP.grav(params)
             S_q_totʲ_prev_level =
                 Fields.field_values(Fields.level(ᶜS_q_totʲ, i - 1))
+            tke_prev_level = Fields.field_values(Fields.level(ᶜtke⁰, i - 1))
             @. entrʲ_prev_level = entrainment(
                 params,
                 z_prev_level,
                 z_sfc_halflevel,
                 p_prev_level,
                 ρ_prev_level,
-                buoyancy_flux_sfc_halflevel,
                 draft_area(ρaʲ_prev_level, ρʲ_prev_level),
                 get_physical_w(
                     u³ʲ_prev_halflevel,
@@ -362,6 +360,7 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_do_integral!(
                 ),
                 TD.relative_humidity(thermo_params, ts_prev_level),
                 FT(0),
+                tke_prev_level,
                 p.atmos.edmfx_entr_model,
             )
 
@@ -490,7 +489,6 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_do_integral!(
                 z_sfc_halflevel,
                 p_prev_level,
                 ρ_prev_level,
-                buoyancy_flux_sfc_halflevel,
                 ρaʲ_prev_level,
                 draft_area(ρaʲ_prev_level, ρʲ_prev_level),
                 get_physical_w(
@@ -508,6 +506,7 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_do_integral!(
                 entrʲ_prev_level,
                 vert_div_level,
                 FT(0), # mass flux divergence is not implemented for diagnostic edmf
+                tke_prev_level,
                 p.atmos.edmfx_detr_model,
             )
 

diff --git a/src/cache/prognostic_edmf_precomputed_quantities.jl b/src/cache/prognostic_edmf_precomputed_quantities.jl
@@ -188,7 +188,7 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_closures!(
     FT = eltype(params)
     n = n_mass_flux_subdomains(turbconv_model)
 
-    (; ᶜtke⁰, ᶜu, ᶜp, ᶜρa⁰, ᶠu³⁰, ᶜts⁰, ᶜρ⁰, ᶜq_tot⁰) = p.precomputed
+    (; ᶜtke⁰, ᶜu, ᶜp, ᶜρa⁰, ᶠu³⁰, ᶜts⁰, ᶜq_tot⁰) = p.precomputed
     (;
         ᶜmixing_length,
         ᶜlinear_buoygrad,
@@ -198,7 +198,7 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_closures!(
         ρatke_flux,
     ) = p.precomputed
     (; ᶜuʲs, ᶜtsʲs, ᶠu³ʲs, ᶜρʲs, ᶜentrʲs, ᶜdetrʲs) = p.precomputed
-    (; ustar, obukhov_length, buoyancy_flux) = p.precomputed.sfc_conditions
+    (; ustar, obukhov_length) = p.precomputed.sfc_conditions
 
     ᶜz = Fields.coordinate_field(Y.c).z
     z_sfc = Fields.level(Fields.coordinate_field(Y.f).z, Fields.half)
@@ -215,14 +215,14 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_closures!(
             z_sfc,
             ᶜp,
             Y.c.ρ,
-            buoyancy_flux,
             draft_area(Y.c.sgsʲs.:($$j).ρa, ᶜρʲs.:($$j)),
             get_physical_w(ᶜuʲs.:($$j), ᶜlg),
             TD.relative_humidity(thermo_params, ᶜtsʲs.:($$j)),
             ᶜphysical_buoyancy(params, Y.c.ρ, ᶜρʲs.:($$j)),
             get_physical_w(ᶜu, ᶜlg),
             TD.relative_humidity(thermo_params, ᶜts⁰),
             FT(0),
+            max(ᶜtke⁰, 0),
             p.atmos.edmfx_entr_model,
         )
         @. ᶜentrʲs.:($$j) = limit_entrainment(
@@ -239,7 +239,6 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_closures!(
             z_sfc,
             ᶜp,
             Y.c.ρ,
-            buoyancy_flux,
             Y.c.sgsʲs.:($$j).ρa,
             draft_area(Y.c.sgsʲs.:($$j).ρa, ᶜρʲs.:($$j)),
             get_physical_w(ᶜuʲs.:($$j), ᶜlg),
@@ -251,6 +250,7 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_closures!(
             ᶜentrʲs.:($$j),
             ᶜvert_div,
             ᶜmassflux_vert_div,
+            ᶜtke⁰,
             p.atmos.edmfx_detr_model,
         )
         @. ᶜdetrʲs.:($$j) = limit_detrainment(

diff --git a/src/prognostic_equations/edmfx_entr_detr.jl b/src/prognostic_equations/edmfx_entr_detr.jl
@@ -2,9 +2,6 @@
 ##### EDMF entrainment detrainment
 #####
 
-# return a harmonic mean of (a, 1-a)
-hm_limiter(a) = 2 * a * (1 - a)
-
 """
    Return entrainment rate [1/s].
 
@@ -25,14 +22,14 @@ function entrainment(
     z_sfc::FT,
     ᶜp::FT,
     ᶜρ::FT,
-    buoy_flux_surface::FT,
     ᶜaʲ::FT,
     ᶜwʲ::FT,
     ᶜRHʲ::FT,
     ᶜbuoyʲ::FT,
     ᶜw⁰::FT,
     ᶜRH⁰::FT,
     ᶜbuoy⁰::FT,
+    ᶜtke⁰::FT,
     ::NoEntrainment,
 ) where {FT}
     return FT(0)
@@ -44,38 +41,36 @@ function entrainment(
     z_sfc::FT,
     ᶜp::FT,
     ᶜρ::FT,
-    buoy_flux_surface::FT,
     ᶜaʲ::FT,
     ᶜwʲ::FT,
     ᶜRHʲ::FT,
     ᶜbuoyʲ::FT,
     ᶜw⁰::FT,
     ᶜRH⁰::FT,
     ᶜbuoy⁰::FT,
+    ᶜtke⁰::FT,
     ::PiGroupsEntrainment,
 ) where {FT}
     if ᶜaʲ <= FT(0)
         return FT(0)
     else
         g = CAP.grav(params)
-
-        # pressure scale height (height where pressure drops by 1/e)
         ref_H = ᶜp / (ᶜρ * g)
-        # convective velocity
-        w_star = get_wstar(buoy_flux_surface)
 
         # non-dimensional pi-groups
-        # TODO - using Π₁ blows things up
-        Π₁ =
-            (ᶜz - z_sfc) * (ᶜbuoyʲ - ᶜbuoy⁰) /
-            ((ᶜwʲ - ᶜw⁰)^2 + w_star^2 + eps(FT))
+        Π₁ = (ᶜz - z_sfc) * (ᶜbuoyʲ - ᶜbuoy⁰) / ((ᶜwʲ - ᶜw⁰)^2 + eps(FT)) / 100
+        Π₂ = ᶜtke⁰ / ((ᶜwʲ - ᶜw⁰)^2 + eps(FT)) / 2
         Π₃ = sqrt(ᶜaʲ)
         Π₄ = ᶜRHʲ - ᶜRH⁰
-        Π₆ = (ᶜz - z_sfc) / ref_H
+        Π₅ = (ᶜz - z_sfc) / ref_H
+        # Π₁, Π₂ are unbounded, so clip values that blow up
+        Π₁ = min(max(Π₁, -1), 1)
+        Π₂ = min(max(Π₂, -1), 1)
         entr =
-            abs(ᶜwʲ) / (ᶜz - z_sfc) * (
-                -4.013288 - 0.000968 * Π₁ + 0.356974 * Π₃ - 0.403124 * Π₄ +
-                1.503261 * Π₆
+            abs(ᶜwʲ - ᶜw⁰) / (ᶜz - z_sfc) * (
+                -0.32332 + 4.79372 * Π₁ + 3.98108 * Π₂ - 21.64173 * Π₃ +
+                18.395 * Π₄ +
+                1.12799 * Π₅
             )
 
         return entr
@@ -88,14 +83,14 @@ function entrainment(
     z_sfc::FT,
     ᶜp::FT,
     ᶜρ::FT,
-    buoy_flux_surface::FT,
     ᶜaʲ::FT,
     ᶜwʲ::FT,
     ᶜRHʲ::FT,
     ᶜbuoyʲ::FT,
     ᶜw⁰::FT,
     ᶜRH⁰::FT,
     ᶜbuoy⁰::FT,
+    ᶜtke⁰::FT,
     ::GeneralizedEntrainment,
 ) where {FT}
     turbconv_params = CAP.turbconv_params(params)
@@ -109,7 +104,9 @@ function entrainment(
         min_area_limiter_scale *
         exp(-min_area_limiter_power * (max(ᶜaʲ, 0) - a_min))
     entr =
-        entr_inv_tau + entr_coeff * abs(ᶜwʲ) / (ᶜz - z_sfc) + min_area_limiter
+        entr_inv_tau +
+        entr_coeff * abs(ᶜwʲ - ᶜw⁰) / (ᶜz - z_sfc) +
+        min_area_limiter
 
     return entr
 end
@@ -134,7 +131,6 @@ function detrainment(
     z_sfc::FT,
     ᶜp::FT,
     ᶜρ::FT,
-    buoy_flux_surface::FT,
     ᶜρaʲ::FT,
     ᶜaʲ::FT,
     ᶜwʲ::FT,
@@ -146,6 +142,7 @@ function detrainment(
     ᶜentr::FT,
     ᶜvert_div::FT,
     ᶜmassflux_vert_div::FT,
+    ᶜtke⁰::FT,
     ::NoDetrainment,
 ) where {FT}
     return FT(0)
@@ -157,7 +154,6 @@ function detrainment(
     z_sfc::FT,
     ᶜp::FT,
     ᶜρ::FT,
-    buoy_flux_surface::FT,
     ᶜρaʲ::FT,
     ᶜaʲ::FT,
     ᶜwʲ::FT,
@@ -169,34 +165,29 @@ function detrainment(
     ᶜentr::FT,
     ᶜvert_div::FT,
     ᶜmassflux_vert_div::FT,
+    ᶜtke⁰::FT,
     ::PiGroupsDetrainment,
 ) where {FT}
 
     if ᶜaʲ <= FT(0)
         return FT(0)
     else
         g = CAP.grav(params)
-        turbconv_params = CAP.turbconv_params(params)
-        ᶜaʲ_max = CAP.max_area(turbconv_params)
-
-        max_area_limiter = FT(0.1) * exp(-10 * (ᶜaʲ_max - ᶜaʲ))
-        # pressure scale height (height where pressure drops by 1/e)
         ref_H = ᶜp / (ᶜρ * g)
-        # convective velocity
-        w_star = get_wstar(buoy_flux_surface)
 
         # non-dimensional pi-groups
-        # TODO - using Π₁ blows things up
-        Π₁ =
-            (ᶜz - z_sfc) * (ᶜbuoyʲ - ᶜbuoy⁰) /
-            ((ᶜwʲ - ᶜw⁰)^2 + w_star^2 + eps(FT))
+        Π₁ = (ᶜz - z_sfc) * (ᶜbuoyʲ - ᶜbuoy⁰) / ((ᶜwʲ - ᶜw⁰)^2 + eps(FT)) / 100
+        Π₂ = ᶜtke⁰ / ((ᶜwʲ - ᶜw⁰)^2 + eps(FT)) / 2
         Π₃ = sqrt(ᶜaʲ)
         Π₄ = ᶜRHʲ - ᶜRH⁰
-        Π₆ = (ᶜz - z_sfc) / ref_H
+        Π₅ = (ᶜz - z_sfc) / ref_H
+        # Π₁, Π₂ are unbounded, so clip values that blow up
+        Π₁ = min(max(Π₁, -1), 1)
+        Π₂ = min(max(Π₂, -1), 1)
         detr =
-            abs(ᶜwʲ) * (
-                3.535208 + 0.598496 * Π₁ + 1.583348 * Π₃ + 0.046275 * Π₄ -
-                0.344836 * Π₆ + max_area_limiter
+            -min(ᶜmassflux_vert_div, 0) / ᶜρaʲ * (
+                0.3410 - 0.56153 * Π₁ - 0.53411 * Π₂ + 6.01925 * Π₃ -
+                1.47516 * Π₄ - 3.85788 * Π₅
             )
         return detr
     end
@@ -208,7 +199,6 @@ function detrainment(
     z_sfc::FT,
     ᶜp::FT,
     ᶜρ::FT,
-    buoy_flux_surface::FT,
     ᶜρaʲ::FT,
     ᶜaʲ::FT,
     ᶜwʲ::FT,
@@ -220,6 +210,7 @@ function detrainment(
     ᶜentr::FT,
     ᶜvert_div::FT,
     ᶜmassflux_vert_div::FT,
+    ᶜtke⁰::FT,
     ::GeneralizedDetrainment,
 ) where {FT}
     turbconv_params = CAP.turbconv_params(params)
@@ -258,7 +249,6 @@ function detrainment(
     z_sfc::FT,
     ᶜp::FT,
     ᶜρ::FT,
-    buoy_flux_surface::FT,
     ᶜρaʲ::FT,
     ᶜaʲ::FT,
     ᶜwʲ::FT,