implicit vertical smagorinsky

Author: haakon-e

post_processing/ci_plots.jl

@@ -1321,6 +1321,7 @@ function make_plots(
         "cl",
         "clw",
         "cli",
+        "Dh_smag", "Dv_smag", "dtmaxh_smag", "dtmaxv_smag", "strainrateh_smag", "strainratev_smag"
     ]
     short_names = short_names ∩ collect(keys(simdirs[1].vars))
 

src/cache/precomputed_quantities.jl

@@ -254,6 +254,8 @@ function precomputed_quantities(Y, atmos)
                 ᶠ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),
+                ᶜstrain_rate_norm_h = similar(Y.c, FT),
+                ᶜstrain_rate_norm_v = similar(Y.c, FT),
             )
         else
             (;)

src/diagnostics/core_diagnostics.jl

@@ -347,21 +347,59 @@ add_diagnostic_variable!(
 # Smagorinski Lilly diffusivity
 ###
 add_diagnostic_variable!(
-    short_name = "D_smag_h",
+    short_name = "Dh_smag",
     long_name = "Horizontal smagorinski diffusivity",
     units = "m^2 s^-1",
     compute! = (out, _, cache, _) -> begin
         isnothing(out) ? copy(cache.precomputed.ᶜD_h) : (out .= cache.precomputed.ᶜD_h)
     end,
 )
 add_diagnostic_variable!(
-    short_name = "D_smag_v",
+    short_name = "Dv_smag",
     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)
     end,
 )
+add_diagnostic_variable!(
+    short_name = "dtmaxh_smag",
+    long_name = "Max horizontal diffusion dt",
+    units = "s",
+    compute! = (out, state, cache, _) -> begin
+        Δh = Spaces.node_horizontal_length_scale(Spaces.horizontal_space(axes(state.f)))
+        dtmax = @. lazy(Δh^2 / cache.precomputed.ᶜD_h)
+        isnothing(out) ? copy(dtmax) : (out .= dtmax)
+    end,
+)
+add_diagnostic_variable!(
+    short_name = "dtmaxv_smag",
+    long_name = "Max diffusion dt",
+    units = "s",
+    compute! = (out, state, cache, _) -> begin
+        ᶜΔ_z = Fields.Δz_field(state.c)
+        dtmax = @. lazy(ᶠinterp(ᶜΔ_z^2) / cache.precomputed.ᶠD_v)
+        isnothing(out) ? copy(dtmax) : (out .= dtmax)
+    end,
+)
+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)
+    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)
+    end,
+)
 
 ###
 # Relative humidity (3d)

src/parameterized_tendencies/les_sgs_models/smagorinsky_lilly.jl

@@ -26,33 +26,33 @@ 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)
+    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
+    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
 
     # 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))
-    # ᶜ∇ᵥθ = @. ᶜtemp_scalar_2 = Geometry.WVector(ᶜgradᵥ(ᶠinterp(ᶜθ_v))).components.data.:1
-    # ᶜN² = @. ᶜtemp_scalar = grav / ᶜθ_v * ᶜ∇ᵥθ
-    # ᶜS_norm = @. ᶜtemp_scalar_2 = √(2 * norm_sqr(ᶜS))
+    ᶜθ_v = @. lazy(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 * norm_sqr(ᶜS))
 
-    # ᶜ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))
+    ᶜ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))
 
     # filter scale
     h_space = Spaces.horizontal_space(axes(Y.c))
     Δ_h = Spaces.node_horizontal_length_scale(h_space)
     ᶜΔ_z = Fields.Δz_field(Y.c)
 
-    @. ᶜL_v = c_smag * ᶜΔ_z #* ᶜfb
+    @. ᶜL_v = c_smag * ᶜΔ_z * ᶜfb
     @. ᶜL_h = c_smag * Δ_h
 
     nothing
@@ -61,15 +61,21 @@ 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
+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_h = @. lazy(Geometry.project((Geometry.UVAxis(),), ᶜS, (Geometry.UVAxis(),)))
-    ᶜS_norm = @. lazy(√(2 * norm_sqr(ᶜS_h)))
+    ᶜS_norm = projected_strain_rate_norm(ᶜS, Geometry.UVAxis())
+    @. p.precomputed.ᶜstrain_rate_norm_h = ᶜS_norm  # save to diagnostics
 
     # Smagorinsky eddy viscosity
     ᶜνₜ_h = @. lazy(ᶜL_h^2 * ᶜS_norm)
@@ -128,8 +134,8 @@ function vertical_smagorinsky_lilly_tendency!(Yₜ, Y, p, t, ::SmagorinskyLilly)
     )
 
     ## Momentum tendencies
-    ᶜS_v = @. lazy(Geometry.project((Geometry.WAxis(),), ᶜS, (Geometry.WAxis(),)))
-    ᶜS_norm = @. lazy(√(2 * norm_sqr(ᶜS_v)))
+    ᶜ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)

src/prognostic_equations/implicit/manual_sparse_jacobian.jl

@@ -591,6 +591,13 @@ 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 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
+            ᶜK_h = @. lazy(ᶜK_u / Pr_t)  # Smagorinsky eddy diffusivity, ᶜD_smag
         elseif turbconv_model isa AbstractEDMF
             (; ᶜlinear_buoygrad, ᶜstrain_rate_norm) = p.precomputed
             ᶜρa⁰ =