From a82a1075047e9ef4e3bece2a568e122127fca76d Mon Sep 17 00:00:00 2001 From: Haakon Ludvig Langeland Ervik <45243236+haakon-e@users.noreply.github.com> Date: Fri, 24 Oct 2025 10:16:01 -0700 Subject: [PATCH] rft: smagorinsky-lilly for future h-v splitting --- config/default_configs/default_config.yml | 4 +- .../box_density_current_test.yml | 2 +- config/model_configs/les_isdac_box.yml | 2 +- post_processing/ci_plots.jl | 21 +- src/cache/precomputed_quantities.jl | 14 +- src/diagnostics/core_diagnostics.jl | 40 +++ .../les_sgs_models/smagorinsky_lilly.jl | 235 +++++++++--------- .../implicit/manual_sparse_jacobian.jl | 6 +- src/solver/model_getters.jl | 16 +- src/solver/types.jl | 43 +++- 10 files changed, 233 insertions(+), 150 deletions(-) diff --git a/config/default_configs/default_config.yml b/config/default_configs/default_config.yml index 6acbda7c9a..374198faec 100644 --- a/config/default_configs/default_config.yml +++ b/config/default_configs/default_config.yml @@ -170,8 +170,8 @@ c_amd: help: "Model coefficient for AMD-LES closure (TODO: Move to parameters.toml)" value: 0.29 smagorinsky_lilly: - help: "Smagorinsky-Lilly diffusive closure [`false` (default), `true`]" - value: false + help: "Smagorinsky-Lilly diffusive closure [`nothing` (default), `UVW`]" + value: ~ bubble: help: "Enable bubble correction for more accurate surface areas" value: false diff --git a/config/model_configs/box_density_current_test.yml b/config/model_configs/box_density_current_test.yml index 42fdf89e49..24ab27003a 100644 --- a/config/model_configs/box_density_current_test.yml +++ b/config/model_configs/box_density_current_test.yml @@ -1,7 +1,7 @@ reference_job_id: "les_box" initial_condition: "DryDensityCurrentProfile" config: "box" -smagorinsky_lilly: true +smagorinsky_lilly: "UVW" discrete_hydrostatic_balance: true hyperdiff: "false" x_max: 51200.0 diff --git a/config/model_configs/les_isdac_box.yml b/config/model_configs/les_isdac_box.yml index dd67a1bc5e..c318a1f681 100644 --- a/config/model_configs/les_isdac_box.yml +++ b/config/model_configs/les_isdac_box.yml @@ -14,7 +14,7 @@ implicit_diffusion: false approximate_linear_solve_iters: 2 hyperdiff: "false" apply_limiter: false -smagorinsky_lilly: true +smagorinsky_lilly: "UVW" # time- and spatial discretization x_elem: 10 x_max: 3.2e3 diff --git a/post_processing/ci_plots.jl b/post_processing/ci_plots.jl index 43cb4ecc21..828b02ba9f 100644 --- a/post_processing/ci_plots.jl +++ b/post_processing/ci_plots.jl @@ -1310,17 +1310,11 @@ function make_plots( reduction = "inst" short_names = [ - "wa", - "ua", - "va", - "ta", - "thetaa", - "ha", - "hus", - "hur", - "cl", - "clw", - "cli", + "wa", "ua", "va", "ta", "thetaa", "ha", + "hus", "hur", "cl", "clw", "cli", "ke", + "Dh_smag", "strainh_smag", # smag horizontal + "Dv_smag", "strainv_smag", # smag vertical + "edt", # DecayWithHeight vertical diffusivity ] short_names = short_names ∩ collect(keys(simdirs[1].vars)) @@ -1333,10 +1327,7 @@ function make_plots( window_end = last(var.dims["time"]) window_start = window_end - 2hours var_window = ClimaAnalysis.window( - var, - "time"; - left = window_start, - right = window_end, + var, "time"; left = window_start, right = window_end, ) var_reduced = horizontal_average(average_time(var_window)) return var_reduced diff --git a/src/cache/precomputed_quantities.jl b/src/cache/precomputed_quantities.jl index 1fca67166e..30e6cc0c54 100644 --- a/src/cache/precomputed_quantities.jl +++ b/src/cache/precomputed_quantities.jl @@ -249,10 +249,12 @@ function precomputed_quantities(Y, atmos) if atmos.smagorinsky_lilly isa SmagorinskyLilly uvw_vec = UVW(FT(0), FT(0), FT(0)) (; - ᶜτ_smag = similar(Y.c, typeof(uvw_vec * uvw_vec')), - ᶠτ_smag = similar(Y.f, typeof(uvw_vec * uvw_vec')), - ᶜD_smag = similar(Y.c, FT), - ᶠD_smag = similar(Y.f, FT), + ᶜS = similar(Y.c, typeof(uvw_vec * uvw_vec')), + ᶠS = similar(Y.f, typeof(uvw_vec * uvw_vec')), + ᶜS_norm_h = similar(Y.c, FT), ᶜS_norm_v = similar(Y.c, FT), + ᶜL_h = similar(Y.c, FT), ᶜL_v = similar(Y.c, FT), + ᶜνₜ_h = similar(Y.c, FT), ᶜνₜ_v = similar(Y.c, FT), + ᶜD_h = similar(Y.c, FT), ᶜD_v = similar(Y.c, FT), ) else (;) @@ -608,9 +610,7 @@ NVTX.@annotate function set_explicit_precomputed_quantities_part2!(Y, p, t) set_cloud_fraction!(Y, p, moisture_model, cloud_model) end - if p.atmos.smagorinsky_lilly isa SmagorinskyLilly - set_smagorinsky_lilly_precomputed_quantities!(Y, p) - end + set_smagorinsky_lilly_precomputed_quantities!(Y, p, p.atmos.smagorinsky_lilly) if p.atmos.amd_les isa AnisotropicMinimumDissipation set_amd_precomputed_quantities!(Y, p) diff --git a/src/diagnostics/core_diagnostics.jl b/src/diagnostics/core_diagnostics.jl index 164e387e24..3942071af0 100644 --- a/src/diagnostics/core_diagnostics.jl +++ b/src/diagnostics/core_diagnostics.jl @@ -343,6 +343,46 @@ add_diagnostic_variable!( end, ) +### +# Smagorinsky Lilly diffusivity +### +add_diagnostic_variable!( + short_name = "Dh_smag", + long_name = "Horizontal smagorinsky diffusivity", + units = "m^2 s^-1", + compute! = (out, _, cache, _) -> begin + (; ᶜD_h) = cache.precomputed + isnothing(out) ? copy(ᶜD_h) : (out .= ᶜD_h) + end, +) +add_diagnostic_variable!( + short_name = "Dv_smag", + long_name = "Vertical smagorinsky diffusivity", + units = "m^2 s^-1", + compute! = (out, _, cache, _) -> begin + (; ᶜD_v) = cache.precomputed + isnothing(out) ? copy(ᶜD_v) : (out .= ᶜD_v) + end, +) +add_diagnostic_variable!( + short_name = "strainh_smag", + long_name = "Horizontal strain rate magnitude (for Smagorinsky)", + units = "s", + compute! = (out, state, cache, _) -> begin + (; ᶜS_norm_h) = cache.precomputed + isnothing(out) ? copy(ᶜS_norm_h) : (out .= ᶜS_norm_h) + end, +) +add_diagnostic_variable!( + short_name = "strainv_smag", + long_name = "Vertical strain rate magnitude (for Smagorinsky)", + units = "s", + compute! = (out, state, cache, _) -> begin + (; ᶜS_norm_v) = cache.precomputed + isnothing(out) ? copy(ᶜS_norm_v) : (out .= ᶜS_norm_v) + end, +) + ### # Relative humidity (3d) ### diff --git a/src/parameterized_tendencies/les_sgs_models/smagorinsky_lilly.jl b/src/parameterized_tendencies/les_sgs_models/smagorinsky_lilly.jl index 962f4fb767..acd2d6e572 100644 --- a/src/parameterized_tendencies/les_sgs_models/smagorinsky_lilly.jl +++ b/src/parameterized_tendencies/les_sgs_models/smagorinsky_lilly.jl @@ -7,179 +7,176 @@ import ClimaCore.Operators as Operators import ClimaCore: Geometry """ - set_smagorinsky_lilly_precomputed_quantities!(Y, p) + lilly_stratification_correction(p, ᶜS) + +Return a lazy representation of the Lilly stratification correction factor + based on the local Richardson number. + +# Arguments +- `p`: The model parameters, e.g. `AtmosCache`. +- `ᶜS`: The cell-centered strain rate tensor. +""" +function lilly_stratification_correction(p, ᶜS) + (; ᶜts) = p.precomputed + (; ᶜtemp_scalar) = p.scratch + grav = CAP.grav(p.params) + Pr_t = CAP.Prandtl_number_0(CAP.turbconv_params(p.params)) + thermo_params = CAP.thermodynamics_params(p.params) + FT = eltype(Pr_t) + # Stratification correction + ᶜθ_v = @. lazy(TD.virtual_pottemp(thermo_params, ᶜts)) + ᶜ∇ᵥθ = @. ᶜtemp_scalar = Geometry.WVector(ᶜgradᵥ(ᶠinterp(ᶜθ_v))).components.data.:1 + ᶜN² = @. lazy(grav / ᶜθ_v * ᶜ∇ᵥθ) + ᶜS_norm = strain_rate_norm(ᶜS, Geometry.WAxis()) + + ᶜRi = @. lazy(ᶜN² / (ᶜS_norm^2 + eps(FT))) # Ri = N² / |S|² + ᶜfb = @. lazy(ifelse(ᶜRi ≤ 0, FT(1), max(0, 1 - ᶜRi / Pr_t)^(1 // 4))) +end -Compute the Smagorinsky-Lilly diffusivity tensors, `ᶜτ_smag`, `ᶠτ_smag`, `ᶜD_smag`, and `ᶠD_smag`. -Store in the precomputed quantities `p.precomputed`. +""" + set_smagorinsky_lilly_precomputed_quantities!(Y, p) -The subgrid-scale momentum flux tensor is defined by `τ = -2 νₜ ∘ S`, where `νₜ` is the Smagorinsky-Lilly eddy viscosity -and `S` is the strain rate tensor. +Compute the Smagorinsky-Lilly horizontal and vertical quantities needed for + subgrid-scale diffusive tendencies -The turbulent diffusivity is defined as `D = νₜ / Pr_t`, where `Pr_t` is the turbulent Prandtl number for neutral -stratification. +The subgrid-scale momentum flux tensor is defined by `τ = -2 νₜ ∘ S`, +where `νₜ` is the Smagorinsky-Lilly eddy viscosity and `S` is the strain rate tensor. -These quantities are computed for both cell centers and faces, with prefixes `ᶜ` and `ᶠ`, respectively. +The turbulent diffusivity is defined as `D = νₜ / Pr_t`, +where `Pr_t` is the turbulent Prandtl number for neutral stratification. +This method precomputes and stores in `p.precomputed` the following quantities: +- strain on centers and faces: `ᶜS`, `ᶠS` +- horizontal and vertical strain rate norm, eddy viscosities, and diffusivities, on centers: + - `ᶜS_norm_h`, `ᶜS_norm_v`, `ᶜνₜ_h`, `ᶜνₜ_v`, `ᶜD_h`, `ᶜD_v` # Arguments - `Y`: The model state. - `p`: The model parameters, e.g. `AtmosCache`. +- `model`: The Smagorinsky model type """ -function set_smagorinsky_lilly_precomputed_quantities!(Y, p) - - (; atmos, precomputed, scratch, params) = p - c_smag = CAP.c_smag(params) - Pr_t = CAP.Prandtl_number_0(CAP.turbconv_params(params)) - (; ᶜu, ᶠu³, ᶜts, ᶜτ_smag, ᶠτ_smag, ᶜD_smag, ᶠD_smag) = precomputed - FT = eltype(Y) - grav = CAP.grav(params) - thermo_params = CAP.thermodynamics_params(params) - (; ᶜtemp_UVWxUVW, ᶠtemp_UVWxUVW, ᶜtemp_strain, ᶠtemp_strain) = scratch - (; ᶜtemp_scalar, ᶜtemp_scalar_2, ᶠtemp_scalar, ᶜtemp_UVW, ᶠtemp_UVW) = - scratch - - ∇ᵥuvw_boundary = Geometry.outer(Geometry.WVector(0), UVW(0, 0, 0)) - ᶠgradᵥ_uvw = Operators.GradientC2F( - bottom = Operators.SetGradient(∇ᵥuvw_boundary), - top = Operators.SetGradient(∇ᵥuvw_boundary), - ) - axis_uvw = (Geometry.UVWAxis(),) - - # Compute UVW velocities - ᶜu_uvw = @. ᶜtemp_UVW = UVW(ᶜu) - ᶠu_uvw = @. ᶠtemp_UVW = UVW(ᶠinterp(Y.c.uₕ)) + UVW(ᶠu³) - - # Gradients - ## cell centers - ∇ᶜu_uvw = @. ᶜtemp_UVWxUVW = Geometry.project(axis_uvw, ᶜgradᵥ(ᶠu_uvw)) # vertical component - @. ∇ᶜu_uvw += Geometry.project(axis_uvw, gradₕ(ᶜu_uvw)) # horizontal component - ## cell faces - ∇ᶠu_uvw = @. ᶠtemp_UVWxUVW = Geometry.project(axis_uvw, ᶠgradᵥ_uvw(ᶜu_uvw)) # vertical component - @. ∇ᶠu_uvw += Geometry.project(axis_uvw, gradₕ(ᶠu_uvw)) # horizontal component - - # Strain rate tensor - ᶜS = @. ᶜtemp_strain = (∇ᶜu_uvw + adjoint(∇ᶜu_uvw)) / 2 - ᶠS = @. ᶠtemp_strain = (∇ᶠu_uvw + adjoint(∇ᶠu_uvw)) / 2 - - # Stratification correction - ᶜθ_v = @. ᶜtemp_scalar = TD.virtual_pottemp(thermo_params, ᶜts) - ᶜ∇ᵥθ = @. ᶜtemp_scalar_2 = - Geometry.WVector(ᶜgradᵥ(ᶠinterp(ᶜθ_v))).components.data.:1 - ᶜN² = @. ᶜtemp_scalar = grav / ᶜθ_v * ᶜ∇ᵥθ - ᶜS_norm = @. ᶜtemp_scalar_2 = √(2 * CA.norm_sqr(ᶜS)) +function set_smagorinsky_lilly_precomputed_quantities!(Y, p, model) + (; ᶜu, ᶠu, ᶜS, ᶠS, ᶜL_h, ᶜL_v, ᶜS_norm_h, ᶜS_norm_v, ᶜνₜ_h, ᶜνₜ_v, ᶜD_h, ᶜD_v) = + p.precomputed + (; ᶜtemp_scalar) = p.scratch + c_smag = CAP.c_smag(p.params) - ᶜRi = @. ᶜtemp_scalar = ᶜN² / (ᶜS_norm^2 + eps(FT)) # Ri = N² / |S|² - ᶜfb = @. ᶜtemp_scalar = ifelse(ᶜRi ≤ 0, 1, max(0, 1 - ᶜRi / Pr_t)^(1 / 4)) + # Precompute 3D strain rate tensor + compute_strain_rate_center_full!(ᶜS, ᶜu, ᶠu) + compute_strain_rate_face_full!(ᶠS, ᶜu, ᶠu) # filter scale h_space = Spaces.horizontal_space(axes(Y.c)) - Δ_xy = Spaces.node_horizontal_length_scale(h_space)^2 # Δ_x * Δ_y - ᶜΔ_z = Fields.Δz_field(Y.c) - ᶜΔ = @. ᶜtemp_scalar = ∛(Δ_xy * ᶜΔ_z) * ᶜfb + Δx = Δy = Spaces.node_horizontal_length_scale(h_space) + ᶜΔz = Fields.Δz_field(Y.c) + ax_xy = is_smagorinsky_UVW_coupled(model) ? Geometry.UVWAxis() : Geometry.UVAxis() + ax_z = is_smagorinsky_UVW_coupled(model) ? Geometry.UVWAxis() : Geometry.WAxis() + + ᶜfb = lilly_stratification_correction(p, ᶜS) + if is_smagorinsky_UVW_coupled(model) + ᶜL_h = ᶜL_v = @. lazy(c_smag * cbrt(Δx * Δy * ᶜΔz) * ᶜfb) + end - # Smagorinsky-Lilly eddy viscosity - ᶜνₜ = @. ᶜtemp_scalar = c_smag^2 * ᶜΔ^2 * ᶜS_norm - ᶠνₜ = @. ᶠtemp_scalar = ᶠinterp(ᶜνₜ) + # Cache strain rate norms for diagnostics + ᶜS_norm_h .= strain_rate_norm(ᶜS, ax_xy) + ᶜS_norm_v .= strain_rate_norm(ᶜS, ax_z) - # Subgrid-scale momentum flux tensor, `τ = -2 νₜ ∘ S` - @. ᶜτ_smag = -2 * ᶜνₜ * ᶜS - @. ᶠτ_smag = -2 * ᶠνₜ * ᶠS + # Smagorinsky eddy viscosity + @. ᶜνₜ_h = ᶜL_h^2 * ᶜS_norm_h + @. ᶜνₜ_v = ᶜL_v^2 * ᶜS_norm_v # Turbulent diffusivity - @. ᶜD_smag = ᶜνₜ / Pr_t - @. ᶠD_smag = ᶠνₜ / Pr_t + Pr_t = CAP.Prandtl_number_0(CAP.turbconv_params(p.params)) + @. ᶜD_h = ᶜνₜ_h / Pr_t + @. ᶜD_v = ᶜνₜ_v / Pr_t nothing end +set_smagorinsky_lilly_precomputed_quantities!(Y, p, ::Nothing) = nothing horizontal_smagorinsky_lilly_tendency!(Yₜ, Y, p, t, ::Nothing) = nothing vertical_smagorinsky_lilly_tendency!(Yₜ, Y, p, t, ::Nothing) = nothing -function horizontal_smagorinsky_lilly_tendency!(Yₜ, Y, p, t, ::SmagorinskyLilly) - (; ᶜτ_smag, ᶠτ_smag, ᶜD_smag, ᶜts) = p.precomputed +function horizontal_smagorinsky_lilly_tendency!(Yₜ, Y, p, t, model::SmagorinskyLilly) + is_smagorinsky_horizontal(model) || return nothing + (; ᶜts, ᶜS, ᶠS, ᶜνₜ_h, ᶜD_h) = p.precomputed + (; ᶜtemp_UVWxUVW, ᶠtemp_UVWxUVW, ᶜtemp_scalar, ᶠtemp_scalar) = p.scratch thermo_params = CAP.thermodynamics_params(p.params) + ᶜρ = Y.c.ρ + ᶠρ = @. ᶠtemp_scalar = ᶠinterp(ᶜρ) + + # Subgrid-scale momentum flux tensor, `τ = -2 νₜ ∘ S` + ᶠνₜ_h = @. lazy(ᶠinterp(ᶜνₜ_h)) + ᶜτ_smag = @. ᶜtemp_UVWxUVW = -2 * ᶜνₜ_h * ᶜS # TODO: Lazify once we can mix lazy horizontal & vertical operations + ᶠτ_smag = @. ᶠtemp_UVWxUVW = -2 * ᶠνₜ_h * ᶠS - ## Momentum tendencies - ᶠρ = @. p.scratch.ᶠtemp_scalar = ᶠinterp(Y.c.ρ) - @. Yₜ.c.uₕ -= C12(wdivₕ(Y.c.ρ * ᶜτ_smag) / Y.c.ρ) + # Apply to tendencies + ## Horizontal momentum tendency + @. Yₜ.c.uₕ -= C12(wdivₕ(ᶜρ * ᶜτ_smag) / ᶜρ) + ## Vertical momentum tendency @. Yₜ.f.u₃ -= C3(wdivₕ(ᶠρ * ᶠτ_smag) / ᶠρ) ## Total energy tendency - ᶜh_tot = @. lazy( - TD.total_specific_enthalpy( - thermo_params, - ᶜts, - specific(Y.c.ρe_tot, Y.c.ρ), - ), - ) - @. Yₜ.c.ρe_tot += wdivₕ(Y.c.ρ * ᶜD_smag * gradₕ(ᶜh_tot)) + ᶜe_tot = @. lazy(specific(Y.c.ρe_tot, ᶜρ)) + ᶜh_tot = @. lazy(TD.total_specific_enthalpy(thermo_params, ᶜts, ᶜe_tot)) + @. Yₜ.c.ρe_tot += wdivₕ(ᶜρ * ᶜD_h * gradₕ(ᶜh_tot)) ## Tracer diffusion and associated mass changes foreach_gs_tracer(Yₜ, Y) do ᶜρχₜ, ᶜρχ, ρχ_name - ᶜχ = @. lazy(specific(ᶜρχ, Y.c.ρ)) - ᶜρχₜ_diffusion = @. lazy(wdivₕ(Y.c.ρ * ᶜD_smag * gradₕ(ᶜχ))) - @. ᶜρχₜ += ᶜρχₜ_diffusion + ᶜχ = @. lazy(specific(ᶜρχ, ᶜρ)) + ᶜ∇ₕρD∇χₜ = @. lazy(wdivₕ(ᶜρ * ᶜD_h * gradₕ(ᶜχ))) + @. ᶜρχₜ += ᶜ∇ₕρD∇χₜ # Rain and snow does not affect the mass if ρχ_name == @name(ρq_tot) - @. Yₜ.c.ρ += ᶜρχₜ_diffusion + @. Yₜ.c.ρ += ᶜ∇ₕρD∇χₜ end end - end -import UnrolledUtilities as UU - -function vertical_smagorinsky_lilly_tendency!(Yₜ, Y, p, t, ::SmagorinskyLilly) +function vertical_smagorinsky_lilly_tendency!(Yₜ, Y, p, t, model::SmagorinskyLilly) + is_smagorinsky_vertical(model) || return nothing FT = eltype(Y) - (; sfc_temp_C3, ᶠtemp_scalar) = p.scratch - (; ᶜτ_smag, ᶠτ_smag, ᶠD_smag, sfc_conditions) = p.precomputed - (; ρ_flux_uₕ, ρ_flux_h_tot) = sfc_conditions - (; ᶜts) = p.precomputed + (; ᶜts, ᶜS, ᶠS, ᶜνₜ_v) = p.precomputed + (; ᶜtemp_UVWxUVW, ᶠtemp_UVWxUVW, ᶠtemp_scalar, ᶠtemp_scalar_2) = p.scratch + Pr_t = CAP.Prandtl_number_0(CAP.turbconv_params(p.params)) thermo_params = CAP.thermodynamics_params(p.params) + ᶜρ = Y.c.ρ + ᶠρ = @. ᶠtemp_scalar = ᶠinterp(ᶜρ) # Define operators ᶠgradᵥ = Operators.GradientC2F() # apply BCs to ᶜdivᵥ, which wraps ᶠgradᵥ - ᶜdivᵥ_uₕ = Operators.DivergenceF2C( - top = Operators.SetValue(C3(FT(0)) ⊗ C12(FT(0), FT(0))), - bottom = Operators.SetValue(C3(FT(0)) ⊗ C12(FT(0), FT(0))), - ) - ᶠdivᵥ = Operators.DivergenceC2F( - bottom = Operators.SetDivergence(FT(0)), - top = Operators.SetDivergence(FT(0)), - ) - ᶜdivᵥ_ρe_tot = Operators.DivergenceF2C(; - top = Operators.SetValue(C3(FT(0))), - bottom = Operators.SetValue(C3(FT(0))), - ) + divᵥ_uₕ_bc = Operators.SetValue(C3(FT(0)) ⊗ C12(FT(0), FT(0))) + ᶜdivᵥ_uₕ = Operators.DivergenceF2C(; bottom = divᵥ_uₕ_bc, top = divᵥ_uₕ_bc) + divᵥ_bc = Operators.SetDivergence(FT(0)) + ᶠdivᵥ = Operators.DivergenceC2F(; bottom = divᵥ_bc, top = divᵥ_bc) + divᵥ_ρχ_bc = Operators.SetValue(C3(FT(0))) + ᶜdivᵥ_ρχ = Operators.DivergenceF2C(; bottom = divᵥ_ρχ_bc, top = divᵥ_ρχ_bc) + + # Subgrid-scale momentum flux tensor, `τ = -2 νₜ ∘ S` + ᶠνₜ_v = @. lazy(ᶠinterp(ᶜνₜ_v)) + ᶜτ_smag = @. ᶜtemp_UVWxUVW = -2 * ᶜνₜ_v * ᶜS + ᶠτ_smag = @. ᶠtemp_UVWxUVW = -2 * ᶠνₜ_v * ᶠS + + # Turbulent diffusivity + ᶠD_smag = @. lazy(ᶠνₜ_v / Pr_t) # Apply to tendencies ## Horizontal momentum tendency - ᶠρ = @. ᶠtemp_scalar = ᶠinterp(Y.c.ρ) - @. Yₜ.c.uₕ -= C12(ᶜdivᵥ(ᶠρ * ᶠτ_smag) / Y.c.ρ) + @. Yₜ.c.uₕ -= C12(ᶜdivᵥ(ᶠρ * ᶠτ_smag) / ᶜρ) ## Apply boundary condition for momentum flux - @. Yₜ.c.uₕ -= ᶜdivᵥ_uₕ(-(FT(0) * ᶠgradᵥ(Y.c.uₕ))) / Y.c.ρ + @. Yₜ.c.uₕ -= ᶜdivᵥ_uₕ(-(FT(0) * ᶠgradᵥ(Y.c.uₕ))) / ᶜρ ## Vertical momentum tendency - @. Yₜ.f.u₃ -= C3(ᶠdivᵥ(Y.c.ρ * ᶜτ_smag) / ᶠρ) + @. Yₜ.f.u₃ -= C3(ᶠdivᵥ(ᶜρ * ᶜτ_smag) / ᶠρ) ## Total energy tendency - ᶜh_tot = @. lazy( - TD.total_specific_enthalpy( - thermo_params, - ᶜts, - specific(Y.c.ρe_tot, Y.c.ρ), - ), - ) - @. Yₜ.c.ρe_tot -= ᶜdivᵥ_ρe_tot(-(ᶠρ * ᶠD_smag * ᶠgradᵥ(ᶜh_tot))) + ᶜe_tot = @. lazy(specific(Y.c.ρe_tot, ᶜρ)) + ᶜh_tot = @. lazy(TD.total_specific_enthalpy(thermo_params, ᶜts, ᶜe_tot)) + @. Yₜ.c.ρe_tot -= ᶜdivᵥ_ρχ(-(ᶠρ * ᶠD_smag * ᶠgradᵥ(ᶜh_tot))) ## Tracer diffusion and associated mass changes - ᶜdivᵥ_ρχ = Operators.DivergenceF2C(; - top = Operators.SetValue(C3(FT(0))), - bottom = Operators.SetValue(C3(FT(0))), - ) - foreach_gs_tracer(Yₜ, Y) do ᶜρχₜ, ᶜρχ, ρχ_name - ᶜ∇ᵥρD∇χₜ = - @. lazy(ᶜdivᵥ_ρχ(-(ᶠρ * ᶠD_smag * ᶠgradᵥ(specific(ᶜρχ, Y.c.ρ))))) + ᶜχ = @. lazy(specific(ᶜρχ, ᶜρ)) + ᶜ∇ᵥρD∇χₜ = @. lazy(ᶜdivᵥ_ρχ(-(ᶠρ * ᶠD_smag * ᶠgradᵥ(ᶜχ)))) @. ᶜρχₜ -= ᶜ∇ᵥρD∇χₜ # Rain and snow does not affect the mass if ρχ_name == @name(ρq_tot) diff --git a/src/prognostic_equations/implicit/manual_sparse_jacobian.jl b/src/prognostic_equations/implicit/manual_sparse_jacobian.jl index 2142bc9b88..cb53a98013 100644 --- a/src/prognostic_equations/implicit/manual_sparse_jacobian.jl +++ b/src/prognostic_equations/implicit/manual_sparse_jacobian.jl @@ -583,7 +583,7 @@ function update_jacobian!(alg::ManualSparseJacobian, cache, Y, p, dtγ, t) (; turbconv_model) = p.atmos turbconv_params = CAP.turbconv_params(params) FT = eltype(params) - (; vertical_diffusion) = p.atmos + (; vertical_diffusion, smagorinsky_lilly) = p.atmos (; ᶜp) = p.precomputed if vertical_diffusion isa DecayWithHeightDiffusion ᶜK_h = ᶜcompute_eddy_diffusivity_coefficient(Y.c.ρ, vertical_diffusion) @@ -591,6 +591,10 @@ function update_jacobian!(alg::ManualSparseJacobian, cache, Y, p, dtγ, t) elseif vertical_diffusion isa VerticalDiffusion ᶜK_h = ᶜcompute_eddy_diffusivity_coefficient(Y.c.uₕ, ᶜp, vertical_diffusion) ᶜK_u = ᶜK_h + elseif is_smagorinsky_vertical(smagorinsky_lilly) + set_smagorinsky_lilly_precomputed_quantities!(Y, p, smagorinsky_lilly) + ᶜK_u = p.precomputed.ᶜνₜ_v + ᶜK_h = p.precomputed.ᶜD_v elseif turbconv_model isa AbstractEDMF (; ᶜlinear_buoygrad, ᶜstrain_rate_norm) = p.precomputed ᶜρa⁰ = diff --git a/src/solver/model_getters.jl b/src/solver/model_getters.jl index 7a76bf840a..ac454ae987 100644 --- a/src/solver/model_getters.jl +++ b/src/solver/model_getters.jl @@ -169,10 +169,20 @@ function get_viscous_sponge_model(parsed_args, params, ::Type{FT}) where {FT} end end +""" + get_smagorinsky_lilly_model(parsed_args) + +Get the Smagorinsky-Lilly turbulence model based on `parsed_args["smagorinsky_lilly"]` + +The possible model configurations flags are: +- `UVW`: Applies the model to all spatial directions. +""" function get_smagorinsky_lilly_model(parsed_args) - is_model_active = parsed_args["smagorinsky_lilly"] - @assert is_model_active in (true, false) - return is_model_active ? SmagorinskyLilly() : nothing + smag = parsed_args["smagorinsky_lilly"] + isnothing(smag) && return nothing + smag_symbol = Symbol(smag) + @assert smag_symbol in (:UVW,) + return SmagorinskyLilly{smag_symbol}() end function get_amd_les_model(parsed_args, ::Type{FT}) where {FT} diff --git a/src/solver/types.jl b/src/solver/types.jl index f4fd58891b..f0f27a644c 100644 --- a/src/solver/types.jl +++ b/src/solver/types.jl @@ -239,7 +239,48 @@ Base.@kwdef struct ViscousSponge{FT} <: AbstractSponge end abstract type AbstractEddyViscosityModel end -struct SmagorinskyLilly <: AbstractEddyViscosityModel end +""" + SmagorinskyLilly{AXES} + +Smagorinsky-Lilly eddy viscosity model. + +`AXES` is a symbol indicating along which axes the model is applied. It can be +- `:UVW` (all axes) +- `:UV` (horizontal axes) +- `:W` (vertical axis) +""" +struct SmagorinskyLilly{AXES} <: AbstractEddyViscosityModel end + +""" + is_smagorinsky_UVW_coupled(model) + +Check if the Smagorinsky model is coupled in all directions. +""" +is_smagorinsky_UVW_coupled(::SmagorinskyLilly{AXES}) where {AXES} = AXES == :UVW +is_smagorinsky_UVW_coupled(::Nothing) = false + +""" + is_smagorinsky_vertical(model) + +Check if the Smagorinsky model is applied in the vertical direction. + +See also [`is_smagorinsky_horizontal`](@ref). +""" +is_smagorinsky_vertical(::SmagorinskyLilly{AXES}) where {AXES} = + AXES == :UVW +is_smagorinsky_vertical(::Nothing) = false + +""" + is_smagorinsky_horizontal(model) + +Check if the Smagorinsky model is applied in the horizontal directions. + +See also [`is_smagorinsky_vertical`](@ref). +""" +is_smagorinsky_horizontal(::SmagorinskyLilly{AXES}) where {AXES} = + AXES == :UVW +is_smagorinsky_horizontal(::Nothing) = false + struct AnisotropicMinimumDissipation{FT} <: AbstractEddyViscosityModel c_amd::FT end