fix: compute full strain

Author: haakon-e

src/cache/precomputed_quantities.jl

@@ -248,7 +248,12 @@ function precomputed_quantities(Y, atmos)
     smagorinsky_quantities =
         if atmos.smagorinsky_horizontal isa SmagorinskyLilly ||
            atmos.smagorinsky_vertical isa SmagorinskyLilly
-            (; ᶜL_h = similar(Y.c, FT), ᶜL_v = similar(Y.c, FT))
+           uvw_vec = UVW(FT(0), FT(0), FT(0))
+            (;
+                ᶜ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),
+            )
         else
             (;)
         end

src/parameterized_tendencies/les_sgs_models/smagorinsky_lilly.jl

@@ -26,16 +26,17 @@ These quantities are computed for both cell centers and faces, with prefixes `
 - `p`: The model parameters, e.g. `AtmosCache`.
 """
 function set_smagorinsky_lilly_precomputed_quantities!(Y, p)
-    FT = eltype(Y)
-    (; ᶠu, ᶜts, ᶜL_h, ᶜL_v) = p.precomputed
+    # FT = eltype(Y)
+    (; ᶜu, ᶠu, ᶜts, ᶜL_h, ᶜL_v, ᶠS, ᶜS) = p.precomputed
     c_smag = CAP.c_smag(p.params)
     # grav = CAP.grav(p.params)
     # Pr_t = CAP.Prandtl_number_0(CAP.turbconv_params(p.params))
     # thermo_params = CAP.thermodynamics_params(p.params)
     # (; ᶜtemp_scalar, ᶜtemp_scalar_2) = p.scratch
 
-    # Strain rate tensor
-    # ᶜS = compute_strain_rate_center(ᶠu)
+    # Precompute full strain rate tensor
+    compute_strain_rate_center_full!(ᶜS, ᶜu, ᶠu)
+    compute_strain_rate_face_full!(ᶠS, ᶜu, ᶠu)
 
     # Stratification correction
     # ᶜθ_v = @. lazy(TD.virtual_pottemp(thermo_params, ᶜts))
@@ -61,14 +62,12 @@ 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)
-    (; ᶜu, ᶠu, ᶜts, ᶜL_h) = p.precomputed
+    (; ᶜ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 = compute_strain_rate_center(ᶠu)
-    ᶠS = compute_strain_rate_face(ᶜu)
     ᶜS_h = @. lazy(Geometry.project((Geometry.UVAxis(),), ᶜS, (Geometry.UVAxis(),)))
     ᶜS_norm = @. lazy(√(2 * norm_sqr(ᶜS_h)))
 
@@ -88,9 +87,8 @@ function horizontal_smagorinsky_lilly_tendency!(Yₜ, Y, p, t, ::SmagorinskyLill
     @. 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.ρ)),
-    )
+    ᶜe_tot = @. lazy(specific(Y.c.ρe_tot, Y.c.ρ))
+    ᶜh_tot = @. lazy(TD.total_specific_enthalpy(thermo_params, ᶜts, ᶜe_tot))
     @. Yₜ.c.ρe_tot += wdivₕ(Y.c.ρ * ᶜD_smag * gradₕ(ᶜh_tot))
 
     ## Tracer diffusion and associated mass changes
@@ -108,7 +106,7 @@ end
 
 function vertical_smagorinsky_lilly_tendency!(Yₜ, Y, p, t, ::SmagorinskyLilly)
     FT = eltype(Y)
-    (; ᶜu, ᶠu, ᶜts, ᶜL_v) = p.precomputed
+    (; ᶜts, ᶜL_v, ᶠS, ᶜS) = p.precomputed
     (; ᶜ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))
@@ -129,8 +127,6 @@ function vertical_smagorinsky_lilly_tendency!(Yₜ, Y, p, t, ::SmagorinskyLilly)
     )
 
     ## Momentum tendencies
-    ᶜS = compute_strain_rate_center(ᶠu)
-    ᶠS = compute_strain_rate_face(ᶜu)
     ᶜS_v = @. lazy(Geometry.project((Geometry.WAxis(),), ᶜS, (Geometry.WAxis(),)))
     ᶜS_norm = @. lazy(√(2 * norm_sqr(ᶜS_v)))
 
@@ -155,9 +151,8 @@ function vertical_smagorinsky_lilly_tendency!(Yₜ, Y, p, t, ::SmagorinskyLilly)
     @. Yₜ.f.u₃ -= C3(ᶠdivᵥ(Y.c.ρ * ᶜτ_smag) / ᶠρ)
 
     ## Total energy tendency
-    ᶜh_tot = @. lazy(
-        TD.total_specific_enthalpy(thermo_params, ᶜts, specific(Y.c.ρe_tot, Y.c.ρ)),
-    )
+    ᶜe_tot = @. lazy(specific(Y.c.ρe_tot, Y.c.ρ))
+    ᶜh_tot = @. lazy(TD.total_specific_enthalpy(thermo_params, ᶜts, ᶜe_tot))
     @. Yₜ.c.ρe_tot -= ᶜdivᵥ_ρe_tot(-(ᶠρ * ᶠD_smag * ᶠgradᵥ(ᶜh_tot)))
 
     ## Tracer diffusion and associated mass changes

src/utils/utilities.jl

@@ -84,7 +84,7 @@ compute_kinetic(Y::Fields.FieldVector) = compute_kinetic(Y.c.uₕ, Y.f.u₃)
 """
     ϵ .= compute_strain_rate_center(u::Field)
 
-Compute the strain_rate at cell centers from velocity at cell faces.
+Compute the vertical strain_rate at cell centers from velocity at cell faces.
 """
 function compute_strain_rate_center(u::Fields.Field)
     axis_uvw = Geometry.UVWAxis()
@@ -99,7 +99,7 @@ end
 """
     ϵ .= compute_strain_rate_face(u::Field)
 
-Compute the strain_rate at cell faces from velocity at cell centers.
+Compute the vertical strain_rate at cell faces from velocity at cell centers.
 """
 function compute_strain_rate_face(u::Fields.Field)
     @assert eltype(u) <: C123
@@ -118,6 +118,27 @@ function compute_strain_rate_face(u::Fields.Field)
     )
 end
 
+function compute_strain_rate_center_full!(ε, ᶜu, ᶠu)
+    axis_uvw = (Geometry.UVWAxis(),)
+    @. ε = Geometry.project(axis_uvw, ᶜgradᵥ(UVW(ᶠu)))  # vertical component
+    @. ε += Geometry.project(axis_uvw, gradₕ(UVW(ᶜu)))  # horizontal component
+    @. ε = (ε + adjoint(ε)) / 2
+    return ε
+end
+
+function compute_strain_rate_face_full!(ε, ᶜu, ᶠu)
+    ∇ᵥuvw_boundary = Geometry.outer(Geometry.WVector(0), UVW(0, 0, 0))
+    ᶠgradᵥ = Operators.GradientC2F(
+        bottom = Operators.SetGradient(∇ᵥuvw_boundary),
+        top = Operators.SetGradient(∇ᵥuvw_boundary),
+    )
+    axis_uvw = (Geometry.UVWAxis(),)
+    @. ε = Geometry.project(axis_uvw, ᶠgradᵥ(UVW(ᶜu)))  # vertical component
+    @. ε += Geometry.project(axis_uvw, gradₕ(UVW(ᶠu)))  # horizontal component
+    @. ε = (ε + adjoint(ε)) / 2
+    return ε
+end
+
 """
     face_density(Y)