Lagrangian Smagorinsky for 2D flows

src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_model.jl

@@ -179,7 +179,7 @@ function HydrostaticFreeSurfaceModel(; grid,
     closure = with_tracers(tracernames(tracers), closure)
 
     # Put CATKE first in the list of closures
-    closure = validate_closure(closure)
+    closure = validate_closure(closure, grid)
 
     # Either check grid-correctness, or construct tuples of fields
     velocities         = hydrostatic_velocity_fields(velocities, grid, clock, boundary_conditions)

src/Models/NonhydrostaticModels/nonhydrostatic_model.jl

@@ -162,7 +162,7 @@ function NonhydrostaticModel(; grid,
         throw(ArgumentError("hydrostatic_pressure_anomaly must be `nothing` or `CenterField(grid)`."))
 
     # We don't support CAKTE for NonhydrostaticModel yet.
-    closure = validate_closure(closure)
+    closure = validate_closure(closure, grid)
     first_closure = closure isa Tuple ? first(closure) : closure
     first_closure isa FlavorOfCATKE &&
         error("CATKEVerticalDiffusivity is not supported for NonhydrostaticModel --- yet!")

src/TurbulenceClosures/TurbulenceClosures.jl

@@ -73,7 +73,7 @@ Abstract supertype for turbulence closures.
 abstract type AbstractTurbulenceClosure{TimeDiscretization, RequiredHalo} end
 
 # Fallbacks
-validate_closure(closure) = closure
+validate_closure(closure, grid) = closure
 closure_summary(closure) = summary(closure)
 with_tracers(tracers, closure::AbstractTurbulenceClosure) = closure
 compute_diffusivities!(K, closure::AbstractTurbulenceClosure, args...; kwargs...) = nothing

src/TurbulenceClosures/turbulence_closure_implementations/Smagorinskys/dynamic_coefficient.jl

@@ -157,6 +157,9 @@ end
 
 const c = Center()
 
+@inline displace_node(node, δ) = node - δ
+@inline displace_node(::Nothing, δ) = zero(δ)
+
 @kernel function _lagrangian_average_LM_MM!(𝒥ᴸᴹ, 𝒥ᴹᴹ, 𝒥ᴸᴹ⁻, 𝒥ᴹᴹ⁻, 𝒥ᴸᴹ_min, Σ, Σ̄, grid, Δt, u, v, w)
     i, j, k = @index(Global, NTuple)
     LM, MM = LM_and_MM(i, j, k, grid, Σ, Σ̄, u, v, w)
@@ -182,7 +185,7 @@ const c = Center()
         δx = u[i, j, k] * Δt
         δy = v[i, j, k] * Δt
         δz = w[i, j, k] * Δt
-
+        
         # Prevent displacements from getting too big?
         Δx = Δxᶜᶜᶜ(i, j, k, grid)
         Δy = Δyᶜᶜᶜ(i, j, k, grid)
@@ -193,9 +196,9 @@ const c = Center()
         δz = clamp(δz, -Δz, Δz)
 
         # Previous locations
-        x⁻ = x - δx
-        y⁻ = y - δy
-        z⁻ = z - δz
+        x⁻ = displace_node(x, δx)
+        y⁻ = displace_node(y, δy)
+        z⁻ = displace_node(z, δz)
         X⁻ = (x⁻, y⁻, z⁻)
 
         itp_𝒥ᴹᴹ⁻ = interpolate(X⁻, 𝒥ᴹᴹ⁻, (c, c, c), grid)

src/TurbulenceClosures/turbulence_closure_implementations/Smagorinskys/smagorinsky.jl

@@ -1,3 +1,4 @@
+using Oceananigans.Grids: topology
 using Oceananigans.AbstractOperations: Average
 using Oceananigans.Fields: FieldBoundaryConditions
 using Oceananigans.Utils: launch!, IterationInterval
@@ -20,7 +21,8 @@ import ..TurbulenceClosures:
     κᶜᶜᶠ,
     compute_diffusivities!,
     DiffusivityFields,
-    tracer_diffusivities
+    tracer_diffusivities,
+    validate_closure
 
 #####
 ##### The turbulence closure proposed by Smagorinsky and Lilly.
@@ -80,6 +82,15 @@ end
 
 Smagorinsky(FT::DataType; kwargs...) = Smagorinsky(ExplicitTimeDiscretization(), FT; kwargs...)
 
+# Smagorinsky is not correct for 2D Grids
+function validate_closure(closure::Smagorinsky, grid) 
+    if Flat in topology(grid)
+        @warn "A smagorinsky closure relies on assumptions that do not hold for two-dimensional Grids. \n 
+               Consider using a different closure."
+    end
+    return closure
+end
+
 function with_tracers(tracers, closure::Smagorinsky{TD}) where TD
     Pr = tracer_diffusivities(tracers, closure.Pr)
     return Smagorinsky{TD}(closure.coefficient, Pr)

src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/catke_vertical_diffusivity.jl

@@ -127,13 +127,18 @@ function with_tracers(tracer_names, closure::FlavorOfCATKE)
     return closure
 end
 
+
 # For tuples of closures, we need to know _which_ closure is CATKE.
-# Here we take a "simple" approach that sorts the tuple so CATKE is first.
+# Here we take the approach that we validate each individual closure first, then 
+# we sort the tuple so that the CATKE closure is first.
 # This is not sustainable though if multiple closures require this.
 # The two other possibilities are:
 # 1. Recursion to find which closure is CATKE in a compiler-inferrable way
 # 2. Store the "CATKE index" inside CATKE via validate_closure.
-validate_closure(closure_tuple::Tuple) = Tuple(sort(collect(closure_tuple), lt=catke_first))
+function validate_closure(closure_tuple::Tuple, grid) 
+    closure_tuple = Tuple(validate_closure(closure, grid) for closure in closure_tuple)
+    return Tuple(sort(collect(closure_tuple), lt=catke_first))
+end
 
 catke_first(closure1, catke::FlavorOfCATKE) = false
 catke_first(catke::FlavorOfCATKE, closure2) = true