modified: config/model_configs/les_dycoms_box.yml

	modified:   src/parameterized_tendencies/les_sgs_models/smagorinsky_lilly.jl
	modified:   test/parameterized_tendencies/les_closures/smagorinsky_lilly.jl

config/model_configs/les_dycoms_box.yml

@@ -17,11 +17,13 @@ y_elem: 4
 z_elem: 300
 z_stretch: false
 smagorinsky_lilly: true
+tracer_upwinding: third_order
 c_smag: 0.15
 dt: "0.05secs"
 t_end: "14400secs"
 dt_save_state_to_disk: "30mins"
 restart_file: "/home/asridhar/Codes/ClimaAtmos.jl/restart/dycoms/day0_1m.hdf5"
+#restart_file: "/home/asridhar/Codes/ClimaAtmos.jl/restart/dycoms/resolution_10m_day0.0.hdf5"
 dt_cloud_fraction: "10mins"
 rayleigh_sponge: true
 toml: [toml/dycoms_box_rhoe.toml]

src/parameterized_tendencies/les_sgs_models/smagorinsky_lilly.jl

@@ -42,6 +42,7 @@ function horizontal_smagorinsky_lilly_tendency!(Yₜ, Y, p, t, sl::SmagorinskyLi
     FT = eltype(v_t)
     grav = CAP.grav(params)
     ᶜJ = Fields.local_geometry_field(Y.c).J 
+    Pr_t_neutral = FT(1/3)
 
     ᶜS  = zero(p.scratch.ᶜtemp_strain)
     ᶠS  = zero(p.scratch.ᶠtemp_strain)
@@ -68,17 +69,18 @@ function horizontal_smagorinsky_lilly_tendency!(Yₜ, Y, p, t, sl::SmagorinskyLi
     ᶜ∇ = Operators.GradientF2C();
     θc2f = Operators.InterpolateC2F(;top = Operators.Extrapolate(),
                                      bottom = Operators.SetValue(θ_v_sfc));
-    ᶜ∇θ = @. ᶜ∇(θc2f(θ_v))
+    ᶜ∇θ = @. ᶜ∇(ᶠinterp(θ_v))
 
+    ᶜfb .= zero(ᶜfb) .+ FT(1)
     N² = @. grav / θ_v * Geometry.WVector(ᶜ∇θ).components.data.:1
     @. ᶜfb = (max(FT(0), 
-                  1 - 3*(N²) / (CA.norm_sqr(ᶜS) + eps(FT))))^(1/2)
+                  1 - (N²) / Pr_t_neutral / (CA.norm_sqr(ᶜS) + eps(FT))))^(1/2)
     ᶜfb .= ifelse.(N² .<= FT(0), zero(ᶜfb) .+ FT(1),ᶜfb)
 
     ᶠρ = @. ᶠinterp(Y.c.ρ)
     ᶜv_t = @. (Cs * Δ_filter)^2 * sqrt(2 * CA.norm_sqr(ᶜS)) * ᶜfb
     ᶠv_t = @. ᶠinterp(ᶜv_t)
-    ᶜD = @. FT(3) * ᶜv_t
+    ᶜD = @. ᶜv_t / Pr_t_neutral
 
     ᶜτ = @. FT(-2) * ᶜv_t * ᶜS
     ᶠτ = @. FT(-2) * ᶠv_t * ᶠS
@@ -115,9 +117,10 @@ function vertical_smagorinsky_lilly_tendency!(Yₜ, Y, p, t, sl::SmagorinskyLill
     (; ᶜu, ᶠu³, sfc_conditions, ᶜspecific, ᶜts) = p.precomputed 
     ᶠgradᵥ = Operators.GradientC2F() # apply BCs to ᶜdivᵥ, which wraps ᶠgradᵥ
 
+    FT = eltype(v_t)
     Δ_filter = Δᵥ_filter
+    Pr_t_neutral = FT(1/3)
 
-    FT = eltype(v_t)
     grav = CAP.grav(params)
     ᶜJ = Fields.local_geometry_field(Y.c).J 
 
@@ -139,22 +142,21 @@ function vertical_smagorinsky_lilly_tendency!(Yₜ, Y, p, t, sl::SmagorinskyLill
 
     thermo_params = CAP.thermodynamics_params(p.params)
     θ_v = @. TD.virtual_pottemp(thermo_params, ᶜts)
-    ### Interp Checks
     ᶠts_sfc = sfc_conditions.ts
     θ_v_sfc = @. TD.virtual_pottemp(thermo_params, ᶠts_sfc)
     ᶜ∇ = Operators.GradientF2C();
     θc2f = Operators.InterpolateC2F(;top = Operators.Extrapolate(),
                                      bottom = Operators.SetValue(θ_v_sfc));
-    ᶜ∇θ = @. ᶜ∇(θc2f(θ_v))
+    ᶜ∇θ = @. ᶜ∇(ᶠinterp(θ_v))
     N² = @. grav / θ_v * Geometry.WVector(ᶜ∇θ).components.data.:1
+    ᶜfb .= zero(ᶜfb) .+ FT(1)
     @. ᶜfb = (max(FT(0), 
-                  1 - 3*(N²) / (CA.norm_sqr(ᶜS) + eps(FT))))^(1/2)
+                  1 - (N²) / Pr_t_neutral / (CA.norm_sqr(ᶜS) + eps(FT))))^(1/2)
     ᶜfb .= ifelse.(N² .<= FT(0), zero(ᶜfb) .+ FT(1),ᶜfb)
 
-    ᶠρ = @. ᶠinterp(Y.c.ρ)
     ᶜv_t = @. (Cs * Δ_filter)^2 * sqrt(2 * CA.norm_sqr(ᶜS)) * ᶜfb
     ᶠv_t = @. ᶠinterp(ᶜv_t)
-    ᶜD = @. FT(3) * ᶜv_t
+    ᶜD = @. ᶜv_t / Pr_t_neutral
 
     @. Yₜ.c.uₕ -= C12(
         ᶜdivᵥ(

test/parameterized_tendencies/les_closures/smagorinsky_lilly.jl

@@ -26,6 +26,7 @@ import Thermodynamics as TD
     simulation = CA.get_simulation(config);
     (; integrator) = simulation;
     Y = integrator.u;
+    Yₜ = similar(Y);
     p = integrator.p;
     params = p.params;
     cm_params = CAP.microphysics_params(params);
@@ -37,11 +38,11 @@ import Thermodynamics as TD
     f_xyz = CC.Fields.coordinate_field(Y.f)
     x = c_xyz.x
     y = c_xyz.y
-    z = c_xyz.z
+    z = f_xyz.z
 
     u = @. sin(x) * cos(y) 
     v = @. sin(y) * cos(x)
-    w = @. z
+    w = @. Geometry.WVector(z ./ maximum(z))
 
     u_x = @. cos(x) * cos(y)
     v_x = @. -sin(y) * sin(x)
@@ -51,23 +52,12 @@ import Thermodynamics as TD
     w_y =  zeros(axes(Y.c))
 
     vel = @. CC.Geometry.UVVector(u, v)
-    vel = @. CC.Geometry.UVWVector(u,v,w)
 
-    wgrad = CC.Operators.WeakGradient()
-    grad = CC.Operators.Gradient()
-    div = CC.Operators.Divergence()
-    wdiv = CC.Operators.WeakDivergence()
-
-    ∇vel = @. grad(vel)
-    ∇velᵀ = similar(∇vel)
-    @. ∇velᵀ = CC.Geometry.AxisTensor(CC.Geometry.axes(∇vel), transpose(CC.Geometry.components(∇vel)))
-    S = @. FT(1/2) * (∇vel + ∇velᵀ)
-    S₃ = @. CC.Geometry.project(CC.Geometry.UVWAxis(), S)
-    # Then, if ᶜϵ is known, we can compute 
-    # the S₃ complete strain-rate tensor, and then take its 
-    # divergence such that we can add it to the tendency terms.
-    divS = @. wdiv(S₃)
+    Y.c.ρ .= FT(1)
+    Y.c.uₕ .= Geometry.Covariant12Vector(vel)
+    Y.f.u₃ .= Geometry.Covariant3Vector(w)
 
+    horizontal_smagorinsky_lilly_tendency(Yₜ, Y, p, t, SmagorinskyLilly(FT(0.2)))
 
     ### Component test begins here
 end