update syntax

Author: haakon-e

src/cache/precomputed_quantities.jl

@@ -253,7 +253,7 @@ function precomputed_quantities(Y, atmos)
                 ᶜS = similar(Y.c, typeof(uvw_vec * uvw_vec')),
                 ᶠS = similar(Y.f, typeof(uvw_vec * uvw_vec')),
                 ᶜL_h = similar(Y.c, FT), ᶜL_v = similar(Y.c, FT),
-                ᶜD_h = similar(Y.c, FT), ᶠD_v = similar(Y.f, FT),
+                ᶜD_h = similar(Y.c, FT), ᶜD_v = similar(Y.c, FT),
                 ᶜstrain_rate_norm_h = similar(Y.c, FT),
                 ᶜstrain_rate_norm_v = similar(Y.c, FT),
             )

src/diagnostics/core_diagnostics.jl

@@ -359,7 +359,7 @@ add_diagnostic_variable!(
     long_name = "Horizontal smagorinski diffusivity",
     units = "m^2 s^-1",
     compute! = (out, _, cache, _) -> begin
-        isnothing(out) ? copy(cache.precomputed.ᶠD_v) : (out .= cache.precomputed.ᶠD_v)
+        isnothing(out) ? copy(cache.precomputed.ᶜD_v) : (out .= cache.precomputed.ᶜD_v)
     end,
 )
 add_diagnostic_variable!(
@@ -378,26 +378,27 @@ add_diagnostic_variable!(
     units = "s",
     compute! = (out, state, cache, _) -> begin
         ᶜΔ_z = Fields.Δz_field(state.c)
-        dtmax = @. lazy(ᶠinterp(ᶜΔ_z^2) / cache.precomputed.ᶠD_v)
+        dtmax = @. lazy(ᶜΔ_z^2 / cache.precomputed.ᶜD_v)
         isnothing(out) ? copy(dtmax) : (out .= dtmax)
     end,
 )
+import ..projected_strain_rate_norm
 add_diagnostic_variable!(
     short_name = "strainrateh_smag",
     long_name = "Strain rate horizontal",
     units = "s",
     compute! = (out, state, cache, _) -> begin
-        ᶜstrain_rate_h = cache.precomputed.ᶜstrain_rate_norm_h
-        isnothing(out) ? copy(ᶜstrain_rate_h) : (out .= ᶜstrain_rate_h)
+        ᶜS_norm_h = projected_strain_rate_norm(cache.precomputed.ᶜS, Geometry.UVAxis())
+        isnothing(out) ? copy(ᶜS_norm_h) : (out .= ᶜS_norm_h)
     end,
 )
 add_diagnostic_variable!(
     short_name = "strainratev_smag",
     long_name = "Strain rate vertical",
     units = "s",
     compute! = (out, state, cache, _) -> begin
-        ᶜstrain_rate_v = cache.precomputed.ᶜstrain_rate_norm_v
-        isnothing(out) ? copy(ᶜstrain_rate_v) : (out .= ᶜstrain_rate_v)
+        ᶜS_norm_v = projected_strain_rate_norm(cache.precomputed.ᶜS, Geometry.WAxis())
+        isnothing(out) ? copy(ᶜS_norm_v) : (out .= ᶜS_norm_v)
     end,
 )
 

src/parameterized_tendencies/les_sgs_models/smagorinsky_lilly.jl

@@ -28,6 +28,7 @@ These quantities are computed for both cell centers and faces, with prefixes `
 function set_smagorinsky_lilly_precomputed_quantities!(Y, p)
     FT = eltype(Y)
     (; ᶜu, ᶠu, ᶜts, ᶜL_h, ᶜL_v, ᶠS, ᶜS) = p.precomputed
+    (; ᶜstrain_rate_norm_h, ᶜstrain_rate_norm_v) = p.precomputed  # diagnostics
     c_smag = CAP.c_smag(p.params)
     grav = CAP.grav(p.params)
     Pr_t = CAP.Prandtl_number_0(CAP.turbconv_params(p.params))
@@ -61,24 +62,33 @@ end
 horizontal_smagorinsky_lilly_tendency!(Yₜ, Y, p, t, ::Nothing) = nothing
 vertical_smagorinsky_lilly_tendency!(Yₜ, Y, p, t, ::Nothing) = nothing
 
-function projected_strain_rate_norm(ᶜS, axis)
-    ᶜS_proj = @. lazy(Geometry.project((axis,), ᶜS, (axis,)))
-    ᶜS_norm = @. lazy(√(2 * norm_sqr(ᶜS_proj)))
-    return ᶜS_norm
+function projected_strain_rate_norm(S, axis)
+    S_proj = @. lazy(Geometry.project((axis,), S, (axis,)))
+    S_norm = @. lazy(√(2 * norm_sqr(S_proj)))
+    return S_norm
+end
+
+function horizontal_smagorinsky_eddy_viscosity(L, S)
+    S_norm = projected_strain_rate_norm(S, Geometry.UVAxis())
+    νₜ = @. lazy(L^2 * S_norm)
+    return νₜ
+end
+
+function vertical_smagorinsky_eddy_viscosity(L, S)
+    S_norm = projected_strain_rate_norm(S, Geometry.WAxis())
+    νₜ = @. lazy(L^2 * S_norm)
+    return νₜ
 end
 
 function horizontal_smagorinsky_lilly_tendency!(Yₜ, Y, p, t, ::SmagorinskyLilly)
     (; ᶜts, ᶜL_h, ᶠS, ᶜS) = p.precomputed
     (; ᶜtemp_UVWxUVW, ᶠtemp_UVWxUVW, ᶜtemp_scalar, ᶠtemp_scalar) = p.scratch
     thermo_params = CAP.thermodynamics_params(p.params)
     Pr_t = CAP.Prandtl_number_0(CAP.turbconv_params(p.params))
-
-    ## Momentum tendencies
-    ᶜS_norm = projected_strain_rate_norm(ᶜS, Geometry.UVAxis())
-    @. p.precomputed.ᶜstrain_rate_norm_h = ᶜS_norm  # save to diagnostics
+    ᶜρ = Y.c.ρ
 
     # Smagorinsky eddy viscosity
-    ᶜνₜ_h = @. lazy(ᶜL_h^2 * ᶜS_norm)
+    ᶜνₜ_h = horizontal_smagorinsky_eddy_viscosity(ᶜL_h, ᶜS)
     ᶠνₜ_h = @. lazy(ᶠinterp(ᶜνₜ_h))
 
     # Turbulent diffusivity
@@ -90,22 +100,22 @@ function horizontal_smagorinsky_lilly_tendency!(Yₜ, Y, p, t, ::SmagorinskyLill
     ᶠτ_smag = @. ᶠtemp_UVWxUVW = -2 * ᶠνₜ_h * ᶠS
 
     ᶠρ = ᶠtemp_scalar .= face_density(Y)
-    @. Yₜ.c.uₕ -= C12(wdivₕ(Y.c.ρ * ᶜτ_smag) / Y.c.ρ)
+    @. Yₜ.c.uₕ -= C12(wdivₕ(ᶜρ * ᶜτ_smag) / ᶜρ)
     @. Yₜ.f.u₃ -= C3(wdivₕ(ᶠρ * ᶠτ_smag) / ᶠρ)
 
     ## Total energy tendency
-    ᶜe_tot = @. lazy(specific(Y.c.ρe_tot, Y.c.ρ))
+    ᶜ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ₕ(Y.c.ρ * ᶜD_smag * gradₕ(ᶜh_tot))
+    @. Yₜ.c.ρe_tot += wdivₕ(ᶜρ * ᶜD_smag * 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_smag * gradₕ(ᶜχ)))
+        @. ᶜρχₜ += ᶜ∇ₕρD∇χₜ
         # Rain and snow does not affect the mass
         if ρχ_name == @name(ρq_tot)
-            @. Yₜ.c.ρ += ᶜρχₜ_diffusion
+            @. Yₜ.c.ρ += ᶜ∇ₕρD∇χₜ
         end
     end
 
@@ -117,6 +127,7 @@ function vertical_smagorinsky_lilly_tendency!(Yₜ, Y, p, t, ::SmagorinskyLilly)
     (; ᶜtemp_UVWxUVW, ᶠtemp_UVWxUVW, ᶠtemp_scalar, ᶠtemp_scalar_2) = p.scratch
     thermo_params = CAP.thermodynamics_params(p.params)
     Pr_t = CAP.Prandtl_number_0(CAP.turbconv_params(p.params))
+    ᶜρ = Y.c.ρ
 
     # Define operators
     ᶠgradᵥ = Operators.GradientC2F() # apply BCs to ᶜdivᵥ, which wraps ᶠgradᵥ
@@ -125,20 +136,14 @@ function vertical_smagorinsky_lilly_tendency!(Yₜ, Y, p, t, ::SmagorinskyLilly)
         bottom = Operators.SetValue(C3(FT(0)) ⊗ C12(FT(0), FT(0))),
     )
     ᶠdivᵥ = Operators.DivergenceC2F(
-        bottom = Operators.SetDivergence(FT(0)),
-        top = Operators.SetDivergence(FT(0)),
+        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ᵥ_ρχ = Operators.DivergenceF2C(;
+        top = Operators.SetValue(C3(FT(0))), bottom = Operators.SetValue(C3(FT(0))),
     )
 
-    ## Momentum tendencies
-    ᶜS_norm = projected_strain_rate_norm(ᶜS, Geometry.WAxis())
-    @. p.precomputed.ᶜstrain_rate_norm_v = ᶜS_norm  # save to diagnostics
-
     # Smagorinsky eddy viscosity
-    ᶜνₜ_v = @. lazy(ᶜL_v^2 * ᶜS_norm)
+    ᶜνₜ_v = vertical_smagorinsky_eddy_viscosity(ᶜL_v, ᶜS)
     ᶠνₜ_v = @. lazy(ᶠinterp(ᶜνₜ_v))
 
     # Subgrid-scale momentum flux tensor, `τ = -2 νₜ ∘ S`
@@ -147,30 +152,26 @@ function vertical_smagorinsky_lilly_tendency!(Yₜ, Y, p, t, ::SmagorinskyLilly)
 
     # Turbulent diffusivity
     ᶠD_smag = @. ᶠtemp_scalar_2 = ᶠνₜ_v / Pr_t
-    @. p.precomputed.ᶠD_v = ᶠD_smag  # save to diagnostics
+    @. p.precomputed.ᶜD_v = ᶜνₜ_v / Pr_t  # save to diagnostics, on centers
 
     # Apply to tendencies
     ## Horizontal momentum tendency
     ᶠρ = ᶠtemp_scalar .= face_density(Y)
-    @. 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
-    ᶜe_tot = @. lazy(specific(Y.c.ρe_tot, Y.c.ρ))
+    ᶜ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ᵥ_ρe_tot(-(ᶠρ * ᶠD_smag * ᶠgradᵥ(ᶜh_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)

src/prognostic_equations/implicit/manual_sparse_jacobian.jl

@@ -594,9 +594,8 @@ function update_jacobian!(alg::ManualSparseJacobian, cache, Y, p, dtγ, t)
         elseif p.atmos.smagorinsky_vertical isa SmagorinskyLilly
             Pr_t = CAP.Prandtl_number_0(CAP.turbconv_params(p.params))
             set_smagorinsky_lilly_precomputed_quantities!(Y, p)  # sets (ᶜS, ᶜL_v) in p.precomputed
-            (; ᶜS, ᶜL_v) = p.precomputed
-            ᶜS_norm = projected_strain_rate_norm(ᶜS, Geometry.WAxis())
-            ᶜK_u = @. lazy(ᶜL_v^2 * ᶜS_norm)  # Smagorinsky eddy viscosity, ᶜνₜ_v
+            (; ᶜL_v, ᶜS) = p.precomputed
+            ᶜK_u = vertical_smagorinsky_eddy_viscosity(ᶜL_v, ᶜS)  # Smagorinsky eddy viscosity, ᶜνₜ_v
             ᶜK_h = @. lazy(ᶜK_u / Pr_t)  # Smagorinsky eddy diffusivity, ᶜD_smag
         elseif turbconv_model isa AbstractEDMF
             (; ᶜlinear_buoygrad, ᶜstrain_rate_norm) = p.precomputed