Fix clock argument in TriadIsopycnalSkewSymmetricDiffusivity's explicit_R₃₃_∂z_c function (#4780)

Co-authored-by: Simone Silvestri <[email protected]>
Co-authored-by: Navid C. Constantinou <[email protected]>

Author: 3 people

src/TurbulenceClosures/abstract_scalar_diffusivity_closure.jl

@@ -312,6 +312,7 @@ end
 @inline νᶜᶠᶠ(i, j, k, grid, loc, ν::Number, clk, fields) = ν
 @inline νᶠᶠᶜ(i, j, k, grid, loc, ν::Number, clk, fields) = ν
 
+@inline κᶜᶜᶜ(i, j, k, grid, loc, κ::Number, clk, fields) = κ
 @inline κᶠᶜᶜ(i, j, k, grid, loc, κ::Number, clk, fields) = κ
 @inline κᶜᶠᶜ(i, j, k, grid, loc, κ::Number, clk, fields) = κ
 @inline κᶜᶜᶠ(i, j, k, grid, loc, κ::Number, clk, fields) = κ

src/TurbulenceClosures/turbulence_closure_implementations/isopycnal_skew_symmetric_diffusivity_with_triads.jl

@@ -88,6 +88,10 @@ function DiffusivityFields(grid, tracer_names, bcs, ::FlavorOfTISSD{TD}) where T
     return K
 end
 
+# Build closure fields for model initialization
+build_closure_fields(grid, clock, tracer_names, bcs, closure::FlavorOfTISSD) =
+    DiffusivityFields(grid, tracer_names, bcs, closure)
+
 function compute_diffusivities!(diffusivities, closure::FlavorOfTISSD{TD}, model; parameters = :xyz) where TD
 
     arch = model.architecture
@@ -164,15 +168,15 @@ end
 @inline triad_mask_y(i, jy, jz, ky, kz, grid) = 
    !peripheral_node(i, jy, ky, grid, Center(), Face(), Center()) & !peripheral_node(i, jz, kz, grid, Center(), Center(), Face())
 
-@inline ϵκx⁺⁺(i, j, k, grid, loc, κ, clock, sl, b, C) = triad_mask_x(i+1, i, j, k, k+1, grid) * κᶜᶜᶜ(i, j, k, grid, loc, κ, clock) * tapering_factorᶜᶜᶜ(i, j, k, grid, sl, b, C)
-@inline ϵκx⁺⁻(i, j, k, grid, loc, κ, clock, sl, b, C) = triad_mask_x(i+1, i, j, k, k,   grid) * κᶜᶜᶜ(i, j, k, grid, loc, κ, clock) * tapering_factorᶜᶜᶜ(i, j, k, grid, sl, b, C)
-@inline ϵκx⁻⁺(i, j, k, grid, loc, κ, clock, sl, b, C) = triad_mask_x(i,   i, j, k, k+1, grid) * κᶜᶜᶜ(i, j, k, grid, loc, κ, clock) * tapering_factorᶜᶜᶜ(i, j, k, grid, sl, b, C)
-@inline ϵκx⁻⁻(i, j, k, grid, loc, κ, clock, sl, b, C) = triad_mask_x(i,   i, j, k, k,   grid) * κᶜᶜᶜ(i, j, k, grid, loc, κ, clock) * tapering_factorᶜᶜᶜ(i, j, k, grid, sl, b, C)
+@inline ϵκx⁺⁺(i, j, k, grid, loc, κ, clock, sl, b, C) = triad_mask_x(i+1, i, j, k, k+1, grid) * κᶜᶜᶜ(i, j, k, grid, loc, κ, clock, C) * tapering_factorᶜᶜᶜ(i, j, k, grid, sl, b, C)
+@inline ϵκx⁺⁻(i, j, k, grid, loc, κ, clock, sl, b, C) = triad_mask_x(i+1, i, j, k, k,   grid) * κᶜᶜᶜ(i, j, k, grid, loc, κ, clock, C) * tapering_factorᶜᶜᶜ(i, j, k, grid, sl, b, C)
+@inline ϵκx⁻⁺(i, j, k, grid, loc, κ, clock, sl, b, C) = triad_mask_x(i,   i, j, k, k+1, grid) * κᶜᶜᶜ(i, j, k, grid, loc, κ, clock, C) * tapering_factorᶜᶜᶜ(i, j, k, grid, sl, b, C)
+@inline ϵκx⁻⁻(i, j, k, grid, loc, κ, clock, sl, b, C) = triad_mask_x(i,   i, j, k, k,   grid) * κᶜᶜᶜ(i, j, k, grid, loc, κ, clock, C) * tapering_factorᶜᶜᶜ(i, j, k, grid, sl, b, C)
 
-@inline ϵκy⁺⁺(i, j, k, grid, loc, κ, clock, sl, b, C) = triad_mask_y(i, j+1, j, k, k+1, grid) * κᶜᶜᶜ(i, j, k, grid, loc, κ, clock) * tapering_factorᶜᶜᶜ(i, j, k, grid, sl, b, C)
-@inline ϵκy⁺⁻(i, j, k, grid, loc, κ, clock, sl, b, C) = triad_mask_y(i, j+1, j, k, k,   grid) * κᶜᶜᶜ(i, j, k, grid, loc, κ, clock) * tapering_factorᶜᶜᶜ(i, j, k, grid, sl, b, C)
-@inline ϵκy⁻⁺(i, j, k, grid, loc, κ, clock, sl, b, C) = triad_mask_y(i, j,   j, k, k+1, grid) * κᶜᶜᶜ(i, j, k, grid, loc, κ, clock) * tapering_factorᶜᶜᶜ(i, j, k, grid, sl, b, C)
-@inline ϵκy⁻⁻(i, j, k, grid, loc, κ, clock, sl, b, C) = triad_mask_y(i, j,   j, k, k,   grid) * κᶜᶜᶜ(i, j, k, grid, loc, κ, clock) * tapering_factorᶜᶜᶜ(i, j, k, grid, sl, b, C)
+@inline ϵκy⁺⁺(i, j, k, grid, loc, κ, clock, sl, b, C) = triad_mask_y(i, j+1, j, k, k+1, grid) * κᶜᶜᶜ(i, j, k, grid, loc, κ, clock, C) * tapering_factorᶜᶜᶜ(i, j, k, grid, sl, b, C)
+@inline ϵκy⁺⁻(i, j, k, grid, loc, κ, clock, sl, b, C) = triad_mask_y(i, j+1, j, k, k,   grid) * κᶜᶜᶜ(i, j, k, grid, loc, κ, clock, C) * tapering_factorᶜᶜᶜ(i, j, k, grid, sl, b, C)
+@inline ϵκy⁻⁺(i, j, k, grid, loc, κ, clock, sl, b, C) = triad_mask_y(i, j,   j, k, k+1, grid) * κᶜᶜᶜ(i, j, k, grid, loc, κ, clock, C) * tapering_factorᶜᶜᶜ(i, j, k, grid, sl, b, C)
+@inline ϵκy⁻⁻(i, j, k, grid, loc, κ, clock, sl, b, C) = triad_mask_y(i, j,   j, k, k,   grid) * κᶜᶜᶜ(i, j, k, grid, loc, κ, clock, C) * tapering_factorᶜᶜᶜ(i, j, k, grid, sl, b, C)
 
 # Triad diagram key
 # =================
@@ -282,7 +286,7 @@ end
                  ϵκʸ⁻⁺ * Sy⁻⁺(i, j, k-1, grid, b, C) * ∂yᶜᶠᶜ(i, j,   k-1, grid, c) +
                  ϵκʸ⁺⁺ * Sy⁺⁺(i, j, k-1, grid, b, C) * ∂yᶜᶠᶜ(i, j+1, k-1, grid, c)) / 4
 
-    κϵ_R₃₃_∂z_c = explicit_R₃₃_∂z_c(i, j, k, grid, TD(), c, closure, b, C)
+    κϵ_R₃₃_∂z_c = explicit_R₃₃_∂z_c(i, j, k, grid, TD(), clock, c, closure, b, C)
 
     return - κR₃₁_∂x_c - κR₃₂_∂y_c - κϵ_R₃₃_∂z_c
 end
@@ -308,13 +312,13 @@ end
     return ϵκR₃₃
 end
 
-@inline function explicit_R₃₃_∂z_c(i, j, k, grid, ::ExplicitTimeDiscretization, c, closure, b, C) 
+@inline function explicit_R₃₃_∂z_c(i, j, k, grid, ::ExplicitTimeDiscretization, clock, c, closure, b, C) 
     κ  = closure.κ_symmetric
     sl = closure.slope_limiter
     return ϵκR₃₃(i, j, k, grid, κ, clock, sl, b, C) * ∂zᶜᶜᶠ(i, j, k, grid, c)
 end
 
-@inline explicit_R₃₃_∂z_c(i, j, k, grid, ::VerticallyImplicitTimeDiscretization, c, closure, b, C) = zero(grid)
+@inline explicit_R₃₃_∂z_c(i, j, k, grid, ::VerticallyImplicitTimeDiscretization, clock, c, closure, b, C) = zero(grid)
 
 @inline κzᶜᶜᶠ(i, j, k, grid, closure::FlavorOfTISSD, K, ::Val{id}, clock) where id = @inbounds K.ϵκR₃₃[i, j, k]
 

test/runtests.jl

@@ -138,6 +138,7 @@ CUDA.allowscalar() do
     if group == :turbulence_closures || group == :all
         @testset "Turbulence closures tests" begin
             include("test_turbulence_closures.jl")
+            include("test_triad_isopycnal_diffusivity.jl")
             include("test_gm_infinite_slope.jl")
         end
     end

test/test_triad_isopycnal_diffusivity.jl

@@ -0,0 +1,45 @@
+include("dependencies_for_runtests.jl")
+
+using Oceananigans.TurbulenceClosures: TriadIsopycnalSkewSymmetricDiffusivity
+using Oceananigans.TurbulenceClosures: diffusive_flux_x, diffusive_flux_y, diffusive_flux_z,
+                                       ExplicitTimeDiscretization, VerticallyImplicitTimeDiscretization,
+                                       compute_diffusivities!
+
+"""
+Test that TriadIsopycnalSkewSymmetricDiffusivity can be constructed and timestepped
+with both any time discretization.
+"""
+function time_step_with_triad_isopycnal_diffusivity(arch, time_discretization)
+    grid = RectilinearGrid(arch, size=(4, 4, 8), extent=(100, 100, 100))
+
+    closure = TriadIsopycnalSkewSymmetricDiffusivity(time_discretization, Float64,
+                                                     κ_skew = 100.0,
+                                                     κ_symmetric = 100.0)
+
+    # TriadIsopycnalSkewSymmetricDiffusivity only works with HydrostaticFreeSurfaceModel
+    model = HydrostaticFreeSurfaceModel(; grid, closure,
+                                        buoyancy = BuoyancyTracer(),
+                                        tracers = (:b, :c))
+
+    # A constant stratification initial condition
+    set!(model, b=(x, y, z) -> 1e-5 * z)
+
+    # Attempt to time-step
+    time_step!(model, 1)
+
+    return true
+end
+
+
+@testset "TriadIsopycnalSkewSymmetricDiffusivity" begin
+    @info "Testing TriadIsopycnalSkewSymmetricDiffusivity..."
+
+    for arch in archs
+        @testset "Time stepping with TriadIsopycnalSkewSymmetricDiffusivity [$arch]" begin
+            for time_discretization in [ExplicitTimeDiscretization(), VerticallyImplicitTimeDiscretization()]
+                @info "  Time-stepping TriadIsopycnalSkewSymmetricDiffusivity with $(typeof(time_discretization)) on $arch..."
+                @test time_step_with_triad_isopycnal_diffusivity(arch, time_discretization)
+           end
+        end
+    end
+end