From e4f8275afb6559c94a103fb483749a3db800afec Mon Sep 17 00:00:00 2001 From: Zhaoyi Shen <11598433+szy21@users.noreply.github.com> Date: Mon, 13 Nov 2023 10:59:38 -0800 Subject: [PATCH] remove theta as an energy variable --- .buildkite/pipeline.yml | 13 --- config/default_configs/default_config.yml | 3 - .../sphere_held_suarez_rhotheta.yml | 6 -- post_processing/post_processing_funcs.jl | 7 +- post_processing/remap/remap_helpers.jl | 17 +--- src/cache/precomputed_quantities.jl | 62 +++++------- .../prognostic_edmf_precomputed_quantities.jl | 3 +- src/callbacks/callbacks.jl | 41 ++++---- src/diagnostics/Diagnostics.jl | 3 - src/diagnostics/core_diagnostics.jl | 41 ++------ src/diagnostics/default_diagnostics.jl | 12 +-- src/initial_conditions/InitialConditions.jl | 2 - src/initial_conditions/atmos_state.jl | 7 +- .../radiation/held_suarez.jl | 7 +- .../radiation/radiation.jl | 7 +- .../sponge/viscous_sponge.jl | 17 +--- src/prognostic_equations/advection.jl | 8 +- src/prognostic_equations/hyperdiffusion.jl | 10 +- .../implicit/implicit_solver.jl | 99 +++++-------------- .../implicit/implicit_tendency.jl | 2 +- .../vertical_diffusion_boundary_layer.jl | 2 - src/solver/model_getters.jl | 13 --- src/solver/type_getters.jl | 1 - src/solver/types.jl | 7 -- src/surface_conditions/SurfaceConditions.jl | 2 - src/surface_conditions/surface_conditions.jl | 18 +--- src/utils/discrete_hydrostatic_balance.jl | 1 - src/utils/utilities.jl | 2 +- src/utils/variable_manipulations.jl | 4 +- 29 files changed, 99 insertions(+), 318 deletions(-) delete mode 100644 config/model_configs/sphere_held_suarez_rhotheta.yml diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml index 7e3603df2f..a276337816 100644 --- a/.buildkite/pipeline.yml +++ b/.buildkite/pipeline.yml @@ -279,19 +279,6 @@ steps: command: "julia --color=yes --project=examples examples/hybrid/driver.jl --config_file $CONFIG_PATH/sphere_baroclinic_wave_rhoe_equilmoist_conservation_ft64.yml" artifact_paths: "sphere_baroclinic_wave_rhoe_equilmoist_conservation_ft64/*" - - label: ":computer: held suarez (ρθ)" - command: > - julia --color=yes --project=examples examples/hybrid/driver.jl - --config_file $CONFIG_PATH/sphere_held_suarez_rhotheta.yml - - julia --color=yes --project=examples post_processing/remap/remap_pipeline.jl - --data_dir sphere_held_suarez_rhotheta --out_dir sphere_held_suarez_rhotheta - - julia --color=yes --project=examples post_processing/plot/plot_pipeline.jl - --nc_dir sphere_held_suarez_rhotheta --fig_dir sphere_held_suarez_rhotheta - --case_name dry_held_suarez - artifact_paths: "sphere_held_suarez_rhotheta/*" - - label: ":computer: held suarez (ρe) hightop" key: sphere_held_suarez_rhoe_hightop command: > diff --git a/config/default_configs/default_config.yml b/config/default_configs/default_config.yml index c6f5886a17..eec06174f0 100644 --- a/config/default_configs/default_config.yml +++ b/config/default_configs/default_config.yml @@ -138,9 +138,6 @@ forcing: test_dycore_consistency: help: "Test dycore consistency [`false` (default), `true`]" value: false -energy_name: - help: "Energy variable name [`rhoe` (default), `rhotheta`]" - value: "rhoe" dt_save_to_disk: help: "Time between saving to disk. Examples: [`10secs`, `1hours`, `Inf` (do not save)]" value: "Inf" diff --git a/config/model_configs/sphere_held_suarez_rhotheta.yml b/config/model_configs/sphere_held_suarez_rhotheta.yml deleted file mode 100644 index fab665b1fd..0000000000 --- a/config/model_configs/sphere_held_suarez_rhotheta.yml +++ /dev/null @@ -1,6 +0,0 @@ -energy_name: "rhotheta" -dt_save_to_disk: "10days" -t_end: "20days" -forcing: "held_suarez" -dt: "400secs" -job_id: "sphere_held_suarez_rhotheta" diff --git a/post_processing/post_processing_funcs.jl b/post_processing/post_processing_funcs.jl index 1b1eae488b..81811500ed 100644 --- a/post_processing/post_processing_funcs.jl +++ b/post_processing/post_processing_funcs.jl @@ -7,12 +7,7 @@ import ClimaCoreTempestRemap: def_space_coord import ClimaCoreSpectra: power_spectrum_1d, power_spectrum_2d import ClimaCore: Geometry, Fields, Spaces using LinearAlgebra: norm -import ClimaAtmos: - DryModel, - EquilMoistModel, - NonEquilMoistModel, - PotentialTemperature, - TotalEnergy +import ClimaAtmos: DryModel, EquilMoistModel, NonEquilMoistModel function process_name(s::AbstractString) # "c_ρ", "c_ρe", "c_uₕ_1", "c_uₕ_2", "f_w_1" diff --git a/post_processing/remap/remap_helpers.jl b/post_processing/remap/remap_helpers.jl index 2a17dd6b73..a1ca3fcc89 100644 --- a/post_processing/remap/remap_helpers.jl +++ b/post_processing/remap/remap_helpers.jl @@ -110,14 +110,7 @@ function remap2latlon(filein, data_dir, remap_tmpdir, weightfile, nlat, nlon) nc_time = def_time_coord(nc) # define variables for the prognostic states nc_rho = defVar(nc, "rho", FT, cspace, ("time",)) - thermo_var = if :ρe_tot in propertynames(Y.c) - "e_tot" - elseif :ρθ in propertynames(Y.c) - "theta" - else - error("Unfound thermodynamic variable") - end - nc_thermo = defVar(nc, thermo_var, FT, cspace, ("time",)) + nc_thermo = defVar(nc, "e_tot", FT, cspace, ("time",)) nc_u = defVar(nc, "u", FT, cspace, ("time",)) nc_v = defVar(nc, "v", FT, cspace, ("time",)) nc_w = defVar(nc, "w", FT, cspace, ("time",)) @@ -199,11 +192,7 @@ function remap2latlon(filein, data_dir, remap_tmpdir, weightfile, nlat, nlon) # density nc_rho[:, 1] = Y.c.ρ # thermodynamics - if :ρe_tot in propertynames(Y.c) - nc_thermo[:, 1] = Y.c.ρe_tot ./ Y.c.ρ - elseif :ρθ in propertynames(Y.c) - nc_thermo[:, 1] = Y.c.ρθ ./ Y.c.ρ - end + nc_thermo[:, 1] = Y.c.ρe_tot ./ Y.c.ρ # physical horizontal velocity uh_phy = Geometry.transform.(tuple(Geometry.UVAxis()), Y.c.uₕ) nc_u[:, 1] = uh_phy.components.data.:1 @@ -279,7 +268,7 @@ function remap2latlon(filein, data_dir, remap_tmpdir, weightfile, nlat, nlon) joinpath(out_dir, first(splitext(basename(filein))) * ".nc") dry_variables = [ "rho", - thermo_var, + "e_tot", "u", "v", "w", diff --git a/src/cache/precomputed_quantities.jl b/src/cache/precomputed_quantities.jl index 116e60fdf5..b0f9a5c9e3 100644 --- a/src/cache/precomputed_quantities.jl +++ b/src/cache/precomputed_quantities.jl @@ -16,8 +16,6 @@ Allocates and returns the precomputed quantities: - `ᶜts`: the thermodynamic state on cell centers - `ᶜp`: the air pressure on cell centers - `sfc_conditions`: the conditions at the surface (at the bottom cell faces) - -If the `energy_form` is TotalEnergy, there is an additional quantity: - `ᶜh_tot`: the total enthalpy on cell centers If the `turbconv_model` is EDMFX, there also two SGS versions of every quantity @@ -36,14 +34,10 @@ TODO: Rename `ᶜK` to `ᶜκ`. """ function precomputed_quantities(Y, atmos) FT = eltype(Y) - @assert ( - !(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 PrognosticEDMFX) + @assert !(atmos.moisture_model isa DryModel) || + !(atmos.turbconv_model isa DiagnosticEDMFX) + @assert !(atmos.moisture_model isa DryModel) || + !(atmos.turbconv_model isa PrognosticEDMFX) @assert !(atmos.edmfx_detr_model isa ConstantAreaDetrainment) || !(atmos.turbconv_model isa DiagnosticEDMFX) TST = thermo_state_type(atmos.moisture_model, FT) @@ -56,10 +50,7 @@ function precomputed_quantities(Y, atmos) ᶜK = similar(Y.c, FT), ᶜts = similar(Y.c, TST), ᶜp = similar(Y.c, FT), - ( - atmos.energy_form isa TotalEnergy ? - (; ᶜh_tot = similar(Y.c, FT)) : (;) - )..., + ᶜh_tot = similar(Y.c, FT), sfc_conditions = Fields.Field(SCT, Spaces.level(axes(Y.f), half)), ) advective_sgs_quantities = @@ -252,12 +243,8 @@ function thermo_state( return get_ts(ρ, p, θ, e_int, q_tot, q_pt) end -function thermo_vars(energy_form, moisture_model, specific, K, Φ) - energy_var = if energy_form isa PotentialTemperature - (; specific.θ) - elseif energy_form isa TotalEnergy - (; e_int = specific.e_tot - K - Φ) - end +function thermo_vars(moisture_model, specific, K, Φ) + energy_var = (; e_int = specific.e_tot - K - Φ) moisture_var = if moisture_model isa DryModel (;) elseif moisture_model isa EquilMoistModel @@ -269,19 +256,17 @@ function thermo_vars(energy_form, moisture_model, specific, K, Φ) return (; energy_var..., moisture_var...) end -ts_gs(thermo_params, energy_form, moisture_model, specific, K, Φ, ρ) = - thermo_state( - thermo_params; - thermo_vars(energy_form, moisture_model, specific, K, Φ)..., - ρ, - ) +ts_gs(thermo_params, moisture_model, specific, K, Φ, ρ) = thermo_state( + thermo_params; + thermo_vars(moisture_model, specific, K, Φ)..., + ρ, +) -ts_sgs(thermo_params, energy_form, moisture_model, specific, K, Φ, p) = - thermo_state( - thermo_params; - thermo_vars(energy_form, moisture_model, specific, K, Φ)..., - p, - ) +ts_sgs(thermo_params, moisture_model, specific, K, Φ, p) = thermo_state( + thermo_params; + thermo_vars(moisture_model, specific, K, Φ)..., + p, +) """ set_precomputed_quantities!(Y, p, t) @@ -300,10 +285,10 @@ function instead of recomputing the value yourself. Otherwise, it will be difficult to ensure that the duplicated computations are consistent. """ NVTX.@annotate function set_precomputed_quantities!(Y, p, t) - (; energy_form, moisture_model, turbconv_model) = p.atmos + (; moisture_model, turbconv_model) = p.atmos thermo_params = CAP.thermodynamics_params(p.params) n = n_mass_flux_subdomains(turbconv_model) - thermo_args = (thermo_params, energy_form, moisture_model) + thermo_args = (thermo_params, moisture_model) (; ᶜΦ) = p.core (; ᶜspecific, ᶜu, ᶠu³, ᶜK, ᶜts, ᶜp) = p.precomputed ᶠuₕ³ = p.scratch.ᶠtemp_CT3 @@ -323,7 +308,7 @@ NVTX.@annotate function set_precomputed_quantities!(Y, p, t) # quantities of the form ᶜρaχ⁰ / ᶜρ and ᶜρaχʲ / ᶜρ. However, we cannot # compute ᶜρ⁰ and ᶜρʲ without first computing ᶜts⁰ and ᶜtsʲ, both of # which depend on the value of ᶜp, which in turn depends on ᶜts. Since - # ᶜts depends on ᶜK (at least when the energy_form is TotalEnergy), this + # ᶜts depends on ᶜK, this # means that the amount by which ᶜK needs to be incremented is a # function of ᶜK itself. So, unless we run a nonlinear solver here, this # circular dependency will prevent us from computing the exact value of @@ -335,11 +320,8 @@ NVTX.@annotate function set_precomputed_quantities!(Y, p, t) @. ᶜts = ts_gs(thermo_args..., ᶜspecific, ᶜK, ᶜΦ, Y.c.ρ) @. ᶜp = TD.air_pressure(thermo_params, ᶜts) - if energy_form isa TotalEnergy - (; ᶜh_tot) = p.precomputed - @. ᶜh_tot = - TD.total_specific_enthalpy(thermo_params, ᶜts, ᶜspecific.e_tot) - end + (; ᶜh_tot) = p.precomputed + @. ᶜh_tot = TD.total_specific_enthalpy(thermo_params, ᶜts, ᶜspecific.e_tot) if !isnothing(p.sfc_setup) SurfaceConditions.update_surface_conditions!(Y, p, t) diff --git a/src/cache/prognostic_edmf_precomputed_quantities.jl b/src/cache/prognostic_edmf_precomputed_quantities.jl index 1527b1702b..a63578d88f 100644 --- a/src/cache/prognostic_edmf_precomputed_quantities.jl +++ b/src/cache/prognostic_edmf_precomputed_quantities.jl @@ -50,9 +50,8 @@ Updates the draft thermo state and boundary conditions precomputed quantities stored in `p` for edmfx. """ function set_prognostic_edmf_precomputed_quantities_draft_and_bc!(Y, p, ᶠuₕ³, t) - (; energy_form, moisture_model, turbconv_model) = p.atmos + (; 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)) diff --git a/src/callbacks/callbacks.jl b/src/callbacks/callbacks.jl index f0af2081d0..ac8eb18539 100644 --- a/src/callbacks/callbacks.jl +++ b/src/callbacks/callbacks.jl @@ -326,31 +326,26 @@ NVTX.@annotate function compute_diagnostics(integrator) turbulence_convection_diagnostic = NamedTuple() end - if p.atmos.energy_form isa TotalEnergy - sfc_local_geometry = - Fields.level(Fields.local_geometry_field(Y.f), Fields.half) - surface_ct3_unit = - CT3.(unit_basis_vector_data.(CT3, sfc_local_geometry)) - (; ρ_flux_uₕ, ρ_flux_h_tot) = p.precomputed.sfc_conditions - sfc_flux_momentum = - Geometry.UVVector.( - adjoint.(ρ_flux_uₕ ./ Spaces.level(ᶠinterp.(Y.c.ρ), half)) .* - surface_ct3_unit - ) + sfc_local_geometry = + Fields.level(Fields.local_geometry_field(Y.f), Fields.half) + surface_ct3_unit = CT3.(unit_basis_vector_data.(CT3, sfc_local_geometry)) + (; ρ_flux_uₕ, ρ_flux_h_tot) = p.precomputed.sfc_conditions + sfc_flux_momentum = + Geometry.UVVector.( + adjoint.(ρ_flux_uₕ ./ Spaces.level(ᶠinterp.(Y.c.ρ), half)) .* + surface_ct3_unit + ) + vert_diff_diagnostic = (; + sfc_flux_u = sfc_flux_momentum.components.data.:1, + sfc_flux_v = sfc_flux_momentum.components.data.:2, + sfc_flux_energy = dot.(ρ_flux_h_tot, surface_ct3_unit), + ) + if :ρq_tot in propertynames(Y.c) + (; ρ_flux_q_tot) = p.precomputed.sfc_conditions vert_diff_diagnostic = (; - sfc_flux_u = sfc_flux_momentum.components.data.:1, - sfc_flux_v = sfc_flux_momentum.components.data.:2, - sfc_flux_energy = dot.(ρ_flux_h_tot, surface_ct3_unit), + vert_diff_diagnostic..., + sfc_evaporation = dot.(ρ_flux_q_tot, surface_ct3_unit), ) - if :ρq_tot in propertynames(Y.c) - (; ρ_flux_q_tot) = p.precomputed.sfc_conditions - vert_diff_diagnostic = (; - vert_diff_diagnostic..., - sfc_evaporation = dot.(ρ_flux_q_tot, surface_ct3_unit), - ) - end - else - vert_diff_diagnostic = NamedTuple() end if p.atmos.radiation_mode isa RRTMGPI.AbstractRRTMGPMode diff --git a/src/diagnostics/Diagnostics.jl b/src/diagnostics/Diagnostics.jl index 9aee4726de..b9ea4868ac 100644 --- a/src/diagnostics/Diagnostics.jl +++ b/src/diagnostics/Diagnostics.jl @@ -19,9 +19,6 @@ import ..DryModel import ..EquilMoistModel import ..NonEquilMoistModel -# energy_form -import ..TotalEnergy - # precip_model import ..Microphysics0Moment diff --git a/src/diagnostics/core_diagnostics.jl b/src/diagnostics/core_diagnostics.jl index 3595e8a525..24a21b79cf 100644 --- a/src/diagnostics/core_diagnostics.jl +++ b/src/diagnostics/core_diagnostics.jl @@ -418,10 +418,7 @@ add_diagnostic_variable!( ### # Eastward and northward surface drag component (2d) ### -compute_tau!(_, _, _, _, energy_form::T) where {T} = - error_diagnostic_variable("tau", energy_form) - -function compute_tau!(out, state, cache, component, energy_form::TotalEnergy) +function compute_tau!(out, state, cache, component) sfc_local_geometry = Fields.level(Fields.local_geometry_field(state.f), Fields.half) surface_ct3_unit = CT3.(unit_basis_vector_data.(CT3, sfc_local_geometry)) @@ -446,10 +443,8 @@ function compute_tau!(out, state, cache, component, energy_form::TotalEnergy) return end -compute_tauu!(out, state, cache, time) = - compute_tau!(out, state, cache, :1, cache.atmos.energy_form) -compute_tauv!(out, state, cache, time) = - compute_tau!(out, state, cache, :2, cache.atmos.energy_form) +compute_tauu!(out, state, cache, time) = compute_tau!(out, state, cache, :1) +compute_tauv!(out, state, cache, time) = compute_tau!(out, state, cache, :2) add_diagnostic_variable!( short_name = "tauu", @@ -472,12 +467,7 @@ add_diagnostic_variable!( ### # Surface energy flux (2d) - TODO: this may need to be split into sensible and latent heat fluxes ### -compute_hfes!(out, state, cache, time) = - compute_hfes!(out, state, cache, time, cache.atmos.energy_form) -compute_hfes!(_, _, _, _, energy_form::T) where {T} = - error_diagnostic_variable("hfes", energy_form) - -function compute_hfes!(out, state, cache, time, energy_form::TotalEnergy) +function compute_hfes!(out, state, cache, time) (; ρ_flux_h_tot) = cache.precomputed.sfc_conditions sfc_local_geometry = Fields.level(Fields.local_geometry_field(state.f), Fields.half) @@ -500,24 +490,10 @@ add_diagnostic_variable!( ### # Surface evaporation (2d) ### -compute_evspsbl!(out, state, cache, time) = compute_evspsbl!( - out, - state, - cache, - time, - cache.atmos.moisture_model, - cache.atmos.energy_form, -) -compute_evspsbl!( - _, - _, - _, - _, - moisture_model::T1, - energy_form::T2, -) where {T1, T2} = error_diagnostic_variable( - "Can only compute surface_evaporation with energy_form = TotalEnergy() and with a moist model", -) +compute_evspsbl!(out, state, cache, time) = + compute_evspsbl!(out, state, cache, time, cache.atmos.moisture_model) +compute_evspsbl!(_, _, _, _, model::T) where {T} = + error_diagnostic_variable("evspsbl", model) function compute_evspsbl!( out, @@ -525,7 +501,6 @@ function compute_evspsbl!( cache, time, moisture_model::T, - energy_form::TotalEnergy, ) where {T <: Union{EquilMoistModel, NonEquilMoistModel}} (; ρ_flux_q_tot) = cache.precomputed.sfc_conditions sfc_local_geometry = diff --git a/src/diagnostics/default_diagnostics.jl b/src/diagnostics/default_diagnostics.jl index 7bfcea5d70..3a51323e77 100644 --- a/src/diagnostics/default_diagnostics.jl +++ b/src/diagnostics/default_diagnostics.jl @@ -193,7 +193,7 @@ hourly_average(short_names; output_writer = HDF5Writer()) = ######## function core_default_diagnostics() core_diagnostics = - ["ts", "ta", "thetaa", "ha", "pfull", "rhoa", "ua", "va", "wa"] + ["ts", "ta", "thetaa", "ha", "pfull", "rhoa", "ua", "va", "wa", "hfes"] return [ daily_averages(core_diagnostics...)..., @@ -202,16 +202,6 @@ function core_default_diagnostics() ] end -############### -# Energy form # -############### -function default_diagnostics(::TotalEnergy) - total_energy_diagnostics = ["hfes"] - - return [daily_averages(total_energy_diagnostics...)...] -end - - ################## # Moisture model # ################## diff --git a/src/initial_conditions/InitialConditions.jl b/src/initial_conditions/InitialConditions.jl index 78e4c73525..2ee06bf533 100644 --- a/src/initial_conditions/InitialConditions.jl +++ b/src/initial_conditions/InitialConditions.jl @@ -1,8 +1,6 @@ module InitialConditions import ..AtmosModel -import ..PotentialTemperature -import ..TotalEnergy import ..DryModel import ..EquilMoistModel import ..NonEquilMoistModel diff --git a/src/initial_conditions/atmos_state.jl b/src/initial_conditions/atmos_state.jl index 3232d8cf2c..901c2acb13 100644 --- a/src/initial_conditions/atmos_state.jl +++ b/src/initial_conditions/atmos_state.jl @@ -57,15 +57,12 @@ atmos_face_variables(ls, atmos_model) = (; grid_scale_center_variables(ls, atmos_model) = (; ρ = ls.ρ, uₕ = C12(ls.velocity, ls.geometry), - energy_variables(ls, atmos_model.energy_form)..., + energy_variables(ls)..., moisture_variables(ls, atmos_model.moisture_model)..., precip_variables(ls, atmos_model.precip_model, atmos_model.perf_mode)..., ) -# TODO: Rename ρθ to ρθ_liq_ice. -energy_variables(ls, ::PotentialTemperature) = - (; ρθ = ls.ρ * TD.liquid_ice_pottemp(ls.thermo_params, ls.thermo_state)) -energy_variables(ls, ::TotalEnergy) = (; +energy_variables(ls) = (; ρe_tot = ls.ρ * ( TD.internal_energy(ls.thermo_params, ls.thermo_state) + norm_sqr(ls.velocity) / 2 + diff --git a/src/parameterized_tendencies/radiation/held_suarez.jl b/src/parameterized_tendencies/radiation/held_suarez.jl index 03f8309c23..0a1e6db62a 100644 --- a/src/parameterized_tendencies/radiation/held_suarez.jl +++ b/src/parameterized_tendencies/radiation/held_suarez.jl @@ -100,11 +100,6 @@ function forcing_tendency!(Yₜ, Y, p, t, colidx, ::HeldSuarezForcing) ) @. Yₜ.c.uₕ[colidx] -= (k_f * ᶜheight_factor[colidx]) * Y.c.uₕ[colidx] - if :ρθ in propertynames(Y.c) - @. Yₜ.c.ρθ[colidx] -= - ᶜΔρT[colidx] * fast_pow((p_ref_theta / ᶜp[colidx]), κ_d) - elseif :ρe_tot in propertynames(Y.c) - @. Yₜ.c.ρe_tot[colidx] -= ᶜΔρT[colidx] * cv_d - end + @. Yₜ.c.ρe_tot[colidx] -= ᶜΔρT[colidx] * cv_d return nothing end diff --git a/src/parameterized_tendencies/radiation/radiation.jl b/src/parameterized_tendencies/radiation/radiation.jl index 1c23fb9119..719921b32e 100644 --- a/src/parameterized_tendencies/radiation/radiation.jl +++ b/src/parameterized_tendencies/radiation/radiation.jl @@ -229,11 +229,7 @@ end function radiation_tendency!(Yₜ, Y, p, t, colidx, ::RRTMGPI.AbstractRRTMGPMode) (; ᶠradiation_flux) = p.radiation - if :ρθ in propertynames(Y.c) - error("radiation_tendency! not implemented for ρθ") - elseif :ρe_tot in propertynames(Y.c) - @. Yₜ.c.ρe_tot[colidx] -= ᶜdivᵥ(ᶠradiation_flux[colidx]) - end + @. Yₜ.c.ρe_tot[colidx] -= ᶜdivᵥ(ᶠradiation_flux[colidx]) return nothing end @@ -263,7 +259,6 @@ function radiation_tendency!( radiation_mode::RadiationDYCOMS_RF01, ) @assert !(p.atmos.moisture_model isa DryModel) - @assert p.atmos.energy_form isa TotalEnergy (; params) = p (; ᶜspecific, ᶜts) = p.precomputed diff --git a/src/parameterized_tendencies/sponge/viscous_sponge.jl b/src/parameterized_tendencies/sponge/viscous_sponge.jl index 97279acfce..64e9bf53ca 100644 --- a/src/parameterized_tendencies/sponge/viscous_sponge.jl +++ b/src/parameterized_tendencies/sponge/viscous_sponge.jl @@ -25,28 +25,13 @@ function viscous_sponge_cache(Y, viscous_sponge::ViscousSponge) return (; ᶜβ_viscous, ᶠβ_viscous) end -add_viscous_sponge_energy_tendency!(Yₜ, Y, p, t) = - add_viscous_sponge_energy_tendency!(Yₜ, Y, p, t, p.atmos.energy_form) - -function add_viscous_sponge_energy_tendency!(Yₜ, Y, p, t, ::TotalEnergy) +function add_viscous_sponge_energy_tendency!(Yₜ, Y, p, t) (; ᶜβ_viscous) = p.viscous_sponge (; ᶜh_tot) = p.precomputed ᶜρ = Y.c.ρ @. Yₜ.c.ρe_tot += ᶜβ_viscous * wdivₕ(ᶜρ * gradₕ(ᶜh_tot)) end -function add_viscous_sponge_energy_tendency!( - Yₜ, - Y, - p, - t, - ::PotentialTemperature, -) - (; ᶜβ_viscous) = p.viscous_sponge - ᶜρ = Y.c.ρ - @. Yₜ.c.ρθ += ᶜβ_viscous * wdivₕ(ᶜρ * gradₕ(Y.c.ρθ / ᶜρ)) -end - function viscous_sponge_tendency!(Yₜ, Y, p, t, ::ViscousSponge) (; ᶜβ_viscous, ᶠβ_viscous) = p.viscous_sponge ᶜuₕ = Y.c.uₕ diff --git a/src/prognostic_equations/advection.jl b/src/prognostic_equations/advection.jl index 6d400bd4e5..bfe3ee1395 100644 --- a/src/prognostic_equations/advection.jl +++ b/src/prognostic_equations/advection.jl @@ -24,12 +24,8 @@ NVTX.@annotate function horizontal_advection_tendency!(Yₜ, Y, p, t) end end - if :ρθ in propertynames(Y.c) - @. Yₜ.c.ρθ -= wdivₕ(Y.c.ρθ * ᶜu) - elseif :ρe_tot in propertynames(Y.c) - (; ᶜh_tot) = p.precomputed - @. Yₜ.c.ρe_tot -= wdivₕ(Y.c.ρ * ᶜh_tot * ᶜu) - end + (; ᶜh_tot) = p.precomputed + @. Yₜ.c.ρe_tot -= wdivₕ(Y.c.ρ * ᶜh_tot * ᶜu) if p.atmos.turbconv_model isa PrognosticEDMFX for j in 1:n diff --git a/src/prognostic_equations/hyperdiffusion.jl b/src/prognostic_equations/hyperdiffusion.jl index 1e06f37eae..5191da19f6 100644 --- a/src/prognostic_equations/hyperdiffusion.jl +++ b/src/prognostic_equations/hyperdiffusion.jl @@ -84,11 +84,8 @@ NVTX.@annotate function hyperdiffusion_tendency!(Yₜ, Y, p, t) C123(wgradₕ(divₕ(p.precomputed.ᶜu))) - C123(wcurlₕ(C123(curlₕ(p.precomputed.ᶜu)))) - if :θ in propertynames(ᶜspecific) - @. ᶜ∇²specific_energy = wdivₕ(gradₕ(ᶜspecific.θ)) - elseif :e_tot in propertynames(ᶜspecific) - @. ᶜ∇²specific_energy = wdivₕ(gradₕ(ᶜspecific.e_tot + ᶜp / Y.c.ρ)) - end + @. ᶜ∇²specific_energy = wdivₕ(gradₕ(ᶜspecific.e_tot + ᶜp / Y.c.ρ)) + if diffuse_tke @. ᶜ∇²tke⁰ = wdivₕ(gradₕ(ᶜtke⁰)) end @@ -144,8 +141,7 @@ NVTX.@annotate function hyperdiffusion_tendency!(Yₜ, Y, p, t) @. Yₜ.c.uₕ -= κ₄ * C12(ᶜ∇²u) @. Yₜ.f.u₃ -= κ₄ * ᶠwinterp(ᶜJ * Y.c.ρ, C3(ᶜ∇²u)) - ᶜρ_energyₜ = :θ in propertynames(ᶜspecific) ? Yₜ.c.ρθ : Yₜ.c.ρe_tot - @. ᶜρ_energyₜ -= κ₄ * wdivₕ(Y.c.ρ * gradₕ(ᶜ∇²specific_energy)) + @. Yₜ.c.ρe_tot -= κ₄ * wdivₕ(Y.c.ρ * gradₕ(ᶜ∇²specific_energy)) # Sub-grid scale hyperdiffusion continued if (turbconv_model isa PrognosticEDMFX) && diffuse_tke diff --git a/src/prognostic_equations/implicit/implicit_solver.jl b/src/prognostic_equations/implicit/implicit_solver.jl index e7c89ae21c..3ef3ac671f 100644 --- a/src/prognostic_equations/implicit/implicit_solver.jl +++ b/src/prognostic_equations/implicit/implicit_solver.jl @@ -112,8 +112,7 @@ function ImplicitEquationJacobian( is_in_Y(name) = MatrixFields.has_field(Y, name) - energy_name = - MatrixFields.unrolled_findonly(is_in_Y, (@name(c.ρe_tot), @name(c.ρθ))) + energy_name = @name(c.ρe_tot) tracer_names = ( @name(c.ρq_tot), @name(c.ρq_liq), @@ -206,7 +205,7 @@ _linsolve!(x, A, b, update_matrix = false; kwargs...) = ldiv!(x, A, b) function Wfact!(A, Y, p, dtγ, t) NVTX.@range "Wfact!" color = colorant"green" begin # Remove unnecessary values from p to avoid allocations in bycolumn. - (; energy_form, rayleigh_sponge) = p.atmos + (; rayleigh_sponge) = p.atmos p′ = (; p.precomputed.ᶜspecific, p.precomputed.ᶠu³, @@ -220,7 +219,7 @@ function Wfact!(A, Y, p, dtγ, t) p.scratch.ᶜtemp_scalar, p.params, p.atmos, - (energy_form isa TotalEnergy ? (; p.precomputed.ᶜh_tot) : (;))..., + p.precomputed.ᶜh_tot, ( rayleigh_sponge isa RayleighSponge ? (; p.rayleigh_sponge.ᶠβ_rayleigh_w) : (;) @@ -273,17 +272,9 @@ function update_implicit_equation_jacobian!(A, Y, p, dtγ, colidx) ᶜinterp_matrix() ⋅ DiagonalMatrixRow(adjoint(CT3(ᶠu₃[colidx]))) + DiagonalMatrixRow(adjoint(CT3(ᶜuₕ[colidx]))) ⋅ ᶜinterp_matrix() - # We can express the pressure as either - # ᶜp = p_ref_theta * (ᶜρθ * R_d / p_ref_theta)^(1 / (1 - κ_d)) or + # We can express the pressure as # ᶜp = R_d * (ᶜρe_tot / cv_d + ᶜρ * (T_tri - (ᶜK + ᶜΦ) / cv_d)) + O(ᶜq_tot). - # In the first case, we find that - # ∂(ᶜp)/∂(ᶜρ) = 0I, - # ∂(ᶜp)/∂(ᶜK) = 0I, and - # ∂(ᶜp)/∂(ᶜρθ) = - # DiagonalMatrixRow( - # R_d / (1 - κ_d) * (ᶜρθ * R_d / p_ref_theta)^(κ_d / (1 - κ_d)) - # ). - # In the second case, we can ignore all O(ᶜq_tot) terms to approximate + # we can ignore all O(ᶜq_tot) terms to approximate # ∂(ᶜp)/∂(ᶜρ) ≈ DiagonalMatrixRow(R_d * (T_tri - (ᶜK + ᶜΦ) / cv_d)), # ∂(ᶜp)/∂(ᶜK) ≈ DiagonalMatrixRow(-R_d * ᶜρ / cv_d), and # ∂(ᶜp)/∂(ᶜρe_tot) ≈ DiagonalMatrixRow(R_d / cv_d). @@ -323,7 +314,6 @@ function update_implicit_equation_jacobian!(A, Y, p, dtγ, colidx) u³_sign(u³) = sign(u³.u³) # sign of the scalar value stored inside u³ tracer_info = ( (@name(c.ρ), @name(ᶜ1), density_upwinding), - (@name(c.ρθ), @name(ᶜspecific.θ), energy_upwinding), (@name(c.ρe_tot), @name(ᶜh_tot), energy_upwinding), (@name(c.ρq_tot), @name(ᶜspecific.q_tot), tracer_upwinding), (@name(c.ρq_liq), @name(ᶜspecific.q_liq), tracer_upwinding), @@ -435,13 +425,7 @@ function update_implicit_equation_jacobian!(A, Y, p, dtγ, colidx) # (ᶠgradᵥ(ᶜp - ᶜp_ref) - ᶠinterp(ᶜρ_ref) * ᶠgradᵥ_ᶜΦ) / # ᶠinterp(ᶜρ)^2 # ) ⋅ ᶠinterp_matrix(). - # If the pressure is computed using potential temperature, then - # ∂(ᶠu₃ₜ)/∂(ᶜρθ) = - # ∂(ᶠu₃ₜ)/∂(ᶠgradᵥ(ᶜp - ᶜp_ref)) ⋅ ∂(ᶠgradᵥ(ᶜp - ᶜp_ref))/∂(ᶜρθ) = - # DiagonalMatrixRow(-1 / ᶠinterp(ᶜρ)) ⋅ ᶠgradᵥ_matrix() ⋅ - # ∂(ᶜp)/∂(ᶜρθ) and - # ∂(ᶠu₃ₜ)/∂(ᶠu₃) = 0I. - # If the pressure is computed using total energy, then + # The pressure is computed using total energy, therefore # ∂(ᶠu₃ₜ)/∂(ᶜρe_tot) = # ∂(ᶠu₃ₜ)/∂(ᶠgradᵥ(ᶜp - ᶜp_ref)) ⋅ ∂(ᶠgradᵥ(ᶜp - ᶜp_ref))/∂(ᶜρe_tot) = # DiagonalMatrixRow(-1 / ᶠinterp(ᶜρ)) ⋅ ᶠgradᵥ_matrix() ⋅ @@ -462,63 +446,34 @@ function update_implicit_equation_jacobian!(A, Y, p, dtγ, colidx) ∂ᶠu₃_err_∂ᶜρ = matrix[@name(f.u₃), @name(c.ρ)] ∂ᶠu₃_err_∂ᶠu₃ = matrix[@name(f.u₃), @name(f.u₃)] I_u₃ = DiagonalMatrixRow(one_C3xACT3) - if MatrixFields.has_field(Y, @name(c.ρθ)) - ᶜρθ = Y.c.ρθ - ∂ᶠu₃_err_∂ᶜρθ = matrix[@name(f.u₃), @name(c.ρθ)] - @. ∂ᶠu₃_err_∂ᶜρ[colidx] = - dtγ * DiagonalMatrixRow( + ∂ᶠu₃_err_∂ᶜρe_tot = matrix[@name(f.u₃), @name(c.ρe_tot)] + @. ∂ᶠu₃_err_∂ᶜρ[colidx] = + dtγ * ( + DiagonalMatrixRow(-1 / ᶠinterp(ᶜρ[colidx])) ⋅ ᶠgradᵥ_matrix() ⋅ + DiagonalMatrixRow( + R_d * (T_tri - (ᶜK[colidx] + ᶜΦ[colidx]) / cv_d), + ) + + DiagonalMatrixRow( ( ᶠgradᵥ(ᶜp[colidx] - ᶜp_ref[colidx]) - ᶠinterp(ᶜρ_ref[colidx]) * ᶠgradᵥ_ᶜΦ[colidx] ) / abs2(ᶠinterp(ᶜρ[colidx])), ) ⋅ ᶠinterp_matrix() - @. ∂ᶠu₃_err_∂ᶜρθ[colidx] = - dtγ * DiagonalMatrixRow(-1 / ᶠinterp(ᶜρ[colidx])) ⋅ - ᶠgradᵥ_matrix() ⋅ DiagonalMatrixRow( - R_d / (1 - κ_d) * - (ᶜρθ[colidx] * R_d / p_ref_theta)^(κ_d / (1 - κ_d)), - ) - if p.atmos.rayleigh_sponge isa RayleighSponge - @. ∂ᶠu₃_err_∂ᶠu₃[colidx] = - dtγ * - DiagonalMatrixRow(-p.ᶠβ_rayleigh_w[colidx] * (one_C3xACT3,)) - - (I_u₃,) - else - ∂ᶠu₃_err_∂ᶠu₃[colidx] .= (-I_u₃,) - end - elseif MatrixFields.has_field(Y, @name(c.ρe_tot)) - ∂ᶠu₃_err_∂ᶜρe_tot = matrix[@name(f.u₃), @name(c.ρe_tot)] - @. ∂ᶠu₃_err_∂ᶜρ[colidx] = + ) + @. ∂ᶠu₃_err_∂ᶜρe_tot[colidx] = + dtγ * DiagonalMatrixRow(-1 / ᶠinterp(ᶜρ[colidx])) ⋅ ᶠgradᵥ_matrix() * + R_d / cv_d + if p.atmos.rayleigh_sponge isa RayleighSponge + @. ∂ᶠu₃_err_∂ᶠu₃[colidx] = dtγ * ( DiagonalMatrixRow(-1 / ᶠinterp(ᶜρ[colidx])) ⋅ ᶠgradᵥ_matrix() ⋅ - DiagonalMatrixRow( - R_d * (T_tri - (ᶜK[colidx] + ᶜΦ[colidx]) / cv_d), - ) + - DiagonalMatrixRow( - ( - ᶠgradᵥ(ᶜp[colidx] - ᶜp_ref[colidx]) - - ᶠinterp(ᶜρ_ref[colidx]) * ᶠgradᵥ_ᶜΦ[colidx] - ) / abs2(ᶠinterp(ᶜρ[colidx])), - ) ⋅ ᶠinterp_matrix() - ) - @. ∂ᶠu₃_err_∂ᶜρe_tot[colidx] = + DiagonalMatrixRow(-R_d * ᶜρ[colidx] / cv_d) ⋅ ∂ᶜK_∂ᶠu₃[colidx] + + DiagonalMatrixRow(-p.ᶠβ_rayleigh_w[colidx] * (one_C3xACT3,)) + ) - (I_u₃,) + else + @. ∂ᶠu₃_err_∂ᶠu₃[colidx] = dtγ * DiagonalMatrixRow(-1 / ᶠinterp(ᶜρ[colidx])) ⋅ - ᶠgradᵥ_matrix() * R_d / cv_d - if p.atmos.rayleigh_sponge isa RayleighSponge - @. ∂ᶠu₃_err_∂ᶠu₃[colidx] = - dtγ * ( - DiagonalMatrixRow(-1 / ᶠinterp(ᶜρ[colidx])) ⋅ - ᶠgradᵥ_matrix() ⋅ - DiagonalMatrixRow(-R_d * ᶜρ[colidx] / cv_d) ⋅ - ∂ᶜK_∂ᶠu₃[colidx] + DiagonalMatrixRow( - -p.ᶠβ_rayleigh_w[colidx] * (one_C3xACT3,), - ) - ) - (I_u₃,) - else - @. ∂ᶠu₃_err_∂ᶠu₃[colidx] = - dtγ * DiagonalMatrixRow(-1 / ᶠinterp(ᶜρ[colidx])) ⋅ - ᶠgradᵥ_matrix() ⋅ DiagonalMatrixRow(-R_d * ᶜρ[colidx] / cv_d) ⋅ - ∂ᶜK_∂ᶠu₃[colidx] - (I_u₃,) - end + ᶠgradᵥ_matrix() ⋅ DiagonalMatrixRow(-R_d * ᶜρ[colidx] / cv_d) ⋅ + ∂ᶜK_∂ᶠu₃[colidx] - (I_u₃,) end end diff --git a/src/prognostic_equations/implicit/implicit_tendency.jl b/src/prognostic_equations/implicit/implicit_tendency.jl index 323332a983..1eb6b34bcf 100644 --- a/src/prognostic_equations/implicit/implicit_tendency.jl +++ b/src/prognostic_equations/implicit/implicit_tendency.jl @@ -98,7 +98,7 @@ function implicit_vertical_advection_tendency!(Yₜ, Y, p, t, colidx) ᶠu³[colidx], ᶜχ[colidx], dt, - χ_name == :θ ? energy_upwinding : tracer_upwinding, + tracer_upwinding, ) end diff --git a/src/prognostic_equations/vertical_diffusion_boundary_layer.jl b/src/prognostic_equations/vertical_diffusion_boundary_layer.jl index 92cbee09ed..a949eb9617 100644 --- a/src/prognostic_equations/vertical_diffusion_boundary_layer.jl +++ b/src/prognostic_equations/vertical_diffusion_boundary_layer.jl @@ -120,8 +120,6 @@ function vertical_diffusion_boundary_layer_tendency!( χ_name == :e_tot && continue if χ_name == :q_tot @. ρ_flux_χ[colidx] = sfc_conditions.ρ_flux_q_tot[colidx] - elseif χ_name == :θ - @. ρ_flux_χ[colidx] = sfc_conditions.ρ_flux_θ[colidx] else @. ρ_flux_χ[colidx] = C3(FT(0)) end diff --git a/src/solver/model_getters.jl b/src/solver/model_getters.jl index 70a01e9f1d..ea76d26efe 100644 --- a/src/solver/model_getters.jl +++ b/src/solver/model_getters.jl @@ -157,19 +157,6 @@ function get_perf_mode(parsed_args) end end -function get_energy_form(parsed_args, vert_diff) - energy_name = parsed_args["energy_name"] - @assert energy_name in ("rhoe", "rhotheta") - if !isnothing(vert_diff) - @assert energy_name == "rhoe" - end - return if energy_name == "rhoe" - TotalEnergy() - elseif energy_name == "rhotheta" - PotentialTemperature() - end -end - function get_radiation_mode(parsed_args, ::Type{FT}) where {FT} idealized_h2o = parsed_args["idealized_h2o"] @assert idealized_h2o in (true, false) diff --git a/src/solver/type_getters.jl b/src/solver/type_getters.jl index 881d8fce75..362f1928a0 100644 --- a/src/solver/type_getters.jl +++ b/src/solver/type_getters.jl @@ -44,7 +44,6 @@ function get_atmos(config::AtmosConfig, params) moisture_model, model_config, perf_mode = get_perf_mode(parsed_args), - energy_form = get_energy_form(parsed_args, vert_diff), radiation_mode, subsidence = get_subsidence_model(parsed_args, radiation_mode, FT), ls_adv = get_large_scale_advection_model(parsed_args, FT), diff --git a/src/solver/types.jl b/src/solver/types.jl index 8dd5747da7..ab882bcba8 100644 --- a/src/solver/types.jl +++ b/src/solver/types.jl @@ -6,10 +6,6 @@ struct DryModel <: AbstractMoistureModel end struct EquilMoistModel <: AbstractMoistureModel end struct NonEquilMoistModel <: AbstractMoistureModel end -abstract type AbstractEnergyFormulation end -struct PotentialTemperature <: AbstractEnergyFormulation end -struct TotalEnergy <: AbstractEnergyFormulation end - abstract type AbstractPrecipitationModel end struct NoPrecipitation <: AbstractPrecipitationModel end struct Microphysics0Moment <: AbstractPrecipitationModel end @@ -230,7 +226,6 @@ struct GCMSurfaceThermoState <: AbstractSurfaceThermoState end # Define broadcasting for types Base.broadcastable(x::AbstractSurfaceThermoState) = tuple(x) Base.broadcastable(x::AbstractMoistureModel) = tuple(x) -Base.broadcastable(x::AbstractEnergyFormulation) = tuple(x) Base.broadcastable(x::AbstractPrecipitationModel) = tuple(x) Base.broadcastable(x::AbstractForcing) = tuple(x) Base.broadcastable(x::PrognosticEDMFX) = tuple(x) @@ -323,7 +318,6 @@ Base.@kwdef struct AtmosModel{ MC, PEM, MM, - EF, PM, F, S, @@ -350,7 +344,6 @@ Base.@kwdef struct AtmosModel{ model_config::MC = nothing perf_mode::PEM = nothing moisture_model::MM = nothing - energy_form::EF = nothing precip_model::PM = nothing forcing_type::F = nothing subsidence::S = nothing diff --git a/src/surface_conditions/SurfaceConditions.jl b/src/surface_conditions/SurfaceConditions.jl index 8a73cabd07..f998197bc8 100644 --- a/src/surface_conditions/SurfaceConditions.jl +++ b/src/surface_conditions/SurfaceConditions.jl @@ -2,8 +2,6 @@ module SurfaceConditions import ..InitialConditions as ICs import ..Parameters as CAP -import ..PotentialTemperature -import ..TotalEnergy import ..DryModel import ..ZonallyAsymmetricSST import ..ZonallySymmetricSST diff --git a/src/surface_conditions/surface_conditions.jl b/src/surface_conditions/surface_conditions.jl index 7ff2723ce8..0c83b28354 100644 --- a/src/surface_conditions/surface_conditions.jl +++ b/src/surface_conditions/surface_conditions.jl @@ -91,9 +91,7 @@ function set_dummy_surface_conditions!(p) ) @. sfc_conditions.ρ_flux_q_tot = C3(FT(0)) end - if atmos.energy_form isa TotalEnergy - @. sfc_conditions.ρ_flux_h_tot = C3(FT(0)) - end + @. sfc_conditions.ρ_flux_h_tot = C3(FT(0)) end """ @@ -384,11 +382,7 @@ function atmos_surface_conditions( thermo_params = CAP.thermodynamics_params(params) surface_normal = C3(unit_basis_vector_data(C3, surface_local_geometry)) - energy_flux = if atmos.energy_form isa PotentialTemperature - (; ρ_flux_θ = shf / TD.cp_m(thermo_params, ts) * surface_normal) - elseif atmos.energy_form isa TotalEnergy - (; ρ_flux_h_tot = (shf + lhf) * surface_normal) - end + energy_flux = (; ρ_flux_h_tot = (shf + lhf) * surface_normal) moisture_flux = atmos.moisture_model isa DryModel ? (;) : (; ρ_flux_q_tot = evaporation * surface_normal) @@ -415,16 +409,12 @@ function atmos_surface_conditions( end """ - surface_conditions_type(moisture_model, energy_form, FT) + surface_conditions_type(moisture_model, FT) Gets the return type of `surface_conditions` without evaluating the function. """ function surface_conditions_type(atmos, ::Type{FT}) where {FT} - energy_flux_names = if atmos.energy_form isa PotentialTemperature - (:ρ_flux_θ,) - elseif atmos.energy_form isa TotalEnergy - (:ρ_flux_h_tot,) - end + energy_flux_names = (:ρ_flux_h_tot,) moisture_flux_names = atmos.moisture_model isa DryModel ? () : (:ρ_flux_q_tot,) names = ( diff --git a/src/utils/discrete_hydrostatic_balance.jl b/src/utils/discrete_hydrostatic_balance.jl index 8706471804..1f37af47ac 100644 --- a/src/utils/discrete_hydrostatic_balance.jl +++ b/src/utils/discrete_hydrostatic_balance.jl @@ -49,7 +49,6 @@ function set_discrete_hydrostatic_balanced_state!(Y, p) Geometry.UVWVector(Y.c.uₕ) + Geometry.UVWVector(ᶜinterp(Y.f.u₃)), ), - p.atmos.energy_form, ), ) # broadcasting doesn't seem to work with kwargs, so define interim ls() end diff --git a/src/utils/utilities.jl b/src/utils/utilities.jl index b37de4a8c2..d76e74537a 100644 --- a/src/utils/utilities.jl +++ b/src/utils/utilities.jl @@ -5,7 +5,7 @@ import ClimaComms import ClimaCore: Spaces, Topologies, Fields, Geometry import LinearAlgebra: norm_sqr -is_energy_var(symbol) = symbol in (:ρθ, :ρe_tot, :ρaθ, :ρae_tot) +is_energy_var(symbol) = symbol in (:ρe_tot, :ρae_tot) is_momentum_var(symbol) = symbol in (:uₕ, :ρuₕ, :u₃, :ρw) is_turbconv_var(symbol) = symbol in (:turbconv, :sgsʲs, :sgs⁰) is_tracer_var(symbol) = !( diff --git a/src/utils/variable_manipulations.jl b/src/utils/variable_manipulations.jl index a543b63da8..c6711e285c 100644 --- a/src/utils/variable_manipulations.jl +++ b/src/utils/variable_manipulations.jl @@ -249,11 +249,11 @@ u₃⁰(ρaʲs, u₃ʲs, ρ, u₃, turbconv_model) = divide_by_ρa( """ remove_energy_var(specific_state) -Creates a copy of `specific_state` with the energy variable (`θ` or `e_tot`) +Creates a copy of `specific_state` with the energy variable removed, where `specific_state` is the result of calling, e.g., `specific_gs`, `specific_sgsʲs`, or `specific_sgs⁰`. """ remove_energy_var(specific_state::NamedTuple) = - Base.structdiff(specific_state, NamedTuple{(:θ, :e_tot)}) + Base.structdiff(specific_state, NamedTuple{(:e_tot,)}) remove_energy_var(specific_state::Tuple) = map(remove_energy_var, specific_state)