From c761d97e7e472b01114acfccaba77412d772c7f0 Mon Sep 17 00:00:00 2001 From: Zhaoyi Shen <11598433+szy21@users.noreply.github.com> Date: Sat, 7 Oct 2023 19:24:38 -0700 Subject: [PATCH] add advective edmf --- .buildkite/pipeline.yml | 53 ++++ .../longrun_aquaplanet_amip.yml | 2 +- .../advective_edmfx_adv_test_box.yml | 23 ++ .../advective_edmfx_bomex_box.yml | 30 ++ .../advective_edmfx_dycoms_rf01_box.yml | 29 ++ .../advective_edmfx_gabls_box.yml | 29 ++ .../advective_edmfx_rico_box.yml | 30 ++ .../advective_edmfx_trmm_box.yml | 30 ++ .../diagnostic_edmfx_aquaplanet_tke.yml | 2 +- .../diagnostic_edmfx_bomex_box.yml | 2 +- .../diagnostic_edmfx_bomex_stretched_box.yml | 2 +- .../diagnostic_edmfx_dycoms_rf01_box.yml | 2 +- .../diagnostic_edmfx_gabls_box.yml | 2 +- .../diagnostic_edmfx_rico_box.yml | 2 +- .../diagnostic_edmfx_trmm_box.yml | 2 +- .../diagnostic_edmfx_trmm_stretched_box.yml | 2 +- config/model_configs/edmfx_bomex_box.yml | 2 +- .../model_configs/edmfx_dycoms_rf01_box.yml | 6 +- config/model_configs/edmfx_gabls_box.yml | 2 +- config/model_configs/edmfx_rico_box.yml | 2 +- config/model_configs/edmfx_trmm_box.yml | 2 +- ..._baroclinic_wave_rhoe_equilmoist_edmfx.yml | 3 +- .../flame_perf_target_diagnostic_edmfx.yml | 2 +- .../perf_configs/flame_perf_target_edmfx.yml | 2 +- examples/hybrid/driver.jl | 1 + perf/flame.jl | 6 +- src/ClimaAtmos.jl | 1 + .../advective_edmf_precomputed_quantities.jl | 258 ++++++++++++++++++ src/cache/precomputed_quantities.jl | 66 ++++- src/callbacks/callbacks.jl | 31 +++ src/initial_conditions/InitialConditions.jl | 1 + src/initial_conditions/atmos_state.jl | 21 ++ src/parameters/Parameters.jl | 3 +- src/prognostic_equations/advection.jl | 108 ++++++-- src/prognostic_equations/edmfx_closures.jl | 11 +- src/prognostic_equations/edmfx_entr_detr.jl | 88 ++++-- src/prognostic_equations/edmfx_sgs_flux.jl | 43 ++- src/prognostic_equations/hyperdiffusion.jl | 76 +++++- .../implicit/implicit_tendency.jl | 9 +- src/prognostic_equations/pressure_work.jl | 3 + src/prognostic_equations/zero_velocity.jl | 3 +- src/solver/model_getters.jl | 13 +- src/solver/types.jl | 12 +- src/utils/variable_manipulations.jl | 16 ++ toml/advective_edmfx_box.toml | 24 ++ 45 files changed, 944 insertions(+), 113 deletions(-) create mode 100644 config/model_configs/advective_edmfx_adv_test_box.yml create mode 100644 config/model_configs/advective_edmfx_bomex_box.yml create mode 100644 config/model_configs/advective_edmfx_dycoms_rf01_box.yml create mode 100644 config/model_configs/advective_edmfx_gabls_box.yml create mode 100644 config/model_configs/advective_edmfx_rico_box.yml create mode 100644 config/model_configs/advective_edmfx_trmm_box.yml create mode 100644 src/cache/advective_edmf_precomputed_quantities.jl create mode 100644 toml/advective_edmfx_box.toml diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml index 98bca456ec..2bd71bba1f 100644 --- a/.buildkite/pipeline.yml +++ b/.buildkite/pipeline.yml @@ -723,6 +723,58 @@ steps: julia --color=yes --project=examples examples/hybrid/driver.jl --config_file $CONFIG_PATH/sphere_baroclinic_wave_rhoe_equilmoist_edmfx.yml artifact_paths: "sphere_baroclinic_wave_rhoe_equilmoist_edmfx/*" + soft_fail: true + agents: + slurm_mem: 20GB + + - label: ":genie: Advective EDMFX advection test in a box" + command: > + julia --color=yes --project=examples examples/hybrid/driver.jl + --config_file $CONFIG_PATH/advective_edmfx_adv_test_box.yml + artifact_paths: "advective_edmfx_adv_test_box/*" + agents: + slurm_mem: 20GB + + - label: ":genie: Advective EDMFX GABLS in a box" + command: > + julia --color=yes --project=examples examples/hybrid/driver.jl + --config_file $CONFIG_PATH/advective_edmfx_gabls_box.yml + artifact_paths: "advective_edmfx_gabls_box/*" + agents: + slurm_mem: 20GB + + - label: ":genie: Advective EDMFX Bomex in a box" + command: > + julia --color=yes --project=examples examples/hybrid/driver.jl + --config_file $CONFIG_PATH/advective_edmfx_bomex_box.yml + artifact_paths: "advective_edmfx_bomex_box/*" + soft_fail: true + agents: + slurm_mem: 20GB + + - label: ":genie: Advective EDMFX Dycoms RF01 in a box" + command: > + julia --color=yes --project=examples examples/hybrid/driver.jl + --config_file $CONFIG_PATH/advective_edmfx_dycoms_rf01_box.yml + artifact_paths: "advective_edmfx_dycoms_rf01_box/*" + agents: + slurm_mem: 20GB + + - label: ":genie: Advective EDMFX Rico in a box" + command: > + julia --color=yes --project=examples examples/hybrid/driver.jl + --config_file $CONFIG_PATH/advective_edmfx_rico_box.yml + artifact_paths: "advective_edmfx_rico_box/*" + soft_fail: true + agents: + slurm_mem: 20GB + + - label: ":genie: Advective EDMFX TRMM in a box" + command: > + julia --color=yes --project=examples examples/hybrid/driver.jl + --config_file $CONFIG_PATH/advective_edmfx_trmm_box.yml + artifact_paths: "advective_edmfx_trmm_box/*" + soft_fail: true agents: slurm_mem: 20GB @@ -1043,6 +1095,7 @@ steps: julia --color=yes --project=perf perf/flame.jl $PERF_CONFIG_PATH/flame_perf_target_edmfx.yml artifact_paths: "flame_perf_target_edmfx/*" + soft_fail: true agents: slurm_mem: 20GB diff --git a/config/longrun_configs/longrun_aquaplanet_amip.yml b/config/longrun_configs/longrun_aquaplanet_amip.yml index 4251539b34..b7c7d20b9a 100644 --- a/config/longrun_configs/longrun_aquaplanet_amip.yml +++ b/config/longrun_configs/longrun_aquaplanet_amip.yml @@ -14,7 +14,7 @@ turbconv: "diagnostic_edmfx" prognostic_tke: true edmfx_upwinding: "first_order" edmfx_entr_model: "ConstantTimescale" -edmfx_detr_model: "ConstantCoefficient" +edmfx_detr_model: "BOverW" edmfx_nh_pressure: true edmfx_sgs_mass_flux: true edmfx_sgs_diffusive_flux: true diff --git a/config/model_configs/advective_edmfx_adv_test_box.yml b/config/model_configs/advective_edmfx_adv_test_box.yml new file mode 100644 index 0000000000..637640b61c --- /dev/null +++ b/config/model_configs/advective_edmfx_adv_test_box.yml @@ -0,0 +1,23 @@ +job_id: "advective_edmfx_adv_test_box" +initial_condition: "MoistAdiabaticProfileEDMFX" +advection_test: true +turbconv: "advective_edmfx" +ode_algo: "SSP33ShuOsher" +edmfx_upwinding: "first_order" +config: "box" +moist: "equil" +hyperdiff: "true" +kappa_4: 1e8 +x_max: 1e4 +y_max: 1e4 +z_max: 5.5e4 +x_elem: 2 +y_elem: 2 +z_elem: 63 +dz_bottom: 30.0 +dz_top: 3000.0 +dt: "10secs" +t_end: "3600secs" +dt_save_to_disk: "100secs" +FLOAT_TYPE: "Float64" +toml: [toml/edmfx_box_advection.toml] diff --git a/config/model_configs/advective_edmfx_bomex_box.yml b/config/model_configs/advective_edmfx_bomex_box.yml new file mode 100644 index 0000000000..9c611f900e --- /dev/null +++ b/config/model_configs/advective_edmfx_bomex_box.yml @@ -0,0 +1,30 @@ +job_id: "advective_edmfx_bomex_box" +initial_condition: "Bomex" +subsidence: "Bomex" +edmf_coriolis: "Bomex" +ls_adv: "Bomex" +surface_setup: "Bomex" +turbconv: "advective_edmfx" +edmfx_upwinding: first_order +edmfx_entr_model: "ConstantCoefficient" +edmfx_detr_model: "BOverW" +edmfx_sgs_mass_flux: true +edmfx_sgs_diffusive_flux: true +edmfx_nh_pressure: true +prognostic_tke: false +moist: "equil" +config: "box" +hyperdiff: "true" +kappa_4: 1.0e12 +x_max: 1e5 +y_max: 1e5 +z_max: 3e3 +x_elem: 2 +y_elem: 2 +z_elem: 60 +z_stretch: false +perturb_initstate: false +dt: "5secs" +t_end: "6hours" +dt_save_to_disk: "5mins" +toml: [toml/advective_edmfx_box.toml] diff --git a/config/model_configs/advective_edmfx_dycoms_rf01_box.yml b/config/model_configs/advective_edmfx_dycoms_rf01_box.yml new file mode 100644 index 0000000000..3c59eab1b7 --- /dev/null +++ b/config/model_configs/advective_edmfx_dycoms_rf01_box.yml @@ -0,0 +1,29 @@ +job_id: advective_edmfx_dycoms_rf01_box +initial_condition: DYCOMS_RF01 +subsidence: DYCOMS +edmf_coriolis: DYCOMS_RF01 +rad: DYCOMS_RF01 +surface_setup: DYCOMS_RF01 +turbconv: advective_edmfx +edmfx_upwinding: first_order +edmfx_entr_model: "ConstantCoefficient" +edmfx_detr_model: "BOverW" +edmfx_sgs_mass_flux: true +edmfx_sgs_diffusive_flux: true +edmfx_nh_pressure: true +prognostic_tke: false +moist: equil +config: box +hyperdiff: "true" +kappa_4: 1e12 +x_max: 1e5 +y_max: 1e5 +z_max: 1500 +x_elem: 2 +y_elem: 2 +z_elem: 30 +z_stretch: false +dt: 5secs +t_end: 4hours +dt_save_to_disk: 10mins +toml: [toml/advective_edmfx_box.toml] diff --git a/config/model_configs/advective_edmfx_gabls_box.yml b/config/model_configs/advective_edmfx_gabls_box.yml new file mode 100644 index 0000000000..50ba179be3 --- /dev/null +++ b/config/model_configs/advective_edmfx_gabls_box.yml @@ -0,0 +1,29 @@ +job_id: "advective_edmfx_gabls_box" +initial_condition: GABLS +edmf_coriolis: GABLS +surface_setup: GABLS +turbconv: "advective_edmfx" +edmfx_upwinding: "first_order" +edmfx_entr_model: "ConstantCoefficient" +edmfx_detr_model: "BOverW" +edmfx_sgs_mass_flux: true +edmfx_sgs_diffusive_flux: true +edmfx_nh_pressure: true +prognostic_tke: false +moist: "equil" +config: "box" +hyperdiff: "true" +kappa_4: 1e12 +x_max: 1e5 +y_max: 1e5 +z_max: 400 +x_elem: 2 +y_elem: 2 +z_elem: 8 +z_stretch: false +dt: "5secs" +t_end: "9hours" +dt_save_to_disk: "10mins" +perturb_initstate: false +FLOAT_TYPE: "Float64" +toml: [toml/advective_edmfx_box.toml] diff --git a/config/model_configs/advective_edmfx_rico_box.yml b/config/model_configs/advective_edmfx_rico_box.yml new file mode 100644 index 0000000000..0320db52d3 --- /dev/null +++ b/config/model_configs/advective_edmfx_rico_box.yml @@ -0,0 +1,30 @@ +job_id: "advective_edmfx_rico_box" +initial_condition: "Rico" +subsidence: "Rico" +ls_adv: "Rico" +surface_setup: "Rico" +turbconv: "advective_edmfx" +edmfx_upwinding: first_order +edmfx_entr_model: "ConstantCoefficient" +edmfx_detr_model: "BOverW" +edmfx_sgs_mass_flux: true +edmfx_sgs_diffusive_flux: true +edmfx_nh_pressure: true +prognostic_tke: false +moist: "equil" +precip_model: "0M" +config: "box" +hyperdiff: "true" +kappa_4: 1e12 +x_max: 1e5 +y_max: 1e5 +z_max: 4e3 +x_elem: 2 +y_elem: 2 +z_elem: 100 +z_stretch: false +perturb_initstate: false +dt: "5secs" +t_end: "6hours" +dt_save_to_disk: "10mins" +toml: [toml/advective_edmfx_box.toml] diff --git a/config/model_configs/advective_edmfx_trmm_box.yml b/config/model_configs/advective_edmfx_trmm_box.yml new file mode 100644 index 0000000000..fc062af52d --- /dev/null +++ b/config/model_configs/advective_edmfx_trmm_box.yml @@ -0,0 +1,30 @@ +job_id: advective_edmfx_trmm_box +initial_condition: TRMM_LBA +rad: TRMM_LBA +surface_setup: TRMM_LBA +turbconv: advective_edmfx +edmfx_upwinding: first_order +edmfx_entr_model: "ConstantCoefficient" +edmfx_detr_model: "BOverW" +edmfx_sgs_mass_flux: true +edmfx_sgs_diffusive_flux: true +edmfx_nh_pressure: true +prognostic_tke: false +moist: equil +apply_limiter: false +precip_model: "0M" +config: box +hyperdiff: "true" +kappa_4: 1e12 +x_max: 1e5 +y_max: 1e5 +z_max: 16400 +x_elem: 2 +y_elem: 2 +z_elem: 82 +z_stretch: false +dt: 5secs +t_end: 6hours +dt_save_to_disk: 10mins +FLOAT_TYPE: "Float64" +toml: [toml/advective_edmfx_box.toml] diff --git a/config/model_configs/diagnostic_edmfx_aquaplanet_tke.yml b/config/model_configs/diagnostic_edmfx_aquaplanet_tke.yml index c1c2a37bf7..a1fdf42e5d 100644 --- a/config/model_configs/diagnostic_edmfx_aquaplanet_tke.yml +++ b/config/model_configs/diagnostic_edmfx_aquaplanet_tke.yml @@ -6,7 +6,7 @@ turbconv: diagnostic_edmfx prognostic_tke: true edmfx_upwinding: first_order edmfx_entr_model: "ConstantTimescale" -edmfx_detr_model: "ConstantCoefficient" +edmfx_detr_model: "BOverW" edmfx_nh_pressure: true edmfx_sgs_mass_flux: true edmfx_sgs_diffusive_flux: true diff --git a/config/model_configs/diagnostic_edmfx_bomex_box.yml b/config/model_configs/diagnostic_edmfx_bomex_box.yml index 6f172eb58e..0aed107f8f 100644 --- a/config/model_configs/diagnostic_edmfx_bomex_box.yml +++ b/config/model_configs/diagnostic_edmfx_bomex_box.yml @@ -9,7 +9,7 @@ turbconv: "diagnostic_edmfx" prognostic_tke: true edmfx_upwinding: "first_order" edmfx_entr_model: "ConstantTimescale" -edmfx_detr_model: "ConstantCoefficient" +edmfx_detr_model: "BOverW" edmfx_nh_pressure: true edmfx_sgs_mass_flux: true edmfx_sgs_diffusive_flux: true diff --git a/config/model_configs/diagnostic_edmfx_bomex_stretched_box.yml b/config/model_configs/diagnostic_edmfx_bomex_stretched_box.yml index 23ebd68de8..9a50a12ade 100644 --- a/config/model_configs/diagnostic_edmfx_bomex_stretched_box.yml +++ b/config/model_configs/diagnostic_edmfx_bomex_stretched_box.yml @@ -9,7 +9,7 @@ turbconv: "diagnostic_edmfx" prognostic_tke: true edmfx_upwinding: "first_order" edmfx_entr_model: "ConstantTimescale" -edmfx_detr_model: "ConstantCoefficient" +edmfx_detr_model: "BOverW" edmfx_nh_pressure: true edmfx_sgs_mass_flux: true edmfx_sgs_diffusive_flux: true diff --git a/config/model_configs/diagnostic_edmfx_dycoms_rf01_box.yml b/config/model_configs/diagnostic_edmfx_dycoms_rf01_box.yml index de7af5b0f8..a7a26db7c2 100644 --- a/config/model_configs/diagnostic_edmfx_dycoms_rf01_box.yml +++ b/config/model_configs/diagnostic_edmfx_dycoms_rf01_box.yml @@ -8,7 +8,7 @@ turbconv: diagnostic_edmfx prognostic_tke: true edmfx_upwinding: first_order edmfx_entr_model: "ConstantTimescale" -edmfx_detr_model: "ConstantCoefficient" +edmfx_detr_model: "BOverW" edmfx_nh_pressure: true edmfx_sgs_mass_flux: true edmfx_sgs_diffusive_flux: true diff --git a/config/model_configs/diagnostic_edmfx_gabls_box.yml b/config/model_configs/diagnostic_edmfx_gabls_box.yml index 74c38a0b3a..4feb77db8c 100644 --- a/config/model_configs/diagnostic_edmfx_gabls_box.yml +++ b/config/model_configs/diagnostic_edmfx_gabls_box.yml @@ -6,7 +6,7 @@ turbconv: diagnostic_edmfx prognostic_tke: true edmfx_upwinding: first_order edmfx_entr_model: "ConstantTimescale" -edmfx_detr_model: "ConstantCoefficient" +edmfx_detr_model: "BOverW" edmfx_nh_pressure: true edmfx_sgs_mass_flux: true edmfx_sgs_diffusive_flux: true diff --git a/config/model_configs/diagnostic_edmfx_rico_box.yml b/config/model_configs/diagnostic_edmfx_rico_box.yml index ef7d30f2a8..7748464716 100644 --- a/config/model_configs/diagnostic_edmfx_rico_box.yml +++ b/config/model_configs/diagnostic_edmfx_rico_box.yml @@ -8,7 +8,7 @@ turbconv: diagnostic_edmfx prognostic_tke: true edmfx_upwinding: first_order edmfx_entr_model: "ConstantTimescale" -edmfx_detr_model: "ConstantCoefficient" +edmfx_detr_model: "BOverW" edmfx_nh_pressure: true edmfx_sgs_mass_flux: true edmfx_sgs_diffusive_flux: true diff --git a/config/model_configs/diagnostic_edmfx_trmm_box.yml b/config/model_configs/diagnostic_edmfx_trmm_box.yml index 2ccbc4da6c..d35e23799d 100644 --- a/config/model_configs/diagnostic_edmfx_trmm_box.yml +++ b/config/model_configs/diagnostic_edmfx_trmm_box.yml @@ -7,7 +7,7 @@ turbconv: diagnostic_edmfx prognostic_tke: true edmfx_upwinding: first_order edmfx_entr_model: "ConstantTimescale" -edmfx_detr_model: "ConstantCoefficient" +edmfx_detr_model: "BOverW" edmfx_nh_pressure: true edmfx_sgs_mass_flux: true edmfx_sgs_diffusive_flux: true diff --git a/config/model_configs/diagnostic_edmfx_trmm_stretched_box.yml b/config/model_configs/diagnostic_edmfx_trmm_stretched_box.yml index 32f585b6cc..5bbd68dc93 100644 --- a/config/model_configs/diagnostic_edmfx_trmm_stretched_box.yml +++ b/config/model_configs/diagnostic_edmfx_trmm_stretched_box.yml @@ -7,7 +7,7 @@ turbconv: diagnostic_edmfx prognostic_tke: true edmfx_upwinding: first_order edmfx_entr_model: "ConstantTimescale" -edmfx_detr_model: "ConstantCoefficient" +edmfx_detr_model: "BOverW" edmfx_nh_pressure: true edmfx_sgs_mass_flux: true edmfx_sgs_diffusive_flux: true diff --git a/config/model_configs/edmfx_bomex_box.yml b/config/model_configs/edmfx_bomex_box.yml index 8aed4c6a7c..507a8528c3 100644 --- a/config/model_configs/edmfx_bomex_box.yml +++ b/config/model_configs/edmfx_bomex_box.yml @@ -7,7 +7,7 @@ surface_setup: "Bomex" turbconv: "edmfx" edmfx_upwinding: first_order edmfx_entr_model: "ConstantCoefficient" -edmfx_detr_model: "ConstantCoefficient" +edmfx_detr_model: "BOverW" edmfx_sgs_mass_flux: true edmfx_sgs_diffusive_flux: true edmfx_nh_pressure: true diff --git a/config/model_configs/edmfx_dycoms_rf01_box.yml b/config/model_configs/edmfx_dycoms_rf01_box.yml index c59a9f7158..1da607d5ac 100644 --- a/config/model_configs/edmfx_dycoms_rf01_box.yml +++ b/config/model_configs/edmfx_dycoms_rf01_box.yml @@ -7,7 +7,7 @@ surface_setup: DYCOMS_RF01 turbconv: edmfx edmfx_upwinding: first_order edmfx_entr_model: "ConstantCoefficient" -edmfx_detr_model: "ConstantCoefficient" +edmfx_detr_model: "BOverW" edmfx_sgs_mass_flux: true edmfx_sgs_diffusive_flux: true edmfx_nh_pressure: true @@ -23,7 +23,7 @@ x_elem: 2 y_elem: 2 z_elem: 30 z_stretch: false -dt: 2secs +dt: 5secs t_end: 4hours -dt_save_to_disk: 10secs +dt_save_to_disk: 10mins toml: [toml/edmfx_box.toml] diff --git a/config/model_configs/edmfx_gabls_box.yml b/config/model_configs/edmfx_gabls_box.yml index c61d11f7b1..b01eec483b 100644 --- a/config/model_configs/edmfx_gabls_box.yml +++ b/config/model_configs/edmfx_gabls_box.yml @@ -5,7 +5,7 @@ surface_setup: GABLS turbconv: "edmfx" edmfx_upwinding: "first_order" edmfx_entr_model: "ConstantCoefficient" -edmfx_detr_model: "ConstantCoefficient" +edmfx_detr_model: "BOverW" edmfx_sgs_mass_flux: true edmfx_sgs_diffusive_flux: true edmfx_nh_pressure: true diff --git a/config/model_configs/edmfx_rico_box.yml b/config/model_configs/edmfx_rico_box.yml index af02f49bb3..634a1a91d3 100644 --- a/config/model_configs/edmfx_rico_box.yml +++ b/config/model_configs/edmfx_rico_box.yml @@ -6,7 +6,7 @@ surface_setup: "Rico" turbconv: "edmfx" edmfx_upwinding: first_order edmfx_entr_model: "ConstantCoefficient" -edmfx_detr_model: "ConstantCoefficient" +edmfx_detr_model: "BOverW" edmfx_sgs_mass_flux: true edmfx_sgs_diffusive_flux: true edmfx_nh_pressure: true diff --git a/config/model_configs/edmfx_trmm_box.yml b/config/model_configs/edmfx_trmm_box.yml index 3720fec2f3..2b790f5f31 100644 --- a/config/model_configs/edmfx_trmm_box.yml +++ b/config/model_configs/edmfx_trmm_box.yml @@ -5,7 +5,7 @@ surface_setup: TRMM_LBA turbconv: edmfx edmfx_upwinding: first_order edmfx_entr_model: "ConstantCoefficient" -edmfx_detr_model: "ConstantCoefficient" +edmfx_detr_model: "BOverW" edmfx_sgs_mass_flux: true edmfx_sgs_diffusive_flux: true edmfx_nh_pressure: true diff --git a/config/model_configs/sphere_baroclinic_wave_rhoe_equilmoist_edmfx.yml b/config/model_configs/sphere_baroclinic_wave_rhoe_equilmoist_edmfx.yml index 656cb478f7..3df6f13ff4 100644 --- a/config/model_configs/sphere_baroclinic_wave_rhoe_equilmoist_edmfx.yml +++ b/config/model_configs/sphere_baroclinic_wave_rhoe_equilmoist_edmfx.yml @@ -2,7 +2,8 @@ dt_save_to_disk: "30secs" initial_condition: "MoistBaroclinicWaveWithEDMF" dt: "1secs" edmfx_entr_model: "ConstantCoefficient" -edmfx_detr_model: "ConstantCoefficient" +edmfx_detr_model: "BOverW" +edmfx_nh_pressure: true t_end: "6mins" turbconv: "edmfx" ode_algo: "SSP33ShuOsher" diff --git a/config/perf_configs/flame_perf_target_diagnostic_edmfx.yml b/config/perf_configs/flame_perf_target_diagnostic_edmfx.yml index b5f20d92f3..e67d1b10b4 100644 --- a/config/perf_configs/flame_perf_target_diagnostic_edmfx.yml +++ b/config/perf_configs/flame_perf_target_diagnostic_edmfx.yml @@ -6,4 +6,4 @@ prognostic_tke: true edmfx_sgs_flux: true turbconv: "diagnostic_edmfx" job_id: "flame_perf_target_diagnostic_edmfx" -edmfx_detr_model: "ConstantCoefficient" +edmfx_detr_model: "BOverW" diff --git a/config/perf_configs/flame_perf_target_edmfx.yml b/config/perf_configs/flame_perf_target_edmfx.yml index f892c20bd4..0d5a514bc4 100644 --- a/config/perf_configs/flame_perf_target_edmfx.yml +++ b/config/perf_configs/flame_perf_target_edmfx.yml @@ -3,7 +3,7 @@ reference_job_id: "sphere_baroclinic_wave_rhoe_equilmoist" ode_algo: "SSP33ShuOsher" initial_condition: "MoistBaroclinicWaveWithEDMF" edmfx_entr_model: "ConstantCoefficient" -edmfx_detr_model: "ConstantCoefficient" +edmfx_detr_model: "BOverW" edmfx_sgs_flux: true turbconv: "edmfx" job_id: "flame_perf_target_edmfx" diff --git a/examples/hybrid/driver.jl b/examples/hybrid/driver.jl index 6af89053b0..552c94041a 100644 --- a/examples/hybrid/driver.jl +++ b/examples/hybrid/driver.jl @@ -40,6 +40,7 @@ reference_job_id = isnothing(ref_job_id) ? simulation.job_id : ref_job_id is_edmfx = atmos.turbconv_model isa CA.EDMFX || + atmos.turbconv_model isa CA.AdvectiveEDMFX || atmos.turbconv_model isa CA.DiagnosticEDMFX if is_edmfx && config.parsed_args["post_process"] contours_and_profiles(simulation.output_dir, reference_job_id) diff --git a/perf/flame.jl b/perf/flame.jl index 60ffe96a30..f9e83957d5 100644 --- a/perf/flame.jl +++ b/perf/flame.jl @@ -63,11 +63,11 @@ allocs_limit["flame_perf_target"] = 4656 allocs_limit["flame_perf_target_tracers"] = 204288 allocs_limit["flame_perf_target_edmfx"] = 253440 allocs_limit["flame_perf_diagnostics"] = 3016136 -allocs_limit["flame_perf_target_diagnostic_edmfx"] = 862464 -allocs_limit["flame_perf_target_edmf"] = 6382052544 +allocs_limit["flame_perf_target_diagnostic_edmfx"] = 1030080 +allocs_limit["flame_perf_target_edmf"] = 6386476224 allocs_limit["flame_perf_target_threaded"] = 5857808 allocs_limit["flame_perf_target_callbacks"] = 46407936 -allocs_limit["flame_perf_gw"] = 4844393072 +allocs_limit["flame_perf_gw"] = 4844555632 if allocs < allocs_limit[job_id] * buffer @info "TODO: lower `allocs_limit[$job_id]` to: $(allocs)" diff --git a/src/ClimaAtmos.jl b/src/ClimaAtmos.jl index 1f0c7fcac3..3daa52df33 100644 --- a/src/ClimaAtmos.jl +++ b/src/ClimaAtmos.jl @@ -26,6 +26,7 @@ include(joinpath("TurbulenceConvection_deprecated", "TurbulenceConvection.jl")) import .TurbulenceConvection as TC include(joinpath("cache", "edmf_precomputed_quantities.jl")) +include(joinpath("cache", "advective_edmf_precomputed_quantities.jl")) include(joinpath("cache", "diagnostic_edmf_precomputed_quantities.jl")) include(joinpath("cache", "precomputed_quantities.jl")) diff --git a/src/cache/advective_edmf_precomputed_quantities.jl b/src/cache/advective_edmf_precomputed_quantities.jl new file mode 100644 index 0000000000..1d78f834ea --- /dev/null +++ b/src/cache/advective_edmf_precomputed_quantities.jl @@ -0,0 +1,258 @@ +##### +##### Precomputed quantities +##### +import Thermodynamics as TD +import ClimaCore: Spaces, Fields + +""" + set_advective_edmf_precomputed_quantities!(Y, p, t) + +Updates the precomputed quantities stored in `p` for edmfx. +""" +function set_advective_edmf_precomputed_quantities!(Y, p, ᶠuₕ³, t) + (; energy_form, moisture_model, turbconv_model) = p.atmos + #EDMFX BCs only support total energy as state variable + @assert energy_form isa TotalEnergy + @assert !(moisture_model isa DryModel) + + FT = Spaces.undertype(axes(Y.c)) + (; params) = p + (; dt) = p.simulation + thermo_params = CAP.thermodynamics_params(params) + n = n_mass_flux_subdomains(turbconv_model) + thermo_args = (thermo_params, energy_form, moisture_model) + + (; ᶜspecific, ᶜp, ᶜΦ, ᶜh_tot, ᶜρ_ref) = p + (; ᶜtke⁰, ᶜρa⁰, ᶠu₃⁰, ᶜu⁰, ᶠu³⁰, ᶜK⁰, ᶜts⁰, ᶜρ⁰, ᶜh_tot⁰, ᶜq_tot⁰) = p + (; ᶜmixing_length, ᶜlinear_buoygrad, ᶜstrain_rate_norm, ᶜK_u, ᶜK_h) = p + (; ᶜuʲs, ᶠu³ʲs, ᶜKʲs, ᶜtsʲs, ᶜρʲs, ᶜentrʲs, ᶜdetrʲs) = p + (; ustar, obukhov_length, buoyancy_flux) = p.sfc_conditions + + @. ᶜρa⁰ = ρa⁰(Y.c) + @. ᶜtke⁰ = divide_by_ρa(Y.c.sgs⁰.ρatke, ᶜρa⁰, 0, Y.c.ρ, turbconv_model) + @. ᶜh_tot⁰ = divide_by_ρa( + Y.c.ρ * ᶜh_tot - ρah_tot⁺(Y.c.sgsʲs), + ᶜρa⁰, + Y.c.ρ * ᶜh_tot, + Y.c.ρ, + turbconv_model, + ) + @. ᶜq_tot⁰ = divide_by_ρa( + Y.c.ρq_tot - ρaq_tot⁺(Y.c.sgsʲs), + ᶜρa⁰, + Y.c.ρq_tot, + Y.c.ρ, + turbconv_model, + ) + set_sgs_ᶠu₃!(u₃⁰, ᶠu₃⁰, Y, turbconv_model) + set_velocity_quantities!(ᶜu⁰, ᶠu³⁰, ᶜK⁰, ᶠu₃⁰, Y.c.uₕ, ᶠuₕ³) + @. ᶜK⁰ += ᶜtke⁰ + @. ᶜts⁰ = TD.PhaseEquil_phq(thermo_params, ᶜp, ᶜh_tot⁰ - ᶜK⁰ - ᶜΦ, ᶜq_tot⁰) + @. ᶜρ⁰ = TD.air_density(thermo_params, ᶜts⁰) + + for j in 1:n + ᶜuʲ = ᶜuʲs.:($j) + ᶠu³ʲ = ᶠu³ʲs.:($j) + ᶜKʲ = ᶜKʲs.:($j) + ᶠu₃ʲ = Y.f.sgsʲs.:($j).u₃ + ᶜtsʲ = ᶜtsʲs.:($j) + ᶜρʲ = ᶜρʲs.:($j) + ᶜh_totʲ = Y.c.sgsʲs.:($j).h_tot + ᶜq_totʲ = Y.c.sgsʲs.:($j).q_tot + + set_velocity_quantities!(ᶜuʲ, ᶠu³ʲ, ᶜKʲ, ᶠu₃ʲ, Y.c.uₕ, ᶠuₕ³) + @. ᶜtsʲ = + TD.PhaseEquil_phq(thermo_params, ᶜp, ᶜh_totʲ - ᶜKʲ - ᶜΦ, ᶜq_totʲ) + @. ᶜρʲ = TD.air_density(thermo_params, ᶜtsʲ) + + # EDMFX boundary condition: + + # We need field_values everywhere because we are mixing + # information from surface and first interior inside the + # sgs_h/q_tot_first_interior_bc call. + ᶜz_int_val = + Fields.field_values(Fields.level(Fields.coordinate_field(Y.c).z, 1)) + z_sfc_val = Fields.field_values( + Fields.level(Fields.coordinate_field(Y.f).z, Fields.half), + ) + ᶜρ_int_val = Fields.field_values(Fields.level(Y.c.ρ, 1)) + ᶜp_int_val = Fields.field_values(Fields.level(ᶜp, 1)) + (; ρ_flux_h_tot, ρ_flux_q_tot, ustar, obukhov_length) = p.sfc_conditions + buoyancy_flux_val = Fields.field_values(buoyancy_flux) + ρ_flux_h_tot_val = Fields.field_values(ρ_flux_h_tot) + ρ_flux_q_tot_val = Fields.field_values(ρ_flux_q_tot) + ustar_val = Fields.field_values(ustar) + obukhov_length_val = Fields.field_values(obukhov_length) + sfc_local_geometry_val = Fields.field_values( + Fields.local_geometry_field(Fields.level(Y.f, Fields.half)), + ) + + # Based on boundary conditions for updrafts we overwrite + # the first interior point for EDMFX ᶜh_totʲ... + ᶜh_tot_int_val = Fields.field_values(Fields.level(ᶜh_tot, 1)) + ᶜh_totʲ_int_val = Fields.field_values(Fields.level(ᶜh_totʲ, 1)) + @. ᶜh_totʲ_int_val = sgs_scalar_first_interior_bc( + ᶜz_int_val - z_sfc_val, + ᶜρ_int_val, + ᶜh_tot_int_val, + buoyancy_flux_val, + ρ_flux_h_tot_val, + ustar_val, + obukhov_length_val, + sfc_local_geometry_val, + ) + + # ... and the first interior point for EDMFX ᶜq_totʲ. + ᶜq_tot_int_val = Fields.field_values(Fields.level(ᶜspecific.q_tot, 1)) + ᶜq_totʲ_int_val = Fields.field_values(Fields.level(ᶜq_totʲ, 1)) + @. ᶜq_totʲ_int_val = sgs_scalar_first_interior_bc( + ᶜz_int_val - z_sfc_val, + ᶜρ_int_val, + ᶜq_tot_int_val, + buoyancy_flux_val, + ρ_flux_q_tot_val, + ustar_val, + obukhov_length_val, + sfc_local_geometry_val, + ) + + # Then overwrite the prognostic variables at first inetrior point. + ᶜKʲ_int_val = Fields.field_values(Fields.level(ᶜKʲ, 1)) + ᶜΦ_int_val = Fields.field_values(Fields.level(ᶜΦ, 1)) + ᶜtsʲ_int_val = Fields.field_values(Fields.level(ᶜtsʲ, 1)) + @. ᶜtsʲ_int_val = TD.PhaseEquil_phq( + thermo_params, + ᶜp_int_val, + ᶜh_totʲ_int_val - ᶜKʲ_int_val - ᶜΦ_int_val, + ᶜq_totʲ_int_val, + ) + sgsʲs_ρ_int_val = Fields.field_values(Fields.level(ᶜρʲs.:($j), 1)) + sgsʲs_ρa_int_val = + Fields.field_values(Fields.level(Y.c.sgsʲs.:($j).ρa, 1)) + + turbconv_params = CAP.turbconv_params(params) + @. sgsʲs_ρa_int_val = + $(FT(turbconv_params.surface_area)) * + TD.air_density(thermo_params, ᶜtsʲ_int_val) + end + + ᶜz = Fields.coordinate_field(Y.c).z + z_sfc = Fields.level(Fields.coordinate_field(Y.f).z, Fields.half) + ᶜdz = Fields.Δz_field(axes(Y.c)) + ᶜlg = Fields.local_geometry_field(Y.c) + + for j in 1:n + @. ᶜentrʲs.:($$j) = entrainment( + params, + ᶜz, + 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, ᶜρ_ref, ᶜρʲs.:($$j)), + get_physical_w(ᶜu⁰, ᶜlg), + TD.relative_humidity(thermo_params, ᶜts⁰), + ᶜphysical_buoyancy(params, ᶜρ_ref, ᶜρ⁰), + dt, + p.atmos.edmfx_entr_model, + ) + @. ᶜdetrʲs.:($$j) = detrainment( + params, + ᶜz, + 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, ᶜρ_ref, ᶜρʲs.:($$j)), + get_physical_w(ᶜu⁰, ᶜlg), + TD.relative_humidity(thermo_params, ᶜts⁰), + ᶜphysical_buoyancy(params, ᶜρ_ref, ᶜρ⁰), + dt, + p.atmos.edmfx_detr_model, + ) + end + + # First order approximation: Use environmental mean fields. + @. ᶜlinear_buoygrad = buoyancy_gradients( + params, + moisture_model, + EnvBuoyGrad( + BuoyGradMean(), + TD.air_temperature(thermo_params, ᶜts⁰), # t_sat + TD.vapor_specific_humidity(thermo_params, ᶜts⁰), # qv_sat + ᶜq_tot⁰, # qt_sat + TD.dry_pottemp(thermo_params, ᶜts⁰), # θ_sat + TD.liquid_ice_pottemp(thermo_params, ᶜts⁰), # θ_liq_ice_sat + projected_vector_data( + C3, + ᶜgradᵥ(ᶠinterp(TD.virtual_pottemp(thermo_params, ᶜts⁰))), + ᶜlg, + ), # ∂θv∂z_unsat + projected_vector_data(C3, ᶜgradᵥ(ᶠinterp(ᶜq_tot⁰)), ᶜlg), # ∂qt∂z_sat + projected_vector_data( + C3, + ᶜgradᵥ(ᶠinterp(TD.liquid_ice_pottemp(thermo_params, ᶜts⁰))), + ᶜlg, + ), # ∂θl∂z_sat + ᶜp, # p + ifelse(TD.has_condensate(thermo_params, ᶜts⁰), 1, 0), # en_cld_frac + ᶜρ⁰, # ρ + ), + ) + + # TODO: Currently the shear production only includes vertical gradients + ᶠu⁰ = p.ᶠtemp_C123 + @. ᶠu⁰ = C123(ᶠinterp(Y.c.uₕ)) + C123(ᶠu³⁰) + ᶜstrain_rate = p.ᶜtemp_UVWxUVW + compute_strain_rate_center!(ᶜstrain_rate, ᶠu⁰) + @. ᶜstrain_rate_norm = norm_sqr(ᶜstrain_rate) + + ᶜprandtl_nvec = p.ᶜtemp_scalar + @. ᶜprandtl_nvec = turbulent_prandtl_number( + params, + obukhov_length, + ᶜlinear_buoygrad, + ᶜstrain_rate_norm, + ) + ᶜtke_exch = p.ᶜtemp_scalar_2 + @. ᶜtke_exch = 0 + for j in 1:n + @. ᶜtke_exch += + Y.c.sgsʲs.:($$j).ρa * ᶜdetrʲs.:($$j) / ᶜρa⁰ * ( + 1 / 2 * + ( + get_physical_w(ᶜuʲs.:($$j), ᶜlg) - get_physical_w(ᶜu⁰, ᶜlg) + )^2 - ᶜtke⁰ + ) + end + + sfc_tke = Fields.level(ᶜtke⁰, 1) + @. ᶜmixing_length = mixing_length( + p.params, + ustar, + ᶜz, + z_sfc, + ᶜdz, + sfc_tke, + ᶜlinear_buoygrad, + ᶜtke⁰, + obukhov_length, + ᶜstrain_rate_norm, + ᶜprandtl_nvec, + ᶜtke_exch, + ) + + turbconv_params = CAP.turbconv_params(params) + c_m = TCP.tke_ed_coeff(turbconv_params) + @. ᶜK_u = c_m * ᶜmixing_length * sqrt(max(ᶜtke⁰, 0)) + # TODO: add Prantdl number + @. ᶜK_h = ᶜK_u / ᶜprandtl_nvec + + return nothing +end diff --git a/src/cache/precomputed_quantities.jl b/src/cache/precomputed_quantities.jl index b3be3719ea..0349a4b089 100644 --- a/src/cache/precomputed_quantities.jl +++ b/src/cache/precomputed_quantities.jl @@ -40,6 +40,10 @@ function precomputed_quantities(Y, atmos) !(atmos.moisture_model isa DryModel) && atmos.energy_form isa TotalEnergy ) || !(atmos.turbconv_model isa DiagnosticEDMFX) + @assert ( + !(atmos.moisture_model isa DryModel) && + atmos.energy_form isa TotalEnergy + ) || !(atmos.turbconv_model isa AdvectiveEDMFX) TST = thermo_state_type(atmos.moisture_model, FT) SCT = SurfaceConditions.surface_conditions_type(atmos, FT) n = n_mass_flux_subdomains(atmos.turbconv_model) @@ -88,6 +92,32 @@ function precomputed_quantities(Y, atmos) ) : (;) )..., ) : (;) + advective_sgs_quantities = + atmos.turbconv_model isa AdvectiveEDMFX ? + (; + ᶜtke⁰ = similar(Y.c, FT), + ᶜρa⁰ = similar(Y.c, FT), + ᶠu₃⁰ = similar(Y.f, C3{FT}), + ᶜu⁰ = similar(Y.c, C123{FT}), + ᶠu³⁰ = similar(Y.f, CT3{FT}), + ᶜK⁰ = similar(Y.c, FT), + ᶜh_tot⁰ = similar(Y.c, FT), + ᶜq_tot⁰ = similar(Y.c, FT), + ᶜts⁰ = similar(Y.c, TST), + ᶜρ⁰ = similar(Y.c, FT), + ᶜlinear_buoygrad = similar(Y.c, FT), + ᶜstrain_rate_norm = similar(Y.c, FT), + ᶜmixing_length = similar(Y.c, FT), + ᶜK_u = similar(Y.c, FT), + ᶜK_h = similar(Y.c, FT), + ᶜuʲs = similar(Y.c, NTuple{n, C123{FT}}), + ᶠu³ʲs = similar(Y.f, NTuple{n, CT3{FT}}), + ᶜKʲs = similar(Y.c, NTuple{n, FT}), + ᶜtsʲs = similar(Y.c, NTuple{n, TST}), + ᶜρʲs = similar(Y.c, NTuple{n, FT}), + ᶜentrʲs = similar(Y.c, NTuple{n, FT}), + ᶜdetrʲs = similar(Y.c, NTuple{n, FT}), + ) : (;) diagnostic_sgs_quantities = atmos.turbconv_model isa DiagnosticEDMFX ? (; @@ -117,7 +147,12 @@ function precomputed_quantities(Y, atmos) ᶜK_h = similar(Y.c, FT), ρatke_flux = similar(Fields.level(Y.f, half), C3{FT}), ) : (;) - return (; gs_quantities..., sgs_quantities..., diagnostic_sgs_quantities...) + return (; + gs_quantities..., + sgs_quantities..., + advective_sgs_quantities..., + diagnostic_sgs_quantities..., + ) end # Interpolates the third contravariant component of Y.c.uₕ to cell faces. @@ -140,7 +175,7 @@ function set_velocity_at_surface!(Y, ᶠuₕ³, turbconv_model) sfc_uₕ³ = Fields.level(ᶠuₕ³.components.data.:1, half) sfc_g³³ = g³³_field(sfc_u₃) @. sfc_u₃ = -sfc_uₕ³ / sfc_g³³ # u³ = uₕ³ + w³ = uₕ³ + w₃ * g³³ - if turbconv_model isa EDMFX + if turbconv_model isa EDMFX || turbconv_model isa AdvectiveEDMFX for j in 1:n_mass_flux_subdomains(turbconv_model) sfc_u₃ʲ = Fields.level(Y.f.sgsʲs.:($j).u₃.components.data.:1, half) @. sfc_u₃ʲ = sfc_u₃ @@ -160,7 +195,7 @@ function set_velocity_at_top!(Y, turbconv_model) Spaces.nlevels(axes(Y.c)) + half, ) @. top_u₃ = 0 - if turbconv_model isa EDMFX + if turbconv_model isa EDMFX || turbconv_model isa AdvectiveEDMFX for j in 1:n_mass_flux_subdomains(turbconv_model) top_u₃ʲ = Fields.level( Y.f.sgsʲs.:($j).u₃.components.data.:1, @@ -341,6 +376,10 @@ NVTX.@annotate function set_precomputed_quantities!(Y, p, t) set_edmf_precomputed_quantities!(Y, p, ᶠuₕ³, t) end + if turbconv_model isa AdvectiveEDMFX + set_advective_edmf_precomputed_quantities!(Y, p, ᶠuₕ³, t) + end + if turbconv_model isa DiagnosticEDMFX set_diagnostic_edmf_precomputed_quantities!(Y, p, t) end @@ -373,6 +412,27 @@ function output_sgs_quantities(Y, p, t) return (; ᶜspecific⁺, ᶠu₃⁺, ᶜu⁺, ᶠu³⁺, ᶜK⁺, ᶜts⁺, ᶜa⁺, ᶜa⁰) end +""" + output_advective_sgs_quantities(Y, p, t) + +Sets `ᶜu⁺`, `ᶠu³⁺`, `ᶜts⁺` and `ᶜa⁺` to be the same as the +values of the first updraft. +""" +function output_advective_sgs_quantities(Y, p, t) + (; turbconv_model) = p.atmos + thermo_params = CAP.thermodynamics_params(p.params) + (; ᶜp, ᶜρa⁰, ᶜρ⁰, ᶜΦ, ᶜtsʲs) = p + ᶠuₕ³ = p.ᶠtemp_CT3 + set_ᶠuₕ³!(ᶠuₕ³, Y) + (ᶠu₃⁺, ᶜu⁺, ᶠu³⁺, ᶜK⁺) = similar.((p.ᶠu₃⁰, p.ᶜu⁰, p.ᶠu³⁰, p.ᶜK⁰)) + set_sgs_ᶠu₃!(u₃⁺, ᶠu₃⁺, Y, turbconv_model) + set_velocity_quantities!(ᶜu⁺, ᶠu³⁺, ᶜK⁺, ᶠu₃⁺, Y.c.uₕ, ᶠuₕ³) + ᶜts⁺ = ᶜtsʲs.:1 + ᶜa⁺ = @. draft_area(ρa⁺(Y.c), TD.air_density(thermo_params, ᶜts⁺)) + ᶜa⁰ = @. draft_area(ᶜρa⁰, ᶜρ⁰) + return (; ᶠu₃⁺, ᶜu⁺, ᶠu³⁺, ᶜK⁺, ᶜts⁺, ᶜa⁺, ᶜa⁰) +end + """ output_diagnostic_sgs_quantities(Y, p, t) diff --git a/src/callbacks/callbacks.jl b/src/callbacks/callbacks.jl index 4625e50872..4112e8d7e7 100644 --- a/src/callbacks/callbacks.jl +++ b/src/callbacks/callbacks.jl @@ -304,6 +304,37 @@ NVTX.@annotate function compute_diagnostics(integrator) ᶜts⁰, ) .+ ᶜa⁺ .* cloud_fraction.(thermo_params, ᶜts⁺, ᶜa⁺), ) + elseif p.atmos.turbconv_model isa AdvectiveEDMFX + (; ᶜtke⁰, ᶜu⁰, ᶜts⁰, ᶜmixing_length) = p + (; ᶜu⁺, ᶜts⁺, ᶜa⁺, ᶜa⁰) = output_advective_sgs_quantities(Y, p, t) + env_diagnostics = (; + common_diagnostics(p, ᶜu⁰, ᶜts⁰)..., + area = ᶜa⁰, + cloud_fraction = get_cloud_fraction.( + thermo_params, + env_thermo_quad, + ᶜp, + ᶜts⁰, + ), + tke = ᶜtke⁰, + mixing_length = ᶜmixing_length, + ) + draft_diagnostics = (; + common_diagnostics(p, ᶜu⁺, ᶜts⁺)..., + area = ᶜa⁺, + cloud_fraction = cloud_fraction.(thermo_params, ᶜts⁺, ᶜa⁺), + ) + turbulence_convection_diagnostic = (; + add_prefix(env_diagnostics, :env_)..., + add_prefix(draft_diagnostics, :draft_)..., + cloud_fraction = ᶜa⁰ .* + get_cloud_fraction.( + thermo_params, + env_thermo_quad, + ᶜp, + ᶜts⁰, + ) .+ ᶜa⁺ .* cloud_fraction.(thermo_params, ᶜts⁺, ᶜa⁺), + ) elseif p.atmos.turbconv_model isa DiagnosticEDMFX (; ᶜtke⁰, ᶜmixing_length) = p (; ᶜu⁺, ᶜts⁺, ᶜa⁺) = output_diagnostic_sgs_quantities(Y, p, t) diff --git a/src/initial_conditions/InitialConditions.jl b/src/initial_conditions/InitialConditions.jl index 29cef59fd7..081bc892a1 100644 --- a/src/initial_conditions/InitialConditions.jl +++ b/src/initial_conditions/InitialConditions.jl @@ -16,6 +16,7 @@ import ..PrognosticSurfaceTemperature import ..C3 import ..C12 import ..EDMFX +import ..AdvectiveEDMFX import ..DiagnosticEDMFX import ..n_mass_flux_subdomains import ..gs_to_sgs diff --git a/src/initial_conditions/atmos_state.jl b/src/initial_conditions/atmos_state.jl index 84865dd0d2..8caf9e2170 100644 --- a/src/initial_conditions/atmos_state.jl +++ b/src/initial_conditions/atmos_state.jl @@ -129,6 +129,7 @@ function turbconv_center_variables(ls, turbconv_model::TC.EDMFModel, gs_vars) up = ntuple(_ -> (; ρarea = ρa, ρae_tot, ρaq_tot), Val(n)) return (; turbconv = (; en, up)) end + function turbconv_center_variables(ls, turbconv_model::EDMFX, gs_vars) n = n_mass_flux_subdomains(turbconv_model) a_draft = ls.turbconv_state.draft_area @@ -137,6 +138,20 @@ function turbconv_center_variables(ls, turbconv_model::EDMFX, gs_vars) return (; sgs⁰, sgsʲs) end +function turbconv_center_variables(ls, turbconv_model::AdvectiveEDMFX, gs_vars) + n = n_mass_flux_subdomains(turbconv_model) + a_draft = ls.turbconv_state.draft_area + sgs⁰ = (; ρatke = ls.ρ * (1 - a_draft) * ls.turbconv_state.tke) + ρa = ls.ρ * a_draft / n + h_tot = + TD.specific_enthalpy(ls.thermo_params, ls.thermo_state) + + norm_sqr(ls.velocity) / 2 + + CAP.grav(ls.params) * ls.geometry.coordinates.z + q_tot = TD.total_specific_humidity(ls.thermo_params, ls.thermo_state) + sgsʲs = ntuple(_ -> (; ρa = ρa, h_tot = h_tot, q_tot = q_tot), Val(n)) + return (; sgs⁰, sgsʲs) +end + function turbconv_center_variables(ls, turbconv_model::DiagnosticEDMFX, gs_vars) sgs⁰ = (; ρatke = ls.ρ * ls.turbconv_state.tke) return (; sgs⁰) @@ -157,4 +172,10 @@ turbconv_face_variables(ls, turbconv_model::EDMFX) = (; Val(n_mass_flux_subdomains(turbconv_model)), ) ) +turbconv_face_variables(ls, turbconv_model::AdvectiveEDMFX) = (; + sgsʲs = ntuple( + _ -> (; u₃ = C3(ls.turbconv_state.velocity, ls.geometry)), + Val(n_mass_flux_subdomains(turbconv_model)), + ) +) turbconv_face_variables(ls, turbconv_model::DiagnosticEDMFX) = (;) diff --git a/src/parameters/Parameters.jl b/src/parameters/Parameters.jl index 851513c8b2..e3c32b00bf 100644 --- a/src/parameters/Parameters.jl +++ b/src/parameters/Parameters.jl @@ -22,9 +22,10 @@ Base.@kwdef struct ClimaAtmosParameters{FT, TP, RP, IP, MPP, SFP, TCP} <: ACAP planet_radius::FT astro_unit::FT entr_tau::FT - detr_tau::FT entr_coeff::FT + detr_tau::FT detr_coeff::FT + detr_buoy_coeff::FT C_E::FT C_H::FT # Held Suarez diff --git a/src/prognostic_equations/advection.jl b/src/prognostic_equations/advection.jl index 94d4bbfb49..0cb082fad5 100644 --- a/src/prognostic_equations/advection.jl +++ b/src/prognostic_equations/advection.jl @@ -12,12 +12,14 @@ NVTX.@annotate function horizontal_advection_tendency!(Yₜ, Y, p, t) if p.atmos.turbconv_model isa AbstractEDMF (; ᶜu⁰) = p end - if p.atmos.turbconv_model isa EDMFX + if p.atmos.turbconv_model isa EDMFX || + p.atmos.turbconv_model isa AdvectiveEDMFX (; ᶜuʲs) = p end @. Yₜ.c.ρ -= wdivₕ(Y.c.ρ * ᶜu) - if p.atmos.turbconv_model isa EDMFX + if p.atmos.turbconv_model isa EDMFX || + p.atmos.turbconv_model isa AdvectiveEDMFX for j in 1:n @. Yₜ.c.sgsʲs.:($$j).ρa -= wdivₕ(Y.c.sgsʲs.:($$j).ρa * ᶜuʲs.:($$j)) end @@ -42,6 +44,13 @@ NVTX.@annotate function horizontal_advection_tendency!(Yₜ, Y, p, t) end end + if p.atmos.turbconv_model isa AdvectiveEDMFX + for j in 1:n + @. Yₜ.c.sgsʲs.:($$j).h_tot -= + adjoint(CA.CT12(Y.c.uₕ)) * CA.gradₕ(Y.c.sgsʲs.:($$j).h_tot) + end + end + if use_prognostic_tke(p.atmos.turbconv_model) @. Yₜ.c.sgs⁰.ρatke -= wdivₕ(Y.c.sgs⁰.ρatke * ᶜu⁰) end @@ -71,6 +80,14 @@ NVTX.@annotate function horizontal_tracer_advection_tendency!(Yₜ, Y, p, t) end end end + + if p.atmos.turbconv_model isa AdvectiveEDMFX + for j in 1:n + @. Yₜ.c.sgsʲs.:($$j).q_tot -= + adjoint(CA.CT12(Y.c.uₕ)) * CA.gradₕ(Y.c.sgsʲs.:($$j).q_tot) + end + end + return nothing end @@ -84,13 +101,16 @@ NVTX.@annotate function explicit_vertical_advection_tendency!(Yₜ, Y, p, t) ᶜJ = Fields.local_geometry_field(Y.c).J (; ᶜu, ᶠu³, ᶜK, ᶜf) = p (; edmfx_upwinding) = n > 0 || advect_tke ? p : all_nothing - (; ᶜuʲs, ᶠu³ʲs, ᶜKʲs, ᶜρʲs, ᶜspecificʲs) = n > 0 ? p : all_nothing + (; ᶜuʲs, ᶠu³ʲs, ᶜKʲs, ᶜρʲs) = n > 0 ? p : all_nothing + (; ᶜspecificʲs) = turbconv_model isa EDMFX ? p : all_nothing (; ᶜp, ᶜp_ref, ᶜρ_ref, ᶠgradᵥ_ᶜΦ) = n > 0 ? p : all_nothing - (; ᶜh_totʲs) = n > 0 && is_total_energy ? p : all_nothing + (; ᶜh_totʲs) = turbconv_model isa EDMFX && is_total_energy ? p : all_nothing (; ᶠu³⁰) = advect_tke ? p : all_nothing ᶜρa⁰ = advect_tke ? (n > 0 ? p.ᶜρa⁰ : Y.c.ρ) : nothing ᶜρ⁰ = advect_tke ? (n > 0 ? p.ᶜρ⁰ : Y.c.ρ) : nothing - ᶜtke⁰ = advect_tke ? (n > 0 ? p.ᶜspecific⁰.tke : p.ᶜtke⁰) : nothing + ᶜtke⁰ = + advect_tke ? (turbconv_model isa EDMFX ? p.ᶜspecific⁰.tke : p.ᶜtke⁰) : + nothing ᶜa_scalar = p.ᶜtemp_scalar ᶜω³ = p.ᶜtemp_CT3 ᶠω¹² = p.ᶠtemp_CT12 @@ -139,25 +159,12 @@ NVTX.@annotate function explicit_vertical_advection_tendency!(Yₜ, Y, p, t) end # TODO: Move this to implicit_vertical_advection_tendency!. - for j in 1:n - @. ᶜa_scalar[colidx] = - draft_area(Y.c.sgsʲs.:($$j).ρa[colidx], ᶜρʲs.:($$j)[colidx]) - vertical_transport!( - Yₜ.c.sgsʲs.:($j).ρa[colidx], - ᶜJ[colidx], - ᶜρʲs.:($j)[colidx], - ᶠu³ʲs.:($j)[colidx], - ᶜa_scalar[colidx], - dt, - edmfx_upwinding, - ) - - if :ρae_tot in propertynames(Yₜ.c.sgsʲs.:($j)) + if p.atmos.turbconv_model isa EDMFX + for j in 1:n @. ᶜa_scalar[colidx] = - ᶜh_totʲs.:($$j)[colidx] * draft_area(Y.c.sgsʲs.:($$j).ρa[colidx], ᶜρʲs.:($$j)[colidx]) vertical_transport!( - Yₜ.c.sgsʲs.:($j).ρae_tot[colidx], + Yₜ.c.sgsʲs.:($j).ρa[colidx], ᶜJ[colidx], ᶜρʲs.:($j)[colidx], ᶠu³ʲs.:($j)[colidx], @@ -165,16 +172,51 @@ NVTX.@annotate function explicit_vertical_advection_tendency!(Yₜ, Y, p, t) dt, edmfx_upwinding, ) + + if :ρae_tot in propertynames(Yₜ.c.sgsʲs.:($j)) + @. ᶜa_scalar[colidx] = + ᶜh_totʲs.:($$j)[colidx] * draft_area( + Y.c.sgsʲs.:($$j).ρa[colidx], + ᶜρʲs.:($$j)[colidx], + ) + vertical_transport!( + Yₜ.c.sgsʲs.:($j).ρae_tot[colidx], + ᶜJ[colidx], + ᶜρʲs.:($j)[colidx], + ᶠu³ʲs.:($j)[colidx], + ᶜa_scalar[colidx], + dt, + edmfx_upwinding, + ) + end + + for (ᶜρaχʲₜ, ᶜχʲ, χ_name) in + matching_subfields(Yₜ.c.sgsʲs.:($j), ᶜspecificʲs.:($j)) + χ_name == :e_tot && continue + @. ᶜa_scalar[colidx] = + ᶜχʲ[colidx] * draft_area( + Y.c.sgsʲs.:($$j).ρa[colidx], + ᶜρʲs.:($$j)[colidx], + ) + vertical_transport!( + ᶜρaχʲₜ[colidx], + ᶜJ[colidx], + ᶜρʲs.:($j)[colidx], + ᶠu³ʲs.:($j)[colidx], + ᶜa_scalar[colidx], + dt, + edmfx_upwinding, + ) + end end + end - for (ᶜρaχʲₜ, ᶜχʲ, χ_name) in - matching_subfields(Yₜ.c.sgsʲs.:($j), ᶜspecificʲs.:($j)) - χ_name == :e_tot && continue + if p.atmos.turbconv_model isa AdvectiveEDMFX + for j in 1:n @. ᶜa_scalar[colidx] = - ᶜχʲ[colidx] * draft_area(Y.c.sgsʲs.:($$j).ρa[colidx], ᶜρʲs.:($$j)[colidx]) vertical_transport!( - ᶜρaχʲₜ[colidx], + Yₜ.c.sgsʲs.:($j).ρa[colidx], ᶜJ[colidx], ᶜρʲs.:($j)[colidx], ᶠu³ʲs.:($j)[colidx], @@ -182,6 +224,20 @@ NVTX.@annotate function explicit_vertical_advection_tendency!(Yₜ, Y, p, t) dt, edmfx_upwinding, ) + + vertical_advection!( + Yₜ.c.sgsʲs.:($j).h_tot[colidx], + ᶠu³ʲs.:($j)[colidx], + Y.c.sgsʲs.:($j).h_tot[colidx], + edmfx_upwinding, + ) + + vertical_advection!( + Yₜ.c.sgsʲs.:($j).q_tot[colidx], + ᶠu³ʲs.:($j)[colidx], + Y.c.sgsʲs.:($j).q_tot[colidx], + edmfx_upwinding, + ) end end diff --git a/src/prognostic_equations/edmfx_closures.jl b/src/prognostic_equations/edmfx_closures.jl index 31d332414b..ffe8394046 100644 --- a/src/prognostic_equations/edmfx_closures.jl +++ b/src/prognostic_equations/edmfx_closures.jl @@ -76,10 +76,17 @@ end edmfx_nh_pressure_cache(Y, turbconv_model) = (;) edmfx_nh_pressure_tendency!(Yₜ, Y, p, t, colidx, turbconv_model) = nothing -function edmfx_nh_pressure_tendency!(Yₜ, Y, p, t, colidx, turbconv_model::EDMFX) +function edmfx_nh_pressure_tendency!( + Yₜ, + Y, + p, + t, + colidx, + turbconv_model::Union{EDMFX, AdvectiveEDMFX}, +) n = n_mass_flux_subdomains(turbconv_model) - (; params, ᶜρʲs, ᶜρ_ref, ᶜp, ᶜp_ref, ᶠgradᵥ_ᶜΦ, ᶠu₃⁰) = p + (; params, ᶜρʲs, ᶠgradᵥ_ᶜΦ, ᶠu₃⁰) = p FT = eltype(Y) ᶜz = Fields.coordinate_field(Y.c).z z_sfc = Fields.level(Fields.coordinate_field(Y.f).z, Fields.half) diff --git a/src/prognostic_equations/edmfx_entr_detr.jl b/src/prognostic_equations/edmfx_entr_detr.jl index 44ceaf99b5..bdb41b6c43 100644 --- a/src/prognostic_equations/edmfx_entr_detr.jl +++ b/src/prognostic_equations/edmfx_entr_detr.jl @@ -211,7 +211,7 @@ function detrainment( g = CAP.grav(params) turbconv_params = CAP.turbconv_params(params) ᶜaʲ_max = TCP.max_area(turbconv_params) - max_area_limiter = 0.1 * exp(-10 * (ᶜaʲ_max - ᶜaʲ)) + 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) @@ -254,10 +254,18 @@ function detrainment( ᶜRH⁰::FT, ᶜbuoy⁰::FT, dt::FT, - ::ConstantCoefficientDetrainment, + ::BOverWDetrainment, ) where {FT} detr_coeff = CAP.detr_coeff(params) - detr = min(detr_coeff * abs(ᶜwʲ), 1 / dt) + detr_buoy_coeff = CAP.detr_buoy_coeff(params) + detr = min( + abs(ᶜwʲ) * ( + detr_coeff + + detr_buoy_coeff * abs(min(ᶜbuoyʲ - ᶜbuoy⁰, 0)) / + max(eps(FT), (ᶜwʲ - ᶜw⁰) * (ᶜwʲ - ᶜw⁰)) + ), + 1 / dt, + ) return detr end @@ -283,6 +291,28 @@ function detrainment( return detr * FT(2) * hm_limiter(ᶜaʲ) 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, + ::ConstantTimescaleDetrainment, +) where {FT} + detr_tau = CAP.detr_tau(params) + detr = min(1 / detr_tau, 1 / dt) + return detr +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) @@ -315,24 +345,36 @@ function edmfx_entr_detr_tendency!(Yₜ, Y, p, t, colidx, turbconv_model::EDMFX) return nothing 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, - ::ConstantTimescaleDetrainment, -) where {FT} - detr_tau = CAP.detr_tau(params) - detr = min(1 / detr_tau, 1 / dt) - return detr +function edmfx_entr_detr_tendency!( + Yₜ, + Y, + p, + t, + colidx, + turbconv_model::AdvectiveEDMFX, +) + + n = n_mass_flux_subdomains(turbconv_model) + (; ᶜentrʲs, ᶜdetrʲs) = p + (; ᶜq_tot⁰, ᶜh_tot⁰, ᶠu₃⁰) = p + + for j in 1:n + + @. Yₜ.c.sgsʲs.:($$j).ρa[colidx] += + Y.c.sgsʲs.:($$j).ρa[colidx] * + (ᶜentrʲs.:($$j)[colidx] - ᶜdetrʲs.:($$j)[colidx]) + + @. Yₜ.c.sgsʲs.:($$j).h_tot[colidx] += + ᶜentrʲs.:($$j)[colidx] * + (ᶜh_tot⁰[colidx] - Y.c.sgsʲs.:($$j).h_tot[colidx]) + + @. Yₜ.c.sgsʲs.:($$j).q_tot[colidx] += + ᶜentrʲs.:($$j)[colidx] * + (ᶜq_tot⁰[colidx] - Y.c.sgsʲs.:($$j).q_tot[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 diff --git a/src/prognostic_equations/edmfx_sgs_flux.jl b/src/prognostic_equations/edmfx_sgs_flux.jl index 1034f36244..c1c1d7cc82 100644 --- a/src/prognostic_equations/edmfx_sgs_flux.jl +++ b/src/prognostic_equations/edmfx_sgs_flux.jl @@ -10,19 +10,18 @@ function edmfx_sgs_mass_flux_tendency!( p, t, colidx, - turbconv_model::EDMFX, + turbconv_model::Union{EDMFX, AdvectiveEDMFX}, ) FT = Spaces.undertype(axes(Y.c)) n = n_mass_flux_subdomains(turbconv_model) - (; edmfx_upwinding, sfc_conditions) = p + (; edmfx_upwinding) = p (; ᶠu³, ᶜh_tot, ᶜspecific) = p - (; ᶠu³ʲs, ᶜh_totʲs, ᶜspecificʲs) = p - (; ᶜρa⁰, ᶠu³⁰, ᶜu⁰, ᶜh_tot⁰, ᶜspecific⁰) = p - (; ᶜK_u, ᶜK_h) = p + (; ᶠu³ʲs) = p + (; ᶜρa⁰, ᶠu³⁰, ᶜu⁰, ᶜh_tot⁰) = p + ᶜq_tot⁰ = turbconv_model isa EDMFX ? p.ᶜspecific⁰.q_tot : p.ᶜq_tot⁰ (; dt) = p.simulation ᶜJ = Fields.local_geometry_field(Y.c).J - ᶠgradᵥ = Operators.GradientC2F() if p.atmos.edmfx_sgs_mass_flux # energy @@ -30,7 +29,10 @@ function edmfx_sgs_mass_flux_tendency!( ᶜh_tot_diff_colidx = ᶜq_tot_diff_colidx = p.ᶜtemp_scalar[colidx] for j in 1:n @. ᶠu³_diff_colidx = ᶠu³ʲs.:($$j)[colidx] - ᶠu³[colidx] - @. ᶜh_tot_diff_colidx = ᶜh_totʲs.:($$j)[colidx] - ᶜh_tot[colidx] + ᶜh_totʲ_colidx = + turbconv_model isa EDMFX ? p.ᶜh_totʲs.:($j)[colidx] : + Y.c.sgsʲs.:($j).h_tot[colidx] + @. ᶜh_tot_diff_colidx = ᶜh_totʲ_colidx - ᶜh_tot[colidx] vertical_transport!( Yₜ.c.ρe_tot[colidx], ᶜJ[colidx], @@ -57,8 +59,11 @@ function edmfx_sgs_mass_flux_tendency!( # specific humidity for j in 1:n @. ᶠu³_diff_colidx = ᶠu³ʲs.:($$j)[colidx] - ᶠu³[colidx] - @. ᶜq_tot_diff_colidx = - ᶜspecificʲs.:($$j).q_tot[colidx] - ᶜspecific.q_tot[colidx] + ᶜq_totʲ_colidx = + turbconv_model isa EDMFX ? + p.ᶜspecificʲs.:($j).q_tot[colidx] : + Y.c.sgsʲs.:($j).q_tot[colidx] + @. ᶜq_tot_diff_colidx = ᶜq_totʲ_colidx - ᶜspecific.q_tot[colidx] vertical_transport!( Yₜ.c.ρq_tot[colidx], ᶜJ[colidx], @@ -70,8 +75,7 @@ function edmfx_sgs_mass_flux_tendency!( ) end @. ᶠu³_diff_colidx = ᶠu³⁰[colidx] - ᶠu³[colidx] - @. ᶜq_tot_diff_colidx = - ᶜspecific⁰.q_tot[colidx] - ᶜspecific.q_tot[colidx] + @. ᶜq_tot_diff_colidx = ᶜq_tot⁰[colidx] - ᶜspecific.q_tot[colidx] vertical_transport!( Yₜ.c.ρq_tot[colidx], ᶜJ[colidx], @@ -159,18 +163,14 @@ function edmfx_sgs_diffusive_flux_tendency!( p, t, colidx, - turbconv_model::EDMFX, + turbconv_model::Union{EDMFX, AdvectiveEDMFX}, ) FT = Spaces.undertype(axes(Y.c)) - n = n_mass_flux_subdomains(turbconv_model) - (; edmfx_upwinding, sfc_conditions) = p - (; ᶠu³, ᶜh_tot, ᶜspecific) = p - (; ᶠu³ʲs, ᶜh_totʲs, ᶜspecificʲs) = p - (; ᶜρa⁰, ᶠu³⁰, ᶜu⁰, ᶜh_tot⁰, ᶜspecific⁰) = p + (; sfc_conditions) = p + (; ᶜρa⁰, ᶠu³⁰, ᶜu⁰, ᶜh_tot⁰) = p + ᶜq_tot⁰ = turbconv_model isa EDMFX ? p.ᶜspecific⁰.q_tot : p.ᶜq_tot⁰ (; ᶜK_u, ᶜK_h) = p - (; dt) = p.simulation - ᶜJ = Fields.local_geometry_field(Y.c).J ᶠgradᵥ = Operators.GradientC2F() if p.atmos.edmfx_sgs_diffusive_flux @@ -193,9 +193,8 @@ function edmfx_sgs_diffusive_flux_tendency!( sfc_conditions.ρ_flux_q_tot[colidx], ), ) - @. Yₜ.c.ρq_tot[colidx] -= ᶜdivᵥ_ρq_tot( - -(ᶠρaK_h[colidx] * ᶠgradᵥ(ᶜspecific⁰.q_tot[colidx])), - ) + @. Yₜ.c.ρq_tot[colidx] -= + ᶜdivᵥ_ρq_tot(-(ᶠρaK_h[colidx] * ᶠgradᵥ(ᶜq_tot⁰[colidx]))) end # momentum diff --git a/src/prognostic_equations/hyperdiffusion.jl b/src/prognostic_equations/hyperdiffusion.jl index c7a886184d..f81c58fb66 100644 --- a/src/prognostic_equations/hyperdiffusion.jl +++ b/src/prognostic_equations/hyperdiffusion.jl @@ -21,7 +21,10 @@ function hyperdiffusion_cache(Y, atmos, do_dss) # Sub-grid scale quantities ᶜ∇²uʲs = - atmos.turbconv_model isa EDMFX ? similar(Y.c, NTuple{n, C123{FT}}) : (;) + ( + atmos.turbconv_model isa EDMFX || + atmos.turbconv_model isa AdvectiveEDMFX + ) ? similar(Y.c, NTuple{n, C123{FT}}) : (;) sgs_quantities = atmos.turbconv_model isa EDMFX ? (; @@ -33,6 +36,14 @@ function hyperdiffusion_cache(Y, atmos, do_dss) specific_sgsʲs.(Y.c, atmos.turbconv_model) ), ) : + atmos.turbconv_model isa AdvectiveEDMFX ? + (; + ᶜ∇²tke⁰ = similar(Y.c, FT), + ᶜ∇²uₕʲs = similar(Y.c, NTuple{n, C12{FT}}), + ᶜ∇²uᵥʲs = similar(Y.c, NTuple{n, C3{FT}}), + ᶜ∇²h_totʲs = similar(Y.c, NTuple{n, FT}), + ᶜ∇²q_totʲs = similar(Y.c, NTuple{n, FT}), + ) : atmos.turbconv_model isa DiagnosticEDMFX ? (; ᶜ∇²tke⁰ = similar(Y.c, FT)) : (;) quantities = (; gs_quantities..., sgs_quantities...) @@ -71,6 +82,9 @@ NVTX.@annotate function hyperdiffusion_tendency!(Yₜ, Y, p, t) ᶜ∇²specific_energyʲs, ) = p end + if turbconv_model isa AdvectiveEDMFX + (; ᶜρa⁰, ᶜρʲs, ᶜ∇²tke⁰, ᶜtke⁰, ᶜ∇²uₕʲs, ᶜ∇²uᵥʲs, ᶜ∇²uʲs, ᶜ∇²h_totʲs) = p + end if turbconv_model isa DiagnosticEDMFX (; ᶜtke⁰, ᶜ∇²tke⁰) = p end @@ -89,7 +103,9 @@ NVTX.@annotate function hyperdiffusion_tendency!(Yₜ, Y, p, t) if turbconv_model isa EDMFX && diffuse_tke @. ᶜ∇²tke⁰ = wdivₕ(gradₕ(ᶜspecific⁰.tke)) end - if turbconv_model isa DiagnosticEDMFX && diffuse_tke + if ( + turbconv_model isa DiagnosticEDMFX || turbconv_model isa AdvectiveEDMFX + ) && diffuse_tke @. ᶜ∇²tke⁰ = wdivₕ(gradₕ(ᶜtke⁰)) end @@ -110,6 +126,15 @@ NVTX.@annotate function hyperdiffusion_tendency!(Yₜ, Y, p, t) end end + if turbconv_model isa AdvectiveEDMFX + for j in 1:n + @. ᶜ∇²uʲs.:($$j) = + C123(wgradₕ(divₕ(p.ᶜuʲs.:($$j)))) - + C123(wcurlₕ(C123(curlₕ(p.ᶜuʲs.:($$j))))) + @. ᶜ∇²h_totʲs.:($$j) = wdivₕ(gradₕ(Y.c.sgsʲs.:($$j).h_tot)) + end + end + if do_dss NVTX.@range "dss_hyperdiffusion_tendency" color = colorant"green" begin for dss_op! in ( @@ -122,10 +147,11 @@ NVTX.@annotate function hyperdiffusion_tendency!(Yₜ, Y, p, t) # operations do not accept Covariant123Vector types dss_op!(ᶜ∇²u, buffer.ᶜ∇²u) dss_op!(ᶜ∇²specific_energy, buffer.ᶜ∇²specific_energy) - if turbconv_model isa EDMFX && diffuse_tke + if diffuse_tke dss_op!(ᶜ∇²tke⁰, buffer.ᶜ∇²tke⁰) end - if turbconv_model isa EDMFX + if (turbconv_model isa EDMFX) || + (turbconv_model isa AdvectiveEDMFX) # Need to split the DSS computation here, because our DSS # operations do not accept Covariant123Vector types for j in 1:n @@ -138,10 +164,12 @@ NVTX.@annotate function hyperdiffusion_tendency!(Yₜ, Y, p, t) @. ᶜ∇²uʲs.:($$j) = C123(ᶜ∇²uₕʲs.:($$j)) + C123(ᶜ∇²uᵥʲs.:($$j)) end + end + if turbconv_model isa EDMFX dss_op!(ᶜ∇²specific_energyʲs, buffer.ᶜ∇²specific_energyʲs) end - if turbconv_model isa DiagnosticEDMFX && diffuse_tke - dss_op!(ᶜ∇²tke⁰, buffer.ᶜ∇²tke⁰) + if turbconv_model isa AdvectiveEDMFX + dss_op!(ᶜ∇²h_totʲs, buffer.ᶜ∇²h_totʲs) end end end @@ -158,7 +186,8 @@ NVTX.@annotate function hyperdiffusion_tendency!(Yₜ, Y, p, t) @. ᶜρ_energyₜ -= κ₄ * wdivₕ(Y.c.ρ * gradₕ(ᶜ∇²specific_energy)) # Sub-grid scale hyperdiffusion continued - if turbconv_model isa EDMFX && diffuse_tke + if (turbconv_model isa EDMFX || turbconv_model isa AdvectiveEDMFX) && + diffuse_tke @. Yₜ.c.sgs⁰.ρatke -= κ₄ * wdivₕ(ᶜρa⁰ * gradₕ(ᶜ∇²tke⁰)) end if turbconv_model isa EDMFX @@ -177,6 +206,19 @@ NVTX.@annotate function hyperdiffusion_tendency!(Yₜ, Y, p, t) wdivₕ(Y.c.sgsʲs.:($$j).ρa * gradₕ(ᶜ∇²specific_energyʲs.:($$j))) end end + if turbconv_model isa AdvectiveEDMFX + for j in 1:n + if point_type <: Geometry.Abstract3DPoint + # only need curl-curl part + @. ᶜ∇²uᵥʲs.:($$j) = C3(wcurlₕ(C123(curlₕ(ᶜ∇²uʲs.:($$j))))) + @. Yₜ.f.sgsʲs.:($$j).u₃ += + κ₄ * ᶠwinterp(ᶜJ * Y.c.ρ, ᶜ∇²uᵥʲs.:($$j)) + end + # Note: It is more correct to have ρa inside and outside the divergence + @. Yₜ.c.sgsʲs.:($$j).h_tot -= κ₄ * wdivₕ(gradₕ(ᶜ∇²h_totʲs.:($$j))) + end + end + if turbconv_model isa DiagnosticEDMFX && diffuse_tke @. Yₜ.c.sgs⁰.ρatke -= κ₄ * wdivₕ(Y.c.ρ * gradₕ(ᶜ∇²tke⁰)) end @@ -193,6 +235,9 @@ NVTX.@annotate function tracer_hyperdiffusion_tendency!(Yₜ, Y, p, t) if turbconv_model isa EDMFX (; ᶜspecificʲs, ᶜ∇²specific_tracersʲs, ᶜp) = p end + if turbconv_model isa AdvectiveEDMFX + (; ᶜ∇²q_totʲs) = p + end if do_dss buffer = p.hyperdiffusion_ghost_buffer end @@ -209,6 +254,13 @@ NVTX.@annotate function tracer_hyperdiffusion_tendency!(Yₜ, Y, p, t) end end + if turbconv_model isa AdvectiveEDMFX + for j in 1:n + # Note: It is more correct to have ρa inside and outside the divergence + @. ᶜ∇²q_totʲs.:($$j) = wdivₕ(gradₕ(Y.c.sgsʲs.:($$j).q_tot)) + end + end + if do_dss NVTX.@range "dss_hyperdiffusion_tendency" color = colorant"green" begin for dss_op! in ( @@ -220,6 +272,9 @@ NVTX.@annotate function tracer_hyperdiffusion_tendency!(Yₜ, Y, p, t) if turbconv_model isa EDMFX dss_op!(ᶜ∇²specific_tracersʲs, buffer.ᶜ∇²specific_tracersʲs) end + if turbconv_model isa AdvectiveEDMFX + dss_op!(ᶜ∇²q_totʲs, buffer.ᶜ∇²q_totʲs) + end end end end @@ -245,5 +300,12 @@ NVTX.@annotate function tracer_hyperdiffusion_tendency!(Yₜ, Y, p, t) end end end + if turbconv_model isa AdvectiveEDMFX + for j in 1:n + @. Yₜ.c.sgsʲs.:($$j).ρa -= + κ₄ * wdivₕ(Y.c.sgsʲs.:($$j).ρa * gradₕ(ᶜ∇²q_totʲs.:($$j))) + @. Yₜ.c.sgsʲs.:($$j).q_tot -= κ₄ * wdivₕ(gradₕ(ᶜ∇²q_totʲs.:($$j))) + end + end return nothing end diff --git a/src/prognostic_equations/implicit/implicit_tendency.jl b/src/prognostic_equations/implicit/implicit_tendency.jl index c7c97a68f5..ca89f3582b 100644 --- a/src/prognostic_equations/implicit/implicit_tendency.jl +++ b/src/prognostic_equations/implicit/implicit_tendency.jl @@ -52,6 +52,13 @@ vertical_transport!(ᶜρχₜ, ᶜJ, ᶜρ, ᶠu³, ᶜχ, dt, ::Val{:zalesak}) ), )) +vertical_advection!(ᶜρχₜ, ᶠu³, ᶜχ, ::Val{:none}) = + @. ᶜρχₜ -= ᶜadvdivᵥ(ᶠu³ * ᶠinterp(ᶜχ)) - ᶜχ * ᶜadvdivᵥ(ᶠu³) +vertical_advection!(ᶜρχₜ, ᶠu³, ᶜχ, ::Val{:first_order}) = + @. ᶜρχₜ -= ᶜadvdivᵥ(ᶠupwind1(ᶠu³, ᶜχ)) - ᶜχ * ᶜadvdivᵥ(ᶠu³) +vertical_advection!(ᶜρχₜ, ᶠu³, ᶜχ, ::Val{:third_order}) = + @. ᶜρχₜ -= ᶜadvdivᵥ(ᶠupwind3(ᶠu³, ᶜχ)) - ᶜχ * ᶜadvdivᵥ(ᶠu³) + function implicit_vertical_advection_tendency!(Yₜ, Y, p, t, colidx) (; energy_upwinding, tracer_upwinding, density_upwinding, edmfx_upwinding) = p @@ -107,7 +114,7 @@ function implicit_vertical_advection_tendency!(Yₜ, Y, p, t, colidx) if rayleigh_sponge isa RayleighSponge (; ᶠβ_rayleigh_w) = p @. Yₜ.f.u₃[colidx] -= ᶠβ_rayleigh_w[colidx] * Y.f.u₃[colidx] - if turbconv_model isa EDMFX + if turbconv_model isa EDMFX || turbconv_model isa AdvectiveEDMFX for j in 1:n @. Yₜ.f.sgsʲs.:($$j).u₃[colidx] -= ᶠβ_rayleigh_w[colidx] * Y.f.sgsʲs.:($$j).u₃[colidx] diff --git a/src/prognostic_equations/pressure_work.jl b/src/prognostic_equations/pressure_work.jl index a190c2cc5f..6df349c1e0 100644 --- a/src/prognostic_equations/pressure_work.jl +++ b/src/prognostic_equations/pressure_work.jl @@ -14,4 +14,7 @@ function pressure_work_tendency!(Yₜ, Y, p, t, colidx, ::EDMFX) return nothing end +# We ignore the small pressure term in advective EDMF for now +pressure_work_tendency!(Yₜ, Y, p, t, colidx, ::AdvectiveEDMFX) = nothing + pressure_work_tendency!(Yₜ, Y, p, t, colidx, ::Any) = nothing diff --git a/src/prognostic_equations/zero_velocity.jl b/src/prognostic_equations/zero_velocity.jl index f62a79b55c..cb6787ddeb 100644 --- a/src/prognostic_equations/zero_velocity.jl +++ b/src/prognostic_equations/zero_velocity.jl @@ -10,7 +10,8 @@ function zero_velocity_tendency!(Yₜ, Y, p, t, colidx) @. Yₜ.c.uₕ[colidx] = C12(FT(0), FT(0)) @. Yₜ.f.u₃[colidx] = Geometry.Covariant3Vector(FT(0)) - if p.atmos.turbconv_model isa EDMFX + if p.atmos.turbconv_model isa EDMFX || + p.atmos.turbconv_model isa AdvectiveEDMFX for j in 1:n @. Yₜ.f.sgsʲs.:($$j).u₃[colidx] = Geometry.Covariant3Vector(FT(0)) diff --git a/src/solver/model_getters.jl b/src/solver/model_getters.jl index ca0600dfb2..e1c0d7f35e 100644 --- a/src/solver/model_getters.jl +++ b/src/solver/model_getters.jl @@ -341,7 +341,8 @@ function get_turbconv_model( turbconv_params, ) turbconv = parsed_args["turbconv"] - @assert turbconv in (nothing, "edmf", "edmfx", "diagnostic_edmfx") + @assert turbconv in + (nothing, "edmf", "edmfx", "advective_edmfx", "diagnostic_edmfx") return if turbconv == "edmf" TC.EDMFModel( @@ -355,6 +356,10 @@ function get_turbconv_model( N = turbconv_params.updraft_number TKE = parsed_args["prognostic_tke"] EDMFX{N, TKE}(turbconv_params.min_area) + elseif turbconv == "advective_edmfx" + N = turbconv_params.updraft_number + TKE = parsed_args["prognostic_tke"] + AdvectiveEDMFX{N, TKE}(turbconv_params.min_area) elseif turbconv == "diagnostic_edmfx" N = turbconv_params.updraft_number TKE = parsed_args["prognostic_tke"] @@ -387,14 +392,14 @@ function get_detrainment_model(parsed_args) NoDetrainment() elseif detr_model == "PiGroups" PiGroupsDetrainment() - elseif detr_model == "ConstantCoefficient" - ConstantCoefficientDetrainment() + elseif detr_model == "BOverW" + BOverWDetrainment() elseif detr_model == "ConstantCoefficientHarmonics" ConstantCoefficientHarmonicsDetrainment() elseif detr_model == "ConstantTimescale" ConstantTimescaleDetrainment() else - error("Invalid entr_model $(entr_model)") + error("Invalid detr_model $(detr_model)") end end diff --git a/src/solver/types.jl b/src/solver/types.jl index 9673374188..0d4ee59794 100644 --- a/src/solver/types.jl +++ b/src/solver/types.jl @@ -157,6 +157,12 @@ struct EDMFX{N, TKE, FT} <: AbstractEDMF end EDMFX{N, TKE}(a_half::FT) where {N, TKE, FT} = EDMFX{N, TKE, FT}(a_half) +struct AdvectiveEDMFX{N, TKE, FT} <: AbstractEDMF + a_half::FT # WARNING: this should never be used outside of divide_by_ρa +end +AdvectiveEDMFX{N, TKE}(a_half::FT) where {N, TKE, FT} = + AdvectiveEDMFX{N, TKE, FT}(a_half) + struct DiagnosticEDMFX{N, TKE, FT} <: AbstractEDMF a_int::FT # area fraction of the first interior cell above the surface a_half::FT # WARNING: this should never be used outside of divide_by_ρa @@ -165,13 +171,16 @@ DiagnosticEDMFX{N, TKE}(a_int::FT, a_half::FT) where {N, TKE, FT} = DiagnosticEDMFX{N, TKE, FT}(a_int, a_half) n_mass_flux_subdomains(::EDMFX{N}) where {N} = N +n_mass_flux_subdomains(::AdvectiveEDMFX{N}) where {N} = N n_mass_flux_subdomains(::DiagnosticEDMFX{N}) where {N} = N n_mass_flux_subdomains(::Any) = 0 n_prognostic_mass_flux_subdomains(::EDMFX{N}) where {N} = N +n_prognostic_mass_flux_subdomains(::AdvectiveEDMFX{N}) where {N} = N n_prognostic_mass_flux_subdomains(::Any) = 0 use_prognostic_tke(::EDMFX{N, TKE}) where {N, TKE} = TKE +use_prognostic_tke(::AdvectiveEDMFX{N, TKE}) where {N, TKE} = TKE use_prognostic_tke(::DiagnosticEDMFX{N, TKE}) where {N, TKE} = TKE use_prognostic_tke(::Any) = false @@ -186,7 +195,7 @@ abstract type AbstractDetrainmentModel end struct NoDetrainment <: AbstractDetrainmentModel end struct PiGroupsDetrainment <: AbstractDetrainmentModel end -struct ConstantCoefficientDetrainment <: AbstractDetrainmentModel end +struct BOverWDetrainment <: AbstractDetrainmentModel end struct ConstantCoefficientHarmonicsDetrainment <: AbstractDetrainmentModel end struct ConstantTimescaleDetrainment <: AbstractDetrainmentModel end @@ -234,6 +243,7 @@ Base.broadcastable(x::AbstractEnergyFormulation) = tuple(x) Base.broadcastable(x::AbstractPrecipitationModel) = tuple(x) Base.broadcastable(x::AbstractForcing) = tuple(x) Base.broadcastable(x::EDMFX) = tuple(x) +Base.broadcastable(x::AdvectiveEDMFX) = tuple(x) Base.broadcastable(x::DiagnosticEDMFX) = tuple(x) Base.broadcastable(x::AbstractEntrainmentModel) = tuple(x) Base.broadcastable(x::AbstractDetrainmentModel) = tuple(x) diff --git a/src/utils/variable_manipulations.jl b/src/utils/variable_manipulations.jl index 03fe35c6d7..e72ea7051f 100644 --- a/src/utils/variable_manipulations.jl +++ b/src/utils/variable_manipulations.jl @@ -260,6 +260,22 @@ mass-flux subdomain states are stored in `gs.sgsʲs`. """ ρa⁺(gs) = mapreduce(sgsʲ -> sgsʲ.ρa, +, gs.sgsʲs) +""" + ρah_tot⁺(gs) + +Computes the total mass-flux subdomain area-weighted ρh_tot, assuming that the +mass-flux subdomain states are stored in `sgsʲs`. +""" +ρah_tot⁺(sgsʲs) = mapreduce(sgsʲ -> sgsʲ.ρa * sgsʲ.h_tot, +, sgsʲs) + +""" + ρaq_tot⁺(gs) + +Computes the total mass-flux subdomain area-weighted ρq_tot, assuming that the +mass-flux subdomain states are stored in `sgsʲs`. +""" +ρaq_tot⁺(sgsʲs) = mapreduce(sgsʲ -> sgsʲ.ρa * sgsʲ.q_tot, +, sgsʲs) + """ ρa⁰(gs) diff --git a/toml/advective_edmfx_box.toml b/toml/advective_edmfx_box.toml new file mode 100644 index 0000000000..3c8a8e76d5 --- /dev/null +++ b/toml/advective_edmfx_box.toml @@ -0,0 +1,24 @@ +[C_E] +alias = "C_E" +value = 0.044 +type = "float" + +[EDMF_surface_area] +alias = "surface_area" +value = 0.1 +type = "float" + +[EDMF_min_area] +alias = "min_area" +value = 1.0e-3 +type = "float" + +[entr_coeff] +alias = "entr_coeff" +value = 0.1 +type = "float" + +[entr_tau] +alias = "entr_tau" +value = 1800 +type = "float"