Use triads (Griffies et al 1998) to compute the skew flux

src/TurbulenceClosures/TurbulenceClosures.jl

@@ -186,6 +186,7 @@ include("turbulence_closure_implementations/ri_based_vertical_diffusivity.jl")
 # Special non-abstracted diffusivities:
 # TODO: introduce abstract typing for these
 include("turbulence_closure_implementations/isopycnal_skew_symmetric_diffusivity.jl")
+include("turbulence_closure_implementations/isopycnal_skew_symmetric_diffusivity_with_triads.jl")
 include("turbulence_closure_implementations/leith_enstrophy_diffusivity.jl")
 
 using .TKEBasedVerticalDiffusivities: CATKEVerticalDiffusivity, TKEDissipationVerticalDiffusivity
@@ -201,3 +202,4 @@ include("turbulence_closure_diagnostics.jl")
 #####
 
 end # module
+

src/TurbulenceClosures/abstract_scalar_diffusivity_closure.jl

@@ -269,7 +269,6 @@ end
 #####
 
 # Number
-
 @inline νᶜᶜᶜ(i, j, k, grid, loc, ν::Number, args...) = ν
 @inline νᶠᶜᶠ(i, j, k, grid, loc, ν::Number, args...) = ν
 @inline νᶜᶠᶠ(i, j, k, grid, loc, ν::Number, args...) = ν
@@ -278,6 +277,7 @@ end
 @inline κᶠᶜᶜ(i, j, k, grid, loc, κ::Number, args...) = κ
 @inline κᶜᶠᶜ(i, j, k, grid, loc, κ::Number, args...) = κ
 @inline κᶜᶜᶠ(i, j, k, grid, loc, κ::Number, args...) = κ
+@inline κᶜᶜᶜ(i, j, k, grid, loc, κ::Number, args...) = κ
 
 # Array / Field at `Center, Center, Center`
 const Lᶜᶜᶜ = Tuple{Center, Center, Center}
@@ -289,6 +289,7 @@ const Lᶜᶜᶜ = Tuple{Center, Center, Center}
 @inline κᶠᶜᶜ(i, j, k, grid, ::Lᶜᶜᶜ, κ::AbstractArray, args...) = ℑxᶠᵃᵃ(i, j, k, grid, κ)
 @inline κᶜᶠᶜ(i, j, k, grid, ::Lᶜᶜᶜ, κ::AbstractArray, args...) = ℑyᵃᶠᵃ(i, j, k, grid, κ)
 @inline κᶜᶜᶠ(i, j, k, grid, ::Lᶜᶜᶜ, κ::AbstractArray, args...) = ℑzᵃᵃᶠ(i, j, k, grid, κ)
+@inline κᶜᶜᶜ(i, j, k, grid, ::Lᶜᶜᶜ, κ::AbstractArray, args...) = @inbounds κ[i, j, k]
 
 # Array / Field at `Center, Center, Face`
 const Lᶜᶜᶠ = Tuple{Center, Center, Face}
@@ -300,6 +301,7 @@ const Lᶜᶜᶠ = Tuple{Center, Center, Face}
 @inline κᶠᶜᶜ(i, j, k, grid, ::Lᶜᶜᶠ, κ::AbstractArray, args...) = ℑxzᶠᵃᶠ(i, j, k, grid, κ)
 @inline κᶜᶠᶜ(i, j, k, grid, ::Lᶜᶜᶠ, κ::AbstractArray, args...) = ℑyzᵃᶠᶠ(i, j, k, grid, κ)
 @inline κᶜᶜᶠ(i, j, k, grid, ::Lᶜᶜᶠ, κ::AbstractArray, args...) = @inbounds κ[i, j, k]
+@inline κᶜᶜᶜ(i, j, k, grid, ::Lᶜᶜᶠ, κ::AbstractArray, args...) = ℑzᵃᵃᶜ(i, j, k, grid, κ)
 
 # Function
 
@@ -314,6 +316,7 @@ const f = Face()
 @inline κᶠᶜᶜ(i, j, k, grid, loc, κ::F, clock, args...) where F<:Function = κ(node(i, j, k, grid, f, c, c)..., clock.time)
 @inline κᶜᶠᶜ(i, j, k, grid, loc, κ::F, clock, args...) where F<:Function = κ(node(i, j, k, grid, c, f, c)..., clock.time)
 @inline κᶜᶜᶠ(i, j, k, grid, loc, κ::F, clock, args...) where F<:Function = κ(node(i, j, k, grid, c, c, f)..., clock.time)
+@inline κᶜᶜᶜ(i, j, k, grid, loc, κ::F, clock, args...) where F<:Function = κ(node(i, j, k, grid, c, c, c)..., clock.time)
 
 # "DiscreteDiffusionFunction"
 @inline νᶜᶜᶜ(i, j, k, grid, loc, ν::DiscreteDiffusionFunction, clock, fields) = getdiffusivity(ν, i, j, k, grid, (c, c, c), clock, fields)
@@ -324,5 +327,6 @@ const f = Face()
 @inline κᶠᶜᶜ(i, j, k, grid, loc, κ::DiscreteDiffusionFunction, clock, fields) = getdiffusivity(κ, i, j, k, grid, (f, c, c), clock, fields)
 @inline κᶜᶠᶜ(i, j, k, grid, loc, κ::DiscreteDiffusionFunction, clock, fields) = getdiffusivity(κ, i, j, k, grid, (c, f, c), clock, fields)
 @inline κᶜᶜᶠ(i, j, k, grid, loc, κ::DiscreteDiffusionFunction, clock, fields) = getdiffusivity(κ, i, j, k, grid, (c, c, f), clock, fields)
+@inline κᶜᶜᶜ(i, j, k, grid, loc, κ::DiscreteDiffusionFunction, clock, fields) = getdiffusivity(κ, i, j, k, grid, (c, c, c), clock, fields)
 
 

src/TurbulenceClosures/turbulence_closure_implementations/isopycnal_skew_symmetric_diffusivity.jl

@@ -175,7 +175,6 @@ end
 end
 
 # Diffusive fluxes
-
 @inline get_tracer_κ(κ::NamedTuple, tracer_index) = @inbounds κ[tracer_index]
 @inline get_tracer_κ(κ, tracer_index) = κ
 

validation/mesoscale_turbulence/coarse_baroclinic_adjustment.jl → validation/isopycnal_skew_symmetric_diffusion/coarse_baroclinic_adjustment.jl

@@ -4,54 +4,47 @@ using Random
 using Oceananigans
 using Oceananigans.Units
 using GLMakie
+using Oceananigans.TurbulenceClosures: IsopycnalSkewSymmetricDiffusivity
+using Oceananigans.TurbulenceClosures: TriadIsopycnalSkewSymmetricDiffusivity
+using Oceananigans.TurbulenceClosures: FluxTapering
 
-gradient = "x"
+gradient = "y"
 filename = "coarse_baroclinic_adjustment_" * gradient
 
-# Architecture
-architecture = CPU()
-
 # Domain
 Ly = 1000kilometers  # north-south extent [m]
 Lz = 1kilometers     # depth [m]
 Ny = 20
 Nz = 20
 save_fields_interval = 0.5day
 stop_time = 30days
-Δt = 20minutes
+Δt = 10minutes
 
-grid = RectilinearGrid(architecture;
-                       topology = (Bounded, Bounded, Bounded),
-                       size = (Ny, Ny, Nz), 
+grid = RectilinearGrid(topology = (Periodic, Bounded, Bounded),
+                       size = (3, Ny, Nz), 
                        x = (-Ly/2, Ly/2),
                        y = (-Ly/2, Ly/2),
                        z = (-Lz, 0),
                        halo = (3, 3, 3))
 
 coriolis = FPlane(latitude = -45)
 
-Δy = Ly/Ny
-@show κh = νh = Δy^4 / 10days
-vertical_closure = VerticalScalarDiffusivity(ν=1e-2, κ=1e-4)
-horizontal_closure = HorizontalScalarBiharmonicDiffusivity(ν=νh, κ=κh)
-
 gerdes_koberle_willebrand_tapering = FluxTapering(1e-2)
-gent_mcwilliams_diffusivity = IsopycnalSkewSymmetricDiffusivity(κ_skew=1e3,
-                                                                κ_symmetric=1e3,
-                                                                slope_limiter=gerdes_koberle_willebrand_tapering)
+triad_closure = TriadIsopycnalSkewSymmetricDiffusivity(κ_skew = 1e3,
+                                                       κ_symmetric = 1e3,
+                                                       slope_limiter = gerdes_koberle_willebrand_tapering)
 
-closures = (vertical_closure, horizontal_closure, gent_mcwilliams_diffusivity)
+cox_closure = IsopycnalSkewSymmetricDiffusivity(VerticallyImplicitTimeDiscretization(),
+                                                κ_skew = 1e3,
+                                                κ_symmetric = 1e3,
+                                                slope_limiter = gerdes_koberle_willebrand_tapering)
 
 @info "Building a model..."
 
-model = HydrostaticFreeSurfaceModel(grid = grid,
-                                    coriolis = coriolis,
+model = HydrostaticFreeSurfaceModel(; grid, coriolis,
+                                    closure = triad_closure,
                                     buoyancy = BuoyancyTracer(),
-                                    closure = closures,
-                                    tracers = (:b, :c),
-                                    momentum_advection = WENO5(),
-                                    tracer_advection = WENO5(),
-                                    free_surface = ImplicitFreeSurface())
+                                    tracers = (:b, :c))
 
 @info "Built $model."
 
@@ -91,30 +84,25 @@ set!(model, b=bᵢ, c=cᵢ)
 
 simulation = Simulation(model; Δt, stop_time)
 
-# add timestep wizard callback
-wizard = TimeStepWizard(cfl=0.1, max_change=1.1, max_Δt=Δt)
-simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(20))
-
-# add progress callback
-wall_clock = [time_ns()]
+wall_clock = Ref(time_ns())
 
-function print_progress(sim)
+function progress(sim)
     @printf("[%05.2f%%] i: %d, t: %s, wall time: %s, max(u): (%6.3e, %6.3e, %6.3e) m/s, next Δt: %s\n",
             100 * (sim.model.clock.time / sim.stop_time),
             sim.model.clock.iteration,
             prettytime(sim.model.clock.time),
-            prettytime(1e-9 * (time_ns() - wall_clock[1])),
+            prettytime(1e-9 * (time_ns() - wall_clock[])),
             maximum(abs, sim.model.velocities.u),
             maximum(abs, sim.model.velocities.v),
             maximum(abs, sim.model.velocities.w),
             prettytime(sim.Δt))
 
-    wall_clock[1] = time_ns()
+    wall_clock[] = time_ns()
 
     return nothing
 end
 
-simulation.callbacks[:print_progress] = Callback(print_progress, IterationInterval(20))
+add_callback!(simulation, progress, IterationInterval(10))
 
 simulation.output_writers[:fields] = JLD2OutputWriter(model, merge(model.velocities, model.tracers),
                                                       schedule = TimeInterval(save_fields_interval),
@@ -140,7 +128,7 @@ bt = FieldTimeSeries(filepath, "b")
 ct = FieldTimeSeries(filepath, "c")
 
 # Build coordinates, rescaling the vertical coordinate
-x, y, z = nodes((Center, Center, Center), grid)
+x, y, z = nodes(bt)
 
 zscale = 1
 z = z .* zscale

validation/isopycnal_skew_symmetric_diffusion/new_front_relax.jl

@@ -0,0 +1,178 @@
+using Printf
+using Statistics
+using Random
+using Oceananigans
+using Oceananigans.Units
+using Oceananigans.Grids: znode
+using GLMakie
+using Oceananigans.TurbulenceClosures: IsopycnalSkewSymmetricDiffusivity
+using Oceananigans.TurbulenceClosures: TriadIsopycnalSkewSymmetricDiffusivity
+using Oceananigans.TurbulenceClosures: FluxTapering
+
+gradient = "y"
+filename = "new_front_relax_" * gradient
+
+# Domain
+dz=[50, 50, 55, 60, 65, 70, 80, 95, 120, 155, 200, 260, 320, 400, 480]
+Ly = 320kilometers  # north-south extent [m]
+Lz = sum(dz)        # depth [m]
+Nx = 3
+Ny = 32
+Nz = 15
+#save_fields_interval = 1day
+stop_time = 30days
+save_fields_interval = 1days
+Δt = 15minutes
+
+# viscosity
+viscAh=300.
+viscAz=2.e-4
+
+grid = RectilinearGrid(topology = (Flat, Bounded, Bounded),
+                       size = (Ny, Nz), 
+                       y = (-Ly/2, Ly/2),
+                       z = ([0 cumsum(reverse(dz))']'.-Lz),
+                       halo = (3, 3))
+
+#coriolis = FPlane(latitude = -45)
+coriolis = FPlane( f = 1.e-4)
+
+h_visc= HorizontalScalarDiffusivity( ν=viscAh )
+v_visc= VerticalScalarDiffusivity( ν=viscAz )
+
+gerdes_koberle_willebrand_tapering = FluxTapering(1e-2)
+
+triad_closure = TriadIsopycnalSkewSymmetricDiffusivity(VerticallyImplicitTimeDiscretization(), κ_skew = 0e3,
+                                                       κ_symmetric = 1e3,
+                                                       slope_limiter = gerdes_koberle_willebrand_tapering)
+
+cox_closure = IsopycnalSkewSymmetricDiffusivity(κ_skew = 0e3,
+                                                κ_symmetric = 1e3,
+                                                slope_limiter = gerdes_koberle_willebrand_tapering)
+
+@info "Building a model..."
+
+model = HydrostaticFreeSurfaceModel(; grid, coriolis,
+                                    closure = triad_closure, #, h_visc, v_visc),
+                                    buoyancy = BuoyancyTracer(),
+                                    tracers = (:b, :c))
+
+@info "Built $model."
+
+# Parameters
+N² = 4e-6 # [s⁻²] buoyancy frequency / stratification, corresponds to N/f = 20
+M² = 4e-9 # [s⁻²] horizontal buoyancy gradienta, gives a slope of 1.e-3
+
+# tracer patch centered on yC,zC:
+yC = 0
+zC = znode(6, grid, Center())
+Δy = Ly/8
+Δz = 500
+
+#Δc = 2Δy
+Δb = Ly/Ny * M²
+ϵb = 1e-2 * Δb # noise amplitude
+
+bᵢ(y, z) = N² * z - M² * Ly*sin(pi*y/Ly) *exp(-(3*z/Lz)^2)
+cᵢ(y, z) = exp( -((y-yC)/Δy)^2 )*exp( -((z-zC)/Δz).^2);
+
+set!(model, b=bᵢ, c=cᵢ)
+
+#####
+##### Simulation building
+#####
+
+simulation = Simulation(model; Δt, stop_time) #, stop_iteration = 10)
+
+wall_clock = Ref(time_ns())
+
+function progress(sim)
+    @printf("[%05.2f%%] i: %d, t: %s, wall time: %s, max(u): (%6.3e, %6.3e, %6.3e) m/s, next Δt: %s\n",
+            100 * (sim.model.clock.time / sim.stop_time),
+            sim.model.clock.iteration,
+            prettytime(sim.model.clock.time),
+            prettytime(1e-9 * (time_ns() - wall_clock[])),
+            maximum(abs, sim.model.velocities.u),
+            maximum(abs, sim.model.velocities.v),
+            maximum(abs, sim.model.velocities.w),
+            prettytime(sim.Δt))
+
+    wall_clock[] = time_ns()
+
+    return nothing
+end
+
+add_callback!(simulation, progress, IterationInterval(10))
+
+simulation.output_writers[:fields] = JLD2OutputWriter(model, merge(model.velocities, model.tracers),
+                                                      schedule = TimeInterval(save_fields_interval),
+                                                      filename = filename * "_fields",
+                                                      overwrite_existing = true)
+
+@info "Running the simulation..."
+
+run!(simulation)
+
+@info "Simulation completed in " * prettytime(simulation.run_wall_time)
+
+#####
+##### Visualize
+#####
+
+fig = Figure(size=(1400, 700))
+
+filepath = filename * "_fields.jld2"
+
+ut = FieldTimeSeries(filepath, "u")
+bt = FieldTimeSeries(filepath, "b")
+ct = FieldTimeSeries(filepath, "c")
+
+# Build coordinates, rescaling the vertical coordinate
+x, y, z = nodes(bt)
+
+zscale = 1
+z = z .* zscale
+
+#####
+##### Plot buoyancy...
+#####
+
+times = bt.times
+Nt = length(times)
+
+un(n) = interior(mean(ut[n], dims=1), 1, :, :)
+bn(n) = interior(mean(bt[n], dims=1), 1, :, :)
+cn(n) = interior(mean(ct[n], dims=1), 1, :, :)
+
+@show min_c = 0
+@show max_c = 1
+@show max_u = maximum(abs, un(Nt))
+min_u = - max_u
+
+axu = Axis(fig[2, 1], xlabel="$gradient (km)", ylabel="z (km)", title="Zonal velocity")
+axc = Axis(fig[3, 1], xlabel="$gradient (km)", ylabel="z (km)", title="Tracer concentration")
+slider = Slider(fig[4, 1:2], range=1:Nt, startvalue=1)
+n = slider.value
+
+u = @lift un($n)
+b = @lift bn($n)
+c = @lift cn($n)
+
+hm = heatmap!(axu, y * 1e-3, z * 1e-3, u, colorrange=(min_u, max_u), colormap=:balance)
+contour!(axu, y * 1e-3, z * 1e-3, b, levels = 25, color=:black, linewidth=2)
+cb = Colorbar(fig[2, 2], hm)
+
+hm = heatmap!(axc, y * 1e-3, z * 1e-3, c, colorrange=(0, 0.5), colormap=:speed)
+contour!(axc, y * 1e-3, z * 1e-3, b, levels = 25, color=:black, linewidth=2)
+cb = Colorbar(fig[3, 2], hm)
+
+title_str = @lift "Baroclinic adjustment with GM at t = " * prettytime(times[$n])
+ax_t = fig[1, 1:2] = Label(fig, title_str)
+
+display(fig)
+
+record(fig, filename * ".mp4", 1:Nt, framerate=8) do i
+    @info "Plotting frame $i of $Nt"
+    n[] = i
+end
+