diff --git a/src/TurbulenceClosures/TurbulenceClosures.jl b/src/TurbulenceClosures/TurbulenceClosures.jl index 2e988816a4..5253459759 100644 --- a/src/TurbulenceClosures/TurbulenceClosures.jl +++ b/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 + diff --git a/src/TurbulenceClosures/abstract_scalar_diffusivity_closure.jl b/src/TurbulenceClosures/abstract_scalar_diffusivity_closure.jl index 01e322645d..2c0930a9dc 100644 --- a/src/TurbulenceClosures/abstract_scalar_diffusivity_closure.jl +++ b/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) diff --git a/src/TurbulenceClosures/turbulence_closure_implementations/isopycnal_skew_symmetric_diffusivity.jl b/src/TurbulenceClosures/turbulence_closure_implementations/isopycnal_skew_symmetric_diffusivity.jl index 87f6ceff5f..f98a53301f 100644 --- a/src/TurbulenceClosures/turbulence_closure_implementations/isopycnal_skew_symmetric_diffusivity.jl +++ b/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) = κ diff --git a/src/TurbulenceClosures/turbulence_closure_implementations/isopycnal_skew_symmetric_diffusivity_with_triads.jl b/src/TurbulenceClosures/turbulence_closure_implementations/isopycnal_skew_symmetric_diffusivity_with_triads.jl new file mode 100644 index 0000000000..0aacb00910 --- /dev/null +++ b/src/TurbulenceClosures/turbulence_closure_implementations/isopycnal_skew_symmetric_diffusivity_with_triads.jl @@ -0,0 +1,397 @@ +struct TriadIsopycnalSkewSymmetricDiffusivity{TD, K, S, M, L, N} <: AbstractTurbulenceClosure{TD, N} + κ_skew :: K + κ_symmetric :: S + isopycnal_tensor :: M + slope_limiter :: L + + function TriadIsopycnalSkewSymmetricDiffusivity{TD, N}(κ_skew :: K, + κ_symmetric :: S, + isopycnal_tensor :: I, + slope_limiter :: L) where {TD, K, S, I, L, N} + + return new{TD, K, S, I, L, N}(κ_skew, κ_symmetric, isopycnal_tensor, slope_limiter) + end +end + +const TISSD{TD} = TriadIsopycnalSkewSymmetricDiffusivity{TD} where TD +const TISSDVector{TD} = AbstractVector{<:TISSD{TD}} where TD +const FlavorOfTISSD{TD} = Union{TISSD{TD}, TISSDVector{TD}} where TD + +""" + TriadIsopycnalSkewSymmetricDiffusivity([time_disc=VerticallyImplicitTimeDiscretization(), FT=Float64;] + κ_skew = 0, + κ_symmetric = 0, + isopycnal_tensor = SmallSlopeIsopycnalTensor(), + slope_limiter = FluxTapering(1e-2), + required_halo_size::Int = 1) + +Return parameters for an isopycnal skew-symmetric tracer diffusivity with skew diffusivity +`κ_skew` and symmetric diffusivity `κ_symmetric` that uses an `isopycnal_tensor` model for +for calculating the isopycnal slopes, and (optionally) applying a `slope_limiter` to the +calculated isopycnal slope values. + +Both `κ_skew` and `κ_symmetric` may be constants, arrays, fields, or functions of `(x, y, z, t)`. + +The formulation follows Griffies et al. (1998) + +References +========== +* Griffies, S. M., A. Gnanadesikan, R. C. Pacanowski, V. D. Larichev, J. K. Dukowicz, and R. D. Smith (1998) Isoneutral diffusion in a z-coordinate ocean model. _J. Phys. Oceanogr._, **28**, 805–830, doi:10.1175/1520-0485(1998)028<0805:IDIAZC>2.0.CO;2 +""" +function TriadIsopycnalSkewSymmetricDiffusivity(time_disc=ExplicitTimeDiscretization(), FT=Float64; + κ_skew = 0, + κ_symmetric = 0, + isopycnal_tensor = SmallSlopeIsopycnalTensor(), + slope_limiter = FluxTapering(1e-2), + required_halo_size::Int = 1) + + isopycnal_tensor isa SmallSlopeIsopycnalTensor || + error("Only isopycnal_tensor=SmallSlopeIsopycnalTensor() is currently supported.") + + TD = typeof(time_disc) + + return TriadIsopycnalSkewSymmetricDiffusivity{TD, required_halo_size}(convert_diffusivity(FT, κ_skew), + convert_diffusivity(FT, κ_symmetric), + isopycnal_tensor, + slope_limiter) +end + +TriadIsopycnalSkewSymmetricDiffusivity(FT::DataType; kw...) = + TriadIsopycnalSkewSymmetricDiffusivity(VerticallyImplicitTimeDiscretization(), FT; kw...) + +with_tracers(tracers, closure::TISSD{TD, N}) where {TD, N} = + TriadIsopycnalSkewSymmetricDiffusivity{TD, N}(closure.κ_skew, closure.κ_symmetric, closure.isopycnal_tensor, closure.slope_limiter) + +# For ensembles of closures +function with_tracers(tracers, closure_vector::TISSDVector) + arch = architecture(closure_vector) + + if arch isa Architectures.GPU + closure_vector = Vector(closure_vector) + end + + Ex = length(closure_vector) + closure_vector = [with_tracers(tracers, closure_vector[i]) for i=1:Ex] + + return on_architecture(arch, closure_vector) +end + +# Note: computing diffusivities at cell centers for now. +function DiffusivityFields(grid, tracer_names, bcs, ::FlavorOfTISSD{TD}) where TD + if TD() isa VerticallyImplicitTimeDiscretization + # Precompute the _tapered_ 33 component of the isopycnal rotation tensor + K = (; ϵκR₃₃ = ZFaceField(grid)) + else + return nothing + end + + return K +end + +function compute_diffusivities!(diffusivities, closure::FlavorOfTISSD{TD}, model; parameters = :xyz) where TD + + arch = model.architecture + grid = model.grid + clock = model.clock + tracers = model.tracers + buoyancy = model.buoyancy + + if TD() isa VerticallyImplicitTimeDiscretization + launch!(arch, grid, parameters, + triad_compute_tapered_R₃₃!, + diffusivities, grid, closure, clock, buoyancy, tracers) + end + + return nothing +end + +@kernel function triad_compute_tapered_R₃₃!(K, grid, closure, clock, b, C) + i, j, k, = @index(Global, NTuple) + closure = getclosure(i, j, closure) + κ = closure.κ_symmetric + @inbounds K.ϵκR₃₃[i, j, k] = ϵκR₃₃(i, j, k, grid, κ, clock, b, C) +end + +##### +##### _triads_ +##### +##### There are two horizontal slopes: Sx and Sy +##### +##### Both slopes are "located" at tracer cell centers. +##### +##### The slopes are computed by a directional derivative, which lends an +##### "orientation" to the slope. For example, the x-slope `Sx` computed +##### with a "+" directional derivative in x, and a "+" directional derivative +##### in z, is +##### +##### Sx⁺⁺ᵢₖ = Δz / Δx * (bᵢ₊₁ - bᵢ) / (bₖ₊₁ - bₖ) +##### +##### The superscript codes ⁺⁺, ⁺⁻, ⁻⁺, ⁻⁻, denote the direction of the derivative +##### in (h, z). +##### +##### from https://github.com/CliMA/Oceananigans.jl/blob/glw/homogeneous-bounded/src/TurbulenceClosures/turbulence_closure_implementations/isopycnal_potential_vorticity_diffusivity.jl +##### + +@inline function triad_Sx(ix, iz, j, kx, kz, grid, buoyancy, tracers) + bx = ∂x_b(ix, j, kx, grid, buoyancy, tracers) + bz = ∂z_b(iz, j, kz, grid, buoyancy, tracers) + bz = max(bz, zero(grid)) + return ifelse(bz == 0, zero(grid), - bx / bz) +end + +@inline function triad_Sy(i, jy, jz, ky, kz, grid, buoyancy, tracers) + by = ∂y_b(i, jy, ky, grid, buoyancy, tracers) + bz = ∂z_b(i, jz, kz, grid, buoyancy, tracers) + bz = max(bz, zero(grid)) + return ifelse(bz == 0, zero(grid), - by / bz) +end + +@inline Sx⁺⁺(i, j, k, grid, buoyancy, tracers) = triad_Sx(i+1, i, j, k, k+1, grid, buoyancy, tracers) +@inline Sx⁺⁻(i, j, k, grid, buoyancy, tracers) = triad_Sx(i+1, i, j, k, k, grid, buoyancy, tracers) +@inline Sx⁻⁺(i, j, k, grid, buoyancy, tracers) = triad_Sx(i, i, j, k, k+1, grid, buoyancy, tracers) +@inline Sx⁻⁻(i, j, k, grid, buoyancy, tracers) = triad_Sx(i, i, j, k, k, grid, buoyancy, tracers) + +@inline Sy⁺⁺(i, j, k, grid, buoyancy, tracers) = triad_Sy(i, j+1, j, k, k+1, grid, buoyancy, tracers) +@inline Sy⁺⁻(i, j, k, grid, buoyancy, tracers) = triad_Sy(i, j+1, j, k, k, grid, buoyancy, tracers) +@inline Sy⁻⁺(i, j, k, grid, buoyancy, tracers) = triad_Sy(i, j, j, k, k+1, grid, buoyancy, tracers) +@inline Sy⁻⁻(i, j, k, grid, buoyancy, tracers) = triad_Sy(i, j, j, k, k, grid, buoyancy, tracers) + +# We remove triads that live on a boundary (immersed or top / bottom / north / south / east / west) +@inline triad_mask_x(ix, iz, j, kx, kz, grid) = + peripheral_node(ix, j, kx, grid, Face(), Center(), Center()) | peripheral_node(iz, j, kz, grid, Center(), Center(), Face()) + +@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, 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, b, C) +@inline ϵκx⁺⁻(i, j, k, grid, loc, κ, clock, b, C) = triad_mask_x(i+1, i, j, k, k, grid) * κᶜᶜᶜ(i, j, k, grid, loc, κ, clock) * tapering_factorᶜᶜᶜ(i, j, k, grid, b, C) +@inline ϵκx⁻⁺(i, j, k, grid, loc, κ, clock, b, C) = triad_mask_x(i, i, j, k, k+1, grid) * κᶜᶜᶜ(i, j, k, grid, loc, κ, clock) * tapering_factorᶜᶜᶜ(i, j, k, grid, b, C) +@inline ϵκx⁻⁻(i, j, k, grid, loc, κ, clock, b, C) = triad_mask_x(i, i, j, k, k, grid) * κᶜᶜᶜ(i, j, k, grid, loc, κ, clock) * tapering_factorᶜᶜᶜ(i, j, k, grid, b, C) + +@inline ϵκy⁺⁺(i, j, k, grid, loc, κ, clock, 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, b, C) +@inline ϵκy⁺⁻(i, j, k, grid, loc, κ, clock, b, C) = triad_mask_y(i, j+1, j, k, k, grid) * κᶜᶜᶜ(i, j, k, grid, loc, κ, clock) * tapering_factorᶜᶜᶜ(i, j, k, grid, b, C) +@inline ϵκy⁻⁺(i, j, k, grid, loc, κ, clock, b, C) = triad_mask_y(i, j, j, k, k+1, grid) * κᶜᶜᶜ(i, j, k, grid, loc, κ, clock) * tapering_factorᶜᶜᶜ(i, j, k, grid, b, C) +@inline ϵκy⁻⁻(i, j, k, grid, loc, κ, clock, b, C) = triad_mask_y(i, j, j, k, k, grid) * κᶜᶜᶜ(i, j, k, grid, loc, κ, clock) * tapering_factorᶜᶜᶜ(i, j, k, grid, b, C) + +# Triad diagram key +# ================= +# +# * ┗ : Sx⁺⁺ / Sy⁺⁺ +# * ┛ : Sx⁻⁺ / Sy⁻⁺ +# * ┓ : Sx⁻⁻ / Sy⁻⁻ +# * ┏ : Sx⁺⁻ / Sy⁺⁻ +# + +# defined at fcc +@inline function diffusive_flux_x(i, j, k, grid, closure::FlavorOfTISSD, K, ::Val{id}, + c, clock, C, b) where id + + closure = getclosure(i, j, closure) + κ = closure.κ_symmetric + + loc = (Center(), Center(), Center()) + + ϵκ⁺⁺ = ϵκx⁺⁺(i-1, j, k, grid, loc, κ, clock, b, C) + ϵκ⁺⁻ = ϵκx⁺⁻(i-1, j, k, grid, loc, κ, clock, b, C) + ϵκ⁻⁺ = ϵκx⁻⁺(i, j, k, grid, loc, κ, clock, b, C) + ϵκ⁻⁻ = ϵκx⁻⁻(i, j, k, grid, loc, κ, clock, b, C) + + # Small slope approximation + ∂x_c = ∂xᶠᶜᶜ(i, j, k, grid, c) + + # i-1 i + # k+1 ------------- + # | | + # ┏┗ ∘ ┛┓ | k + # | | + # k ------|------| + + Fx = (ϵκ⁺⁺ * (∂x_c + Sx⁺⁺(i-1, j, k, grid, b, C) * ∂zᶜᶜᶠ(i-1, j, k+1, grid, c)) + + ϵκ⁺⁻ * (∂x_c + Sx⁺⁻(i-1, j, k, grid, b, C) * ∂zᶜᶜᶠ(i-1, j, k, grid, c)) + + ϵκ⁻⁺ * (∂x_c + Sx⁻⁺(i, j, k, grid, b, C) * ∂zᶜᶜᶠ(i, j, k+1, grid, c)) + + ϵκ⁻⁻ * (∂x_c + Sx⁻⁻(i, j, k, grid, b, C) * ∂zᶜᶜᶠ(i, j, k, grid, c))) / 4 + + return - Fx +end + +# defined at cfc +@inline function diffusive_flux_y(i, j, k, grid, closure::FlavorOfTISSD, K, ::Val{id}, + c, clock, C, b) where id + + closure = getclosure(i, j, closure) + κ = closure.κ_symmetric + + loc = (Center(), Center(), Center()) + + ∂y_c = ∂yᶜᶠᶜ(i, j, k, grid, c) + + ϵκ⁺⁺ = ϵκy⁺⁺(i, j-1, k, grid, loc, κ, clock, b, C) + ϵκ⁺⁻ = ϵκy⁺⁻(i, j-1, k, grid, loc, κ, clock, b, C) + ϵκ⁻⁺ = ϵκy⁻⁺(i, j, k, grid, loc, κ, clock, b, C) + ϵκ⁻⁻ = ϵκy⁻⁻(i, j, k, grid, loc, κ, clock, b, C) + + Fy = (ϵκ⁺⁺ * (∂y_c + Sy⁺⁺(i, j-1, k, grid, b, C) * ∂zᶜᶜᶠ(i, j-1, k+1, grid, c)) + + ϵκ⁺⁻ * (∂y_c + Sy⁺⁻(i, j-1, k, grid, b, C) * ∂zᶜᶜᶠ(i, j-1, k, grid, c)) + + ϵκ⁻⁺ * (∂y_c + Sy⁻⁺(i, j, k, grid, b, C) * ∂zᶜᶜᶠ(i, j, k+1, grid, c)) + + ϵκ⁻⁻ * (∂y_c + Sy⁻⁻(i, j, k, grid, b, C) * ∂zᶜᶜᶠ(i, j, k, grid, c))) / 4 + + return - Fy +end + +# defined at ccf +@inline function diffusive_flux_z(i, j, k, grid, closure::FlavorOfTISSD{TD}, K, ::Val{id}, + c, clock, C, b) where {TD, id} + + closure = getclosure(i, j, closure) + κ = closure.κ_symmetric + + loc = (Center(), Center(), Center()) + + ϵκˣ⁻⁻ = ϵκx⁻⁻(i, j, k, grid, loc, κ, clock, b, C) + ϵκˣ⁺⁻ = ϵκx⁺⁻(i, j, k, grid, loc, κ, clock, b, C) + ϵκˣ⁻⁺ = ϵκx⁻⁺(i, j, k-1, grid, loc, κ, clock, b, C) + ϵκˣ⁺⁺ = ϵκx⁺⁺(i, j, k-1, grid, loc, κ, clock, b, C) + + ϵκʸ⁻⁻ = ϵκy⁻⁻(i, j, k, grid, loc, κ, clock, b, C) + ϵκʸ⁺⁻ = ϵκy⁺⁻(i, j, k, grid, loc, κ, clock, b, C) + ϵκʸ⁻⁺ = ϵκy⁻⁺(i, j, k-1, grid, loc, κ, clock, b, C) + ϵκʸ⁺⁺ = ϵκy⁺⁺(i, j, k-1, grid, loc, κ, clock, b, C) + + # Triad diagram: + # + # i-1 i i+1 + # ------------------- + # | | | | + # | | ┓ ┏ | k | + # | | | | + # - k -- ∘ -- - + # | | | | + # | | ┛ ┗ | k-1 | + # | | | | + # -------------------- + + κR₃₁_∂x_c = (ϵκˣ⁻⁻ * Sx⁻⁻(i, j, k, grid, b, C) * ∂xᶠᶜᶜ(i, j, k, grid, c) + + ϵκˣ⁺⁻ * Sx⁺⁻(i, j, k, grid, b, C) * ∂xᶠᶜᶜ(i+1, j, k, grid, c) + + ϵκˣ⁻⁺ * Sx⁻⁺(i, j, k-1, grid, b, C) * ∂xᶠᶜᶜ(i, j, k-1, grid, c) + + ϵκˣ⁺⁺ * Sx⁺⁺(i, j, k-1, grid, b, C) * ∂xᶠᶜᶜ(i+1, j, k-1, grid, c)) / 4 + + κR₃₂_∂y_c = (ϵκʸ⁻⁻ * Sy⁻⁻(i, j, k, grid, b, C) * ∂yᶜᶠᶜ(i, j, k, grid, c) + + ϵκʸ⁺⁻ * Sy⁺⁻(i, j, k, grid, b, C) * ∂yᶜᶠᶜ(i, j+1, k, grid, c) + + ϵκʸ⁻⁺ * 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) + + return - κR₃₁_∂x_c - κR₃₂_∂y_c - κϵ_R₃₃_∂z_c +end + +@inline function ϵκR₃₃(i, j, k, grid, κ, clock, b, C) + loc = (Center(), Center(), Center()) + + ϵκˣ⁻⁻ = ϵκx⁻⁻(i, j, k, grid, loc, κ, clock, b, C) + ϵκˣ⁺⁻ = ϵκx⁺⁻(i, j, k, grid, loc, κ, clock, b, C) + ϵκˣ⁻⁺ = ϵκx⁻⁺(i, j, k-1, grid, loc, κ, clock, b, C) + ϵκˣ⁺⁺ = ϵκx⁺⁺(i, j, k-1, grid, loc, κ, clock, b, C) + + ϵκʸ⁻⁻ = ϵκy⁻⁻(i, j, k, grid, loc, κ, clock, b, C) + ϵκʸ⁺⁻ = ϵκy⁺⁻(i, j, k, grid, loc, κ, clock, b, C) + ϵκʸ⁻⁺ = ϵκy⁻⁺(i, j, k-1, grid, loc, κ, clock, b, C) + ϵκʸ⁺⁺ = ϵκy⁺⁺(i, j, k-1, grid, loc, κ, clock, b, C) + + ϵκR₃₃ = (ϵκˣ⁻⁻ * Sx⁻⁻(i, j, k, grid, b, C)^2 + ϵκʸ⁻⁻ * Sy⁻⁻(i, j, k, grid, b, C)^2 + + ϵκˣ⁺⁻ * Sx⁺⁻(i, j, k, grid, b, C)^2 + ϵκʸ⁺⁻ * Sy⁺⁻(i, j, k, grid, b, C)^2 + + ϵκˣ⁻⁺ * Sx⁻⁺(i, j, k-1, grid, b, C)^2 + ϵκʸ⁻⁺ * Sy⁻⁺(i, j, k-1, grid, b, C)^2 + + ϵκˣ⁺⁺ * Sx⁺⁺(i, j, k-1, grid, b, C)^2 + ϵκʸ⁺⁺ * Sy⁺⁺(i, j, k-1, grid, b, C)^2) / 4 + + return ϵκR₃₃ +end + +@inline function explicit_R₃₃_∂z_c(i, j, k, grid, ::ExplicitTimeDiscretization, c, closure, b, C) + κ = closure.κ_symmetric + return ϵκR₃₃(i, j, k, grid, κ, clock, 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 κzᶜᶜᶠ(i, j, k, grid, closure::FlavorOfTISSD, K, ::Val{id}, clock) where id = @inbounds K.ϵκR₃₃[i, j, k] + +@inline viscous_flux_ux(i, j, k, grid, closure::Union{TISSD, TISSDVector}, args...) = zero(grid) +@inline viscous_flux_uy(i, j, k, grid, closure::Union{TISSD, TISSDVector}, args...) = zero(grid) +@inline viscous_flux_uz(i, j, k, grid, closure::Union{TISSD, TISSDVector}, args...) = zero(grid) + +@inline viscous_flux_vx(i, j, k, grid, closure::Union{TISSD, TISSDVector}, args...) = zero(grid) +@inline viscous_flux_vy(i, j, k, grid, closure::Union{TISSD, TISSDVector}, args...) = zero(grid) +@inline viscous_flux_vz(i, j, k, grid, closure::Union{TISSD, TISSDVector}, args...) = zero(grid) + +@inline viscous_flux_wx(i, j, k, grid, closure::Union{TISSD, TISSDVector}, args...) = zero(grid) +@inline viscous_flux_wy(i, j, k, grid, closure::Union{TISSD, TISSDVector}, args...) = zero(grid) +@inline viscous_flux_wz(i, j, k, grid, closure::Union{TISSD, TISSDVector}, args...) = zero(grid) + +##### +##### Show +##### + +Base.summary(closure::TISSD) = string("TriadIsopycnalSkewSymmetricDiffusivity", + "(κ_skew=", + prettysummary(closure.κ_skew), + ", κ_symmetric=", prettysummary(closure.κ_symmetric), ")") + +Base.show(io::IO, closure::TISSD) = + print(io, "TriadIsopycnalSkewSymmetricDiffusivity: " * + "(κ_symmetric=$(closure.κ_symmetric), κ_skew=$(closure.κ_skew), " * + "(isopycnal_tensor=$(closure.isopycnal_tensor), slope_limiter=$(closure.slope_limiter))") + +@inline tapering_factorᶜᶜᶜ(i, j, k, grid, b, C) = one(grid) + +# """ +# taper_factor(i, j, k, grid, closure, tracers, buoyancy) + +# Return the tapering factor `min(1, Sₘₐₓ² / slope²)`, where `slope² = slope_x² + slope_y²` +# that multiplies all components of the isopycnal slope tensor. The tapering factor is calculated on all the +# faces involved in the isopycnal slope tensor calculation. The minimum value of tapering is selected. + +# References +# ========== +# R. Gerdes, C. Koberle, and J. Willebrand. (1991), "The influence of numerical advection schemes +# on the results of ocean general circulation models", Clim. Dynamics, 5 (4), 211–226. +# """ +# @inline function tapering_factor(i, j, k, grid, closure, tracers, buoyancy) +# ϵᶠᶜᶜ = tapering_factorᶠᶜᶜ(i, j, k, grid, closure, tracers, buoyancy) +# ϵᶜᶠᶜ = tapering_factorᶜᶠᶜ(i, j, k, grid, closure, tracers, buoyancy) +# ϵᶜᶜᶠ = tapering_factorᶜᶜᶠ(i, j, k, grid, closure, tracers, buoyancy) +# return min(ϵᶠᶜᶜ, ϵᶜᶠᶜ, ϵᶜᶜᶠ) +# end + +# @inline function tapering_factorᶠᶜᶜ(i, j, k, grid, closure, tracers, buoyancy) +# by = ℑxyᶠᶜᵃ(i, j, k, grid, ∂y_b, buoyancy, tracers) +# bz = ℑxzᶠᵃᶜ(i, j, k, grid, ∂z_b, buoyancy, tracers) +# bx = ∂x_b(i, j, k, grid, buoyancy, tracers) +# return calc_tapering(bx, by, bz, grid, closure.isopycnal_tensor, closure.slope_limiter) +# end + +# @inline function tapering_factorᶜᶠᶜ(i, j, k, grid, closure, tracers, buoyancy) +# bx = ℑxyᶜᶠᵃ(i, j, k, grid, ∂x_b, buoyancy, tracers) +# bz = ℑyzᵃᶠᶜ(i, j, k, grid, ∂z_b, buoyancy, tracers) +# by = ∂y_b(i, j, k, grid, buoyancy, tracers) +# return calc_tapering(bx, by, bz, grid, closure.isopycnal_tensor, closure.slope_limiter) +# end + +# @inline function tapering_factorᶜᶜᶠ(i, j, k, grid, closure, tracers, buoyancy) +# bx = ℑxzᶜᵃᶠ(i, j, k, grid, ∂x_b, buoyancy, tracers) +# by = ℑyzᵃᶜᶠ(i, j, k, grid, ∂y_b, buoyancy, tracers) +# bz = ∂z_b(i, j, k, grid, buoyancy, tracers) +# return calc_tapering(bx, by, bz, grid, closure.isopycnal_tensor, closure.slope_limiter) +# end + +# @inline function calc_tapering(bx, by, bz, grid, slope_model, slope_limiter) + +# bz = max(bz, slope_model.minimum_bz) + +# slope_x = - bx / bz +# slope_y = - by / bz + +# # in case of a stable buoyancy gradient (bz > 0), the slope is set to zero +# slope² = ifelse(bz <= 0, zero(grid), slope_x^2 + slope_y^2) + +# return min(one(grid), slope_limiter.max_slope^2 / slope²) +# end diff --git a/validation/mesoscale_turbulence/coarse_baroclinic_adjustment.jl b/validation/isopycnal_skew_symmetric_diffusion/coarse_baroclinic_adjustment.jl similarity index 71% rename from validation/mesoscale_turbulence/coarse_baroclinic_adjustment.jl rename to validation/isopycnal_skew_symmetric_diffusion/coarse_baroclinic_adjustment.jl index cd26a180dc..e72aa85464 100644 --- a/validation/mesoscale_turbulence/coarse_baroclinic_adjustment.jl +++ b/validation/isopycnal_skew_symmetric_diffusion/coarse_baroclinic_adjustment.jl @@ -4,13 +4,13 @@ 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] @@ -18,11 +18,10 @@ 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), @@ -30,28 +29,22 @@ grid = RectilinearGrid(architecture; 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 diff --git a/validation/mesoscale_turbulence/coarse_lat_lon_baroclinic_adjustment.jl b/validation/isopycnal_skew_symmetric_diffusion/coarse_lat_lon_baroclinic_adjustment.jl similarity index 100% rename from validation/mesoscale_turbulence/coarse_lat_lon_baroclinic_adjustment.jl rename to validation/isopycnal_skew_symmetric_diffusion/coarse_lat_lon_baroclinic_adjustment.jl diff --git a/validation/isopycnal_skew_symmetric_diffusion/new_front_relax.jl b/validation/isopycnal_skew_symmetric_diffusion/new_front_relax.jl new file mode 100644 index 0000000000..5fbae78de7 --- /dev/null +++ b/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 + diff --git a/validation/mesoscale_turbulence/zonally_averaged_baroclinic_adjustment.jl b/validation/isopycnal_skew_symmetric_diffusion/zonally_averaged_baroclinic_adjustment.jl similarity index 100% rename from validation/mesoscale_turbulence/zonally_averaged_baroclinic_adjustment.jl rename to validation/isopycnal_skew_symmetric_diffusion/zonally_averaged_baroclinic_adjustment.jl