From 7e2f82babc5e649f4c00e4aa5f2cb27d8cbaa760 Mon Sep 17 00:00:00 2001 From: Zhaoyi Shen <11598433+szy21@users.noreply.github.com> Date: Fri, 5 Apr 2024 23:04:04 -0700 Subject: [PATCH] update pi groups --- .buildkite/pipeline.yml | 8 +++ ...prognostic_edmfx_bomex_explicit_column.yml | 39 +++++++++++ post_processing/ci_plots.jl | 1 + .../diagnostic_edmf_precomputed_quantities.jl | 9 ++- .../prognostic_edmf_precomputed_quantities.jl | 8 +-- src/prognostic_equations/edmfx_entr_detr.jl | 66 ++++++++----------- toml/prognostic_edmfx_bomex_explicit.toml | 62 +++++++++++++++++ 7 files changed, 146 insertions(+), 47 deletions(-) create mode 100644 config/model_configs/prognostic_edmfx_bomex_explicit_column.yml create mode 100644 toml/prognostic_edmfx_bomex_explicit.toml diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml index fd7e92ef51..cc1954c11c 100644 --- 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 new file mode 100644 index 0000000000..17a93b9003 --- /dev/null +++ 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 index f072752667..a661c0f66f 100644 --- 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 index 860207a493..c5f3971ccd 100644 --- 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 index 6abc5a016a..8542205849 100644 --- 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,7 +215,6 @@ 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)), @@ -223,6 +222,7 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_closures!( 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 index 2cf7ccb3c8..9a42720944 100644 --- 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,7 +22,6 @@ function entrainment( z_sfc::FT, ᶜp::FT, ᶜρ::FT, - buoy_flux_surface::FT, ᶜaʲ::FT, ᶜwʲ::FT, ᶜRHʲ::FT, @@ -33,6 +29,7 @@ function entrainment( ᶜw⁰::FT, ᶜRH⁰::FT, ᶜbuoy⁰::FT, + ᶜtke⁰::FT, ::NoEntrainment, ) where {FT} return FT(0) @@ -44,7 +41,6 @@ function entrainment( z_sfc::FT, ᶜp::FT, ᶜρ::FT, - buoy_flux_surface::FT, ᶜaʲ::FT, ᶜwʲ::FT, ᶜRHʲ::FT, @@ -52,30 +48,29 @@ function entrainment( ᶜ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,7 +83,6 @@ function entrainment( z_sfc::FT, ᶜp::FT, ᶜρ::FT, - buoy_flux_surface::FT, ᶜaʲ::FT, ᶜwʲ::FT, ᶜRHʲ::FT, @@ -96,6 +90,7 @@ function entrainment( ᶜ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,6 +165,7 @@ function detrainment( ᶜentr::FT, ᶜvert_div::FT, ᶜmassflux_vert_div::FT, + ᶜtke⁰::FT, ::PiGroupsDetrainment, ) where {FT} @@ -176,27 +173,21 @@ function detrainment( 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, diff --git a/toml/prognostic_edmfx_bomex_explicit.toml b/toml/prognostic_edmfx_bomex_explicit.toml new file mode 100644 index 0000000000..6c256730fb --- /dev/null +++ b/toml/prognostic_edmfx_bomex_explicit.toml @@ -0,0 +1,62 @@ +[C_E] +value = 0.0176 + +[EDMF_surface_area] +value = 0.1 + +[EDMF_min_area] +value = 1.0e-5 + +[EDMF_max_area] +value = 0.5 + +[entr_inv_tau] +value = 0.0002 + +[entr_coeff] +value = 0.3 + +[min_area_limiter_scale] +value = 0 + +[min_area_limiter_power] +value = 10 + +[detr_inv_tau] +value = 0 + +[detr_coeff] +value = 0 + +[detr_buoy_coeff] +value = 0 + +[detr_vertdiv_coeff] +value = 0 + +[detr_massflux_vertdiv_coeff] +value = 1 + +[max_area_limiter_scale] +value = 0 + +[pressure_normalmode_drag_coeff] +value = 20.0 + +[pressure_normalmode_buoy_coeff1] +value = 0.857 + +[mixing_length_eddy_viscosity_coefficient] +value = 0.041 + +[mixing_length_diss_coeff] +value = 0.174 + +[mixing_length_static_stab_coeff] +value = 0.292 + +[mixing_length_tke_surf_scale] +value = 9.609 + +[mixing_length_Prandtl_number_0] +value = 0.515