From f07c3c10f65d3e8aafe941bc10f12c3c45496604 Mon Sep 17 00:00:00 2001 From: Gregory Wagner Date: Tue, 5 Nov 2024 14:56:07 -0700 Subject: [PATCH 01/15] Get triads working --- mwe.jl | 13 + src/TurbulenceClosures/TurbulenceClosures.jl | 2 + ..._skew_symmetric_diffusivity_with_triads.jl | 473 ++++++++++++++++++ .../coarse_baroclinic_adjustment.jl | 56 +-- .../coarse_lat_lon_baroclinic_adjustment.jl | 0 .../zonally_averaged_baroclinic_adjustment.jl | 0 6 files changed, 509 insertions(+), 35 deletions(-) create mode 100644 mwe.jl create mode 100644 src/TurbulenceClosures/turbulence_closure_implementations/isopycnal_skew_symmetric_diffusivity_with_triads.jl rename validation/{mesoscale_turbulence => isopycnal_skew_symmetric_diffusion}/coarse_baroclinic_adjustment.jl (71%) rename validation/{mesoscale_turbulence => isopycnal_skew_symmetric_diffusion}/coarse_lat_lon_baroclinic_adjustment.jl (100%) rename validation/{mesoscale_turbulence => isopycnal_skew_symmetric_diffusion}/zonally_averaged_baroclinic_adjustment.jl (100%) diff --git a/mwe.jl b/mwe.jl new file mode 100644 index 0000000000..212ffe6a62 --- /dev/null +++ b/mwe.jl @@ -0,0 +1,13 @@ +using Oceananigans +using Oceananigans.Solvers: ConjugateGradientPoissonSolver, fft_poisson_solver +using Oceananigans.Models.NonhydrostaticModels: calculate_pressure_correction! + +N = 2 +x = y = (0, 1) +z = [0, 0.2, 1] +grid = RectilinearGrid(size=(N, N, N); x, y, z, halo=(2, 2, 2), topology=(Bounded, Periodic, Bounded)) +fft_solver = fft_poisson_solver(grid) +pressure_solver = ConjugateGradientPoissonSolver(grid, preconditioner=fft_poisson_solver(grid)) +model = NonhydrostaticModel(; grid, pressure_solver) +set!(model, u=1) + diff --git a/src/TurbulenceClosures/TurbulenceClosures.jl b/src/TurbulenceClosures/TurbulenceClosures.jl index 5b0aabcac3..44bd23d1fe 100644 --- a/src/TurbulenceClosures/TurbulenceClosures.jl +++ b/src/TurbulenceClosures/TurbulenceClosures.jl @@ -175,6 +175,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 @@ -188,3 +189,4 @@ include("turbulence_closure_diagnostics.jl") ##### end # module + 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..a008241e1f --- /dev/null +++ b/src/TurbulenceClosures/turbulence_closure_implementations/isopycnal_skew_symmetric_diffusivity_with_triads.jl @@ -0,0 +1,473 @@ +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 +const issd_coefficient_loc = (Center(), Center(), Center()) + +""" + TriadIsopycnalSkewSymmetricDiffusivity([time_disc=VerticallyImplicitTimeDiscretization(), FT=Float64;] + κ_skew = 0, + κ_symmetric = 0, + isopycnal_tensor = SmallSlopeIsopycnalTensor(), + slope_limiter = FluxTapering(1e-2)) + +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)`. +""" +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...) + +function with_tracers(tracers, closure::TISSD{TD, N}) where {TD, N} + κ_skew = !isa(closure.κ_skew, NamedTuple) ? closure.κ_skew : tracer_diffusivities(tracers, closure.κ_skew) + κ_symmetric = !isa(closure.κ_symmetric, NamedTuple) ? closure.κ_symmetric : tracer_diffusivities(tracers, closure.κ_symmetric) + return TriadIsopycnalSkewSymmetricDiffusivity{TD, N}(κ_skew, κ_symmetric, closure.isopycnal_tensor, closure.slope_limiter) +end + +# 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, closure::FlavorOfTISSD{TD}) where TD + if TD() isa VerticallyImplicitTimeDiscretization + # Precompute the _tapered_ 33 component of the isopycnal rotation tensor + return (; ϵ_R₃₃ = Field((Center, Center, Face), grid)) + else + return nothing + end +end + +function compute_diffusivities!(diffusivities, closure::FlavorOfTISSD{TD}, model; parameters = :xyz) where TD + + arch = model.architecture + grid = model.grid + tracers = model.tracers + buoyancy = model.buoyancy + + if TD() isa VerticallyImplicitTimeDiscretization + launch!(arch, grid, parameters, + triad_compute_tapered_R₃₃!, + diffusivities.ϵ_R₃₃, grid, closure, tracers, buoyancy) + end + + return nothing +end + +@kernel function triad_compute_tapered_R₃₃!(ϵ_R₃₃, grid, closure, tracers, buoyancy) + i, j, k, = @index(Global, NTuple) + closure = getclosure(i, j, closure) + R₃₃ = isopycnal_rotation_tensor_zz_ccf(i, j, k, grid, buoyancy, tracers, closure.isopycnal_tensor) + ϵ = tapering_factorᶜᶜᶠ(i, j, k, grid, closure, tracers, buoyancy) + @inbounds ϵ_R₃₃[i, j, k] = ϵ * R₃₃ +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, jx, kx, iz, jz, kz, grid, buoyancy, tracers) + bx = ∂x_b(ix, jx, kx, grid, buoyancy, tracers) + bz = ∂z_b(iz, jz, kz, grid, buoyancy, tracers) + bz = max(bz, zero(grid)) + return ifelse(bz == 0, zero(grid), - bx / bz) +end + +@inline function triad_Sy(iy, jy, ky, iz, jz, kz, grid, buoyancy, tracers) + by = ∂y_b(iy, jy, ky, grid, buoyancy, tracers) + bz = ∂z_b(iz, 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, j, k, i, j, k+1, grid, buoyancy, tracers) +@inline Sx⁺⁻(i, j, k, grid, buoyancy, tracers) = triad_Sx(i+1, j, k, i, j, k, grid, buoyancy, tracers) +@inline Sx⁻⁺(i, j, k, grid, buoyancy, tracers) = triad_Sx(i, j, k, i, j, k+1, grid, buoyancy, tracers) +@inline Sx⁻⁻(i, j, k, grid, buoyancy, tracers) = triad_Sx(i, j, k, i, j, k, grid, buoyancy, tracers) + +@inline Sy⁺⁺(i, j, k, grid, buoyancy, tracers) = triad_Sy(i, j+1, k, i, j, k+1, grid, buoyancy, tracers) +@inline Sy⁺⁻(i, j, k, grid, buoyancy, tracers) = triad_Sy(i, j+1, k, i, j, k, grid, buoyancy, tracers) +@inline Sy⁻⁺(i, j, k, grid, buoyancy, tracers) = triad_Sy(i, j, k, i, j, k+1, grid, buoyancy, tracers) +@inline Sy⁻⁻(i, j, k, grid, buoyancy, tracers) = triad_Sy(i, j, k, i, j, k, grid, buoyancy, tracers) + +# Tensor components +@inline Sxᶠᶜᶜ(i, j, k, grid, args...) = (Sx⁻⁺(i, j, k, grid, args...) + Sx⁻⁻(i, j, k, grid, args...) + # west cell, z-average + Sx⁻⁺(i-1, j, k, grid, args...) + Sx⁻⁻(i-1, j, k, grid, args...)) / 4 # east cell, z-average + +@inline Syᶜᶠᶜ(i, j, k, grid, args...) = (Sy⁻⁺(i, j, k, grid, args...) + Sy⁻⁻(i, j, k, grid, args...) + # south cell, z-average + Sy⁻⁺(i, j-1, k, grid, args...) + Sy⁻⁻(i, j-1, k, grid, args...)) / 4 # north cell, z-average + +@inline Syᶜᶜᶠ(i, j, k, grid, args...) = (Sy⁺⁻(i, j, k, grid, args...) + Sy⁻⁻(i, j, k, grid, args...) + # top cell, y-average + Sy⁺⁺(i, j, k-1, grid, args...) + Sy⁻⁺(i, j, k-1, grid, args...)) / 4 # bottom cell, y-average + +@inline Sxᶜᶜᶠ(i, j, k, grid, args...) = (Sx⁺⁻(i, j, k, grid, args...) + Sx⁻⁻(i, j, k, grid, args...) + # top cell, x-average + Sx⁺⁺(i, j, k-1, grid, args...) + Sx⁻⁺(i, j, k-1, grid, args...)) / 4 # bottom cell, x-average + +@inline Syᶜᶜᶠ(i, j, k, grid, args...) = (Sy⁺⁻(i, j, k, grid, args...) + Sy⁻⁻(i, j, k, grid, args...) + # top cell, y-average + Sy⁺⁺(i, j, k-1, grid, args...) + Sy⁻⁺(i, j, k-1, grid, args...)) / 4 # bottom cell, y-average + + + + +struct FluxTapering{FT} + max_slope :: FT +end + +""" + 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 + +# Diffusive fluxes + +@inline get_tracer_κ(κ::NamedTuple, tracer_index) = @inbounds κ[tracer_index] +@inline get_tracer_κ(κ, tracer_index) = κ + + +# 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) + + κ = get_tracer_κ(closure.κ_skew, id) + κ = κᶠᶜᶜ(i, j, k, grid, issd_coefficient_loc, κ, clock) + + # Small slope approximation + R₁₁_∂x_c = ∂xᶠᶜᶜ(i, j, k, grid, c) + R₁₂_∂y_c = zero(grid) + + # i-1 i + # k+1 ------------- + # | | + # ┏┗ ∘ ┛┓ | k + # | | + # k ------|------| + + R₁₃_∂z_c = (Sx⁺⁺(i-1, j, k, grid, b, C) * ∂zᶜᶜᶠ(i-1, j, k+1, grid, c) + + Sx⁺⁻(i-1, j, k, grid, b, C) * ∂zᶜᶜᶠ(i-1, j, k, grid, c) + + Sx⁻⁺(i, j, k, grid, b, C) * ∂zᶜᶜᶠ(i, j, k+1, grid, c) + + Sx⁻⁻(i, j, k, grid, b, C) * ∂zᶜᶜᶠ(i, j, k, grid, c)) / 4 + + return - κ * (R₁₁_∂x_c + R₁₂_∂y_c + R₁₃_∂z_c) +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) + + κ = get_tracer_κ(closure.κ_skew, id) + κ = κᶜᶠᶜ(i, j, k, grid, issd_coefficient_loc, κ, clock) + + # Small slope approximation + R₂₁_∂x_c = zero(grid) + R₂₂_∂y_c = ∂yᶜᶠᶜ(i, j, k, grid, c) + + R₂₃_∂z_c = (Sy⁺⁺(i, j-1, k, grid, b, C) * ∂zᶜᶜᶠ(i, j-1, k+1, grid, c) + + Sy⁺⁻(i, j-1, k, grid, b, C) * ∂zᶜᶜᶠ(i, j-1, k, grid, c) + + Sy⁻⁺(i, j, k, grid, b, C) * ∂zᶜᶜᶠ(i, j, k+1, grid, c) + + Sy⁻⁻(i, j, k, grid, b, C) * ∂zᶜᶜᶠ(i, j, k, grid, c)) / 4 + + return - κ * (R₂₁_∂x_c + R₂₂_∂y_c + R₂₃_∂z_c) +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) + + κ = get_tracer_κ(closure.κ_skew, id) + κ = κᶜᶜᶠ(i, j, k, grid, issd_coefficient_loc, κ, clock) + + # 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 explicit_R₃₃_∂z_c(i, j, k, grid, ::ExplicitTimeDiscretization, c, closure, b, C) + ϵ = tapering_factorᶜᶜᶠ(i, j, k, grid, closure, C, b) + ∂z_c = ∂zᶜᶜᶠ(i, j, k, grid, c) + Sx = Sxᶜᶜᶠ(i, j, k, grid, b, C) + Sy = Syᶜᶜᶠ(i, j, k, grid, b, C) + R₃₃ = Sx^2 + Sy^2 + return ϵ * R₃₃ * ∂z_c +end + +#= +# defined at fcc +@inline function diffusive_flux_x(i, j, k, grid, + closure::Union{TISSD, TISSDVector}, diffusivity_fields, ::Val{tracer_index}, + c, clock, fields, buoyancy) where tracer_index + + closure = getclosure(i, j, closure) + + κ_skew = get_tracer_κ(closure.κ_skew, tracer_index) + κ_symmetric = get_tracer_κ(closure.κ_symmetric, tracer_index) + + κ_skewᶠᶜᶜ = κᶠᶜᶜ(i, j, k, grid, issd_coefficient_loc, κ_skew, clock) + κ_symmetricᶠᶜᶜ = κᶠᶜᶜ(i, j, k, grid, issd_coefficient_loc, κ_symmetric, clock) + + ∂x_c = ∂xᶠᶜᶜ(i, j, k, grid, c) + + # Average... of... the gradient! + ∂y_c = ℑxyᶠᶜᵃ(i, j, k, grid, ∂yᶜᶠᶜ, c) + ∂z_c = ℑxzᶠᵃᶜ(i, j, k, grid, ∂zᶜᶜᶠ, c) + + R₁₁ = one(grid) + R₁₂ = zero(grid) + R₁₃ = isopycnal_rotation_tensor_xz_fcc(i, j, k, grid, buoyancy, fields, closure.isopycnal_tensor) + + ϵ = tapering_factorᶠᶜᶜ(i, j, k, grid, closure, fields, buoyancy) + + return - ϵ * ( κ_symmetricᶠᶜᶜ * R₁₁ * ∂x_c + + κ_symmetricᶠᶜᶜ * R₁₂ * ∂y_c + + (κ_symmetricᶠᶜᶜ - κ_skewᶠᶜᶜ) * R₁₃ * ∂z_c) +end + +# defined at cfc +@inline function diffusive_flux_y(i, j, k, grid, + closure::Union{TISSD, TISSDVector}, diffusivity_fields, ::Val{tracer_index}, + c, clock, fields, buoyancy) where tracer_index + + closure = getclosure(i, j, closure) + + κ_skew = get_tracer_κ(closure.κ_skew, tracer_index) + κ_symmetric = get_tracer_κ(closure.κ_symmetric, tracer_index) + + κ_skewᶜᶠᶜ = κᶜᶠᶜ(i, j, k, grid, issd_coefficient_loc, κ_skew, clock) + κ_symmetricᶜᶠᶜ = κᶜᶠᶜ(i, j, k, grid, issd_coefficient_loc, κ_symmetric, clock) + + ∂y_c = ∂yᶜᶠᶜ(i, j, k, grid, c) + + # Average... of... the gradient! + ∂x_c = ℑxyᶜᶠᵃ(i, j, k, grid, ∂xᶠᶜᶜ, c) + ∂z_c = ℑyzᵃᶠᶜ(i, j, k, grid, ∂zᶜᶜᶠ, c) + + R₂₁ = zero(grid) + R₂₂ = one(grid) + R₂₃ = isopycnal_rotation_tensor_yz_cfc(i, j, k, grid, buoyancy, fields, closure.isopycnal_tensor) + + ϵ = tapering_factorᶜᶠᶜ(i, j, k, grid, closure, fields, buoyancy) + + return - ϵ * (κ_symmetricᶜᶠᶜ * R₂₁ * ∂x_c + + κ_symmetricᶜᶠᶜ * R₂₂ * ∂y_c + + (κ_symmetricᶜᶠᶜ - κ_skewᶜᶠᶜ) * R₂₃ * ∂z_c) +end + +# defined at ccf +@inline function diffusive_flux_z(i, j, k, grid, + closure::FlavorOfTISSD{TD}, diffusivity_fields, ::Val{tracer_index}, + c, clock, fields, buoyancy) where {tracer_index, TD} + + closure = getclosure(i, j, closure) + + κ_skew = get_tracer_κ(closure.κ_skew, tracer_index) + κ_symmetric = get_tracer_κ(closure.κ_symmetric, tracer_index) + + κ_skewᶜᶜᶠ = κᶜᶜᶠ(i, j, k, grid, issd_coefficient_loc, κ_skew, clock) + κ_symmetricᶜᶜᶠ = κᶜᶜᶠ(i, j, k, grid, issd_coefficient_loc, κ_symmetric, clock) + + # Average... of... the gradient! + ∂x_c = ℑxzᶜᵃᶠ(i, j, k, grid, ∂xᶠᶜᶜ, c) + ∂y_c = ℑyzᵃᶜᶠ(i, j, k, grid, ∂yᶜᶠᶜ, c) + + R₃₁ = isopycnal_rotation_tensor_xz_ccf(i, j, k, grid, buoyancy, fields, closure.isopycnal_tensor) + R₃₂ = isopycnal_rotation_tensor_yz_ccf(i, j, k, grid, buoyancy, fields, closure.isopycnal_tensor) + + κ_symmetric_∂z_c = explicit_κ_∂z_c(i, j, k, grid, TD(), c, κ_symmetricᶜᶜᶠ, closure, buoyancy, fields) + + ϵ = tapering_factorᶜᶜᶠ(i, j, k, grid, closure, fields, buoyancy) + + return - ϵ * κ_symmetric_∂z_c - ϵ * ((κ_symmetricᶜᶜᶠ + κ_skewᶜᶜᶠ) * R₃₁ * ∂x_c + + (κ_symmetricᶜᶜᶠ + κ_skewᶜᶜᶠ) * R₃₂ * ∂y_c) +end +=# + +@inline function explicit_κ_∂z_c(i, j, k, grid, ::ExplicitTimeDiscretization, κ_symmetricᶜᶜᶠ, closure, buoyancy, tracers) + ∂z_c = ∂zᶜᶜᶠ(i, j, k, grid, c) + R₃₃ = isopycnal_rotation_tensor_zz_ccf(i, j, k, grid, buoyancy, tracers, closure.isopycnal_tensor) + ϵ = tapering_factorᶜᶜᶠ(i, j, k, grid, closure, tracers, buoyancy) + return ϵ * κ_symmetricᶜᶜᶠ * R₃₃ * ∂z_c +end + +@inline explicit_R₃₃_∂z_c(i, j, k, grid, ::VerticallyImplicitTimeDiscretization, args...) = zero(grid) +@inline explicit_R₃₃_∂z_u(i, j, k, grid, ::VerticallyImplicitTimeDiscretization, args...) = zero(grid) +@inline explicit_R₃₃_∂z_v(i, j, k, grid, ::VerticallyImplicitTimeDiscretization, args...) = zero(grid) + +@inline function κzᶜᶜᶠ(i, j, k, grid, closure::FlavorOfTISSD, K, ::Val{id}, clock) where id + closure = getclosure(i, j, closure) + κ_symmetric = get_tracer_κ(closure.κ_symmetric, id) + ϵ_R₃₃ = @inbounds K.ϵ_R₃₃[i, j, k] # tapered 33 component of rotation tensor + return ϵ_R₃₃ * κᶜᶜᶠ(i, j, k, grid, issd_coefficient_loc, κ_symmetric, clock) +end + +@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))") + 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..b12b92ddfc 100644 --- a/validation/mesoscale_turbulence/coarse_baroclinic_adjustment.jl +++ b/validation/isopycnal_skew_symmetric_diffusion/coarse_baroclinic_adjustment.jl @@ -4,13 +4,12 @@ using Random using Oceananigans using Oceananigans.Units using GLMakie +using Oceananigans.TurbulenceClosures: IsopycnalSkewSymmetricDiffusivity +using Oceananigans.TurbulenceClosures: TriadIsopycnalSkewSymmetricDiffusivity -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 +17,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 +28,21 @@ 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(κ_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 +82,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 +126,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/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 From 5e90efaaab29b6320bc4e99731cf60747af79f1f Mon Sep 17 00:00:00 2001 From: Gregory Wagner Date: Tue, 5 Nov 2024 15:00:23 -0700 Subject: [PATCH 02/15] Remove mwe --- mwe.jl | 13 ------------- 1 file changed, 13 deletions(-) delete mode 100644 mwe.jl diff --git a/mwe.jl b/mwe.jl deleted file mode 100644 index 212ffe6a62..0000000000 --- a/mwe.jl +++ /dev/null @@ -1,13 +0,0 @@ -using Oceananigans -using Oceananigans.Solvers: ConjugateGradientPoissonSolver, fft_poisson_solver -using Oceananigans.Models.NonhydrostaticModels: calculate_pressure_correction! - -N = 2 -x = y = (0, 1) -z = [0, 0.2, 1] -grid = RectilinearGrid(size=(N, N, N); x, y, z, halo=(2, 2, 2), topology=(Bounded, Periodic, Bounded)) -fft_solver = fft_poisson_solver(grid) -pressure_solver = ConjugateGradientPoissonSolver(grid, preconditioner=fft_poisson_solver(grid)) -model = NonhydrostaticModel(; grid, pressure_solver) -set!(model, u=1) - From 9cd2684291c3ac8930d5d9416a99dd36507deaf9 Mon Sep 17 00:00:00 2001 From: Jean-Michel Campin Date: Tue, 12 Nov 2024 11:17:27 -0500 Subject: [PATCH 03/15] get FluxTapering --- .../coarse_baroclinic_adjustment.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/validation/isopycnal_skew_symmetric_diffusion/coarse_baroclinic_adjustment.jl b/validation/isopycnal_skew_symmetric_diffusion/coarse_baroclinic_adjustment.jl index b12b92ddfc..b1eeec5693 100644 --- a/validation/isopycnal_skew_symmetric_diffusion/coarse_baroclinic_adjustment.jl +++ b/validation/isopycnal_skew_symmetric_diffusion/coarse_baroclinic_adjustment.jl @@ -6,6 +6,7 @@ using Oceananigans.Units using GLMakie using Oceananigans.TurbulenceClosures: IsopycnalSkewSymmetricDiffusivity using Oceananigans.TurbulenceClosures: TriadIsopycnalSkewSymmetricDiffusivity +using Oceananigans.TurbulenceClosures: FluxTapering gradient = "y" filename = "coarse_baroclinic_adjustment_" * gradient From 9827c4e54154180cc510ef71a77a9b52b905b7b7 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Sun, 17 Nov 2024 11:27:30 +1100 Subject: [PATCH 04/15] add reference --- .../isopycnal_skew_symmetric_diffusivity_with_triads.jl | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) 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 index a008241e1f..c2edc0f778 100644 --- 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 @@ -23,7 +23,8 @@ const issd_coefficient_loc = (Center(), Center(), Center()) κ_skew = 0, κ_symmetric = 0, isopycnal_tensor = SmallSlopeIsopycnalTensor(), - slope_limiter = FluxTapering(1e-2)) + 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 @@ -31,6 +32,12 @@ for calculating the isopycnal slopes, and (optionally) applying a `slope_limiter 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, From c34ad9105e1225ec6f64e15b94cfc55296094fde Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 19 Nov 2024 16:42:07 -0800 Subject: [PATCH 05/15] tapering --- ..._skew_symmetric_diffusivity_with_triads.jl | 74 ++++++++++++------- 1 file changed, 48 insertions(+), 26 deletions(-) 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 index c2edc0f778..fb3de07d90 100644 --- 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 @@ -134,29 +134,49 @@ end ##### 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, jx, kx, iz, jz, kz, grid, buoyancy, tracers) - bx = ∂x_b(ix, jx, kx, grid, buoyancy, tracers) - bz = ∂z_b(iz, jz, kz, grid, buoyancy, tracers) +@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(iy, jy, ky, iz, jz, kz, grid, buoyancy, tracers) - by = ∂y_b(iy, jy, ky, grid, buoyancy, tracers) - bz = ∂z_b(iz, jz, kz, grid, buoyancy, tracers) +@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, j, k, i, j, k+1, grid, buoyancy, tracers) -@inline Sx⁺⁻(i, j, k, grid, buoyancy, tracers) = triad_Sx(i+1, j, k, i, j, k, grid, buoyancy, tracers) -@inline Sx⁻⁺(i, j, k, grid, buoyancy, tracers) = triad_Sx(i, j, k, i, j, k+1, grid, buoyancy, tracers) -@inline Sx⁻⁻(i, j, k, grid, buoyancy, tracers) = triad_Sx(i, j, k, i, j, k, grid, buoyancy, tracers) +@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, k, i, j, k+1, grid, buoyancy, tracers) -@inline Sy⁺⁻(i, j, k, grid, buoyancy, tracers) = triad_Sy(i, j+1, k, i, j, k, grid, buoyancy, tracers) -@inline Sy⁻⁺(i, j, k, grid, buoyancy, tracers) = triad_Sy(i, j, k, i, j, k+1, grid, buoyancy, tracers) -@inline Sy⁻⁻(i, j, k, grid, buoyancy, tracers) = triad_Sy(i, j, k, i, j, k, grid, buoyancy, tracers) +@inline Sy⁺⁺(i, j, k, grid, buoyancy, tracers) = triad_Sy(i, j+1, k, j, k+1, grid, buoyancy, tracers) +@inline Sy⁺⁻(i, j, k, grid, buoyancy, tracers) = triad_Sy(i, j+1, k, j, k, grid, buoyancy, tracers) +@inline Sy⁻⁺(i, j, k, grid, buoyancy, tracers) = triad_Sy(i, j, k, j, k+1, grid, buoyancy, tracers) +@inline Sy⁻⁻(i, j, k, grid, buoyancy, tracers) = triad_Sy(i, j, k, j, k, grid, buoyancy, tracers) + +# We remove triads that have at least one immersed interface +@inline triad_mask_x(ix, iz, j, kx, kz, grid) = one(grid) +@inline triad_mask_y(ix, iz, j, kx, kz, grid) = one(grid) + +@inline triad_mask_x(ix, iz, j, kx, kz, grid::ImmersedBoundaryGrid) = + immersed_peripheral_node(ix, j, kx, grid) | immersed_peripheral_node(iz, j, kz, grid) + +@inline triad_mask_y(i, jy, jz, ky, kz, grid::ImmersedBoundaryGrid) = + immersed_peripheral_node(i, jy, ky, grid) | immersed_peripheral_node(i, jz, kz, grid) + +@inline ϵκx⁺⁺(i, j, k, grid, κ, clock, b, C) = triad_mask_x(i+1, i, j, k, k+1, grid) * κᶜᶜᶜ(i, j, k, grid, κ, clock) * tapering_factorᶜᶜᶜ(i, j, k, grid, b, C) +@inline ϵκx⁺⁻(i, j, k, grid, κ, clock, b, C) = triad_mask_x(i+1, i, j, k, k, grid) * κᶜᶜᶜ(i, j, k, grid, κ, clock) * tapering_factorᶜᶜᶜ(i, j, k, grid, b, C) +@inline ϵκx⁻⁺(i, j, k, grid, κ, clock, b, C) = triad_mask_x(i, i, j, k, k+1, grid) * κᶜᶜᶜ(i, j, k, grid, κ, clock) * tapering_factorᶜᶜᶜ(i, j, k, grid, b, C) +@inline ϵκx⁻⁻(i, j, k, grid, κ, clock, b, C) = triad_mask_x(i, i, j, k, k, grid) * κᶜᶜᶜ(i, j, k, grid, κ, clock) * tapering_factorᶜᶜᶜ(i, j, k, grid, b, C) + +@inline ϵκy⁺⁺(i, j, k, grid, κ, clock, b, C) = triad_mask_y(i, j+1, j, k, k+1, grid) * κᶜᶜᶜ(i, j, k, grid, κ, clock) * tapering_factorᶜᶜᶜ(i, j, k, grid, b, C) +@inline ϵκy⁺⁻(i, j, k, grid, κ, clock, b, C) = triad_mask_y(i, j+1, j, k, k, grid) * κᶜᶜᶜ(i, j, k, grid, κ, clock) * tapering_factorᶜᶜᶜ(i, j, k, grid, b, C) +@inline ϵκy⁻⁺(i, j, k, grid, κ, clock, b, C) = triad_mask_y(i, j, j, k, k+1, grid) * κᶜᶜᶜ(i, j, k, grid, κ, clock) * tapering_factorᶜᶜᶜ(i, j, k, grid, b, C) +@inline ϵκy⁻⁻(i, j, k, grid, κ, clock, b, C) = triad_mask_y(i, j, j, k, k, grid) * κᶜᶜᶜ(i, j, k, grid, κ, clock) * tapering_factorᶜᶜᶜ(i, j, k, grid, b, C) # Tensor components @inline Sxᶠᶜᶜ(i, j, k, grid, args...) = (Sx⁻⁺(i, j, k, grid, args...) + Sx⁻⁻(i, j, k, grid, args...) + # west cell, z-average @@ -174,9 +194,6 @@ end @inline Syᶜᶜᶠ(i, j, k, grid, args...) = (Sy⁺⁻(i, j, k, grid, args...) + Sy⁻⁻(i, j, k, grid, args...) + # top cell, y-average Sy⁺⁺(i, j, k-1, grid, args...) + Sy⁻⁺(i, j, k-1, grid, args...)) / 4 # bottom cell, y-average - - - struct FluxTapering{FT} max_slope :: FT end @@ -256,25 +273,30 @@ end closure = getclosure(i, j, closure) κ = get_tracer_κ(closure.κ_skew, id) - κ = κᶠᶜᶜ(i, j, k, grid, issd_coefficient_loc, κ, clock) + + ϵκ⁺⁺ = ϵκx⁺⁺(i-1, j, k, grid, κ, clock) + ϵκ⁺⁻ = ϵκx⁺⁻(i-1, j, k, grid, κ, clock) + ϵκ⁻⁺ = ϵκx⁻⁺(i, j, k, grid, κ, clock) + ϵκ⁻⁻ = ϵκx⁻⁻(i, j, k, grid, κ, clock) # Small slope approximation - R₁₁_∂x_c = ∂xᶠᶜᶜ(i, j, k, grid, c) - R₁₂_∂y_c = zero(grid) + ∂x_c = ∂xᶠᶜᶜ(i, j, k, grid, c) + + R₁₁_∂x_c = (ϵκ⁺ + ϵκ⁻) * ∂x_c / 2 # i-1 i # k+1 ------------- # | | # ┏┗ ∘ ┛┓ | k # | | - # k ------|------| + # k ------|------| - R₁₃_∂z_c = (Sx⁺⁺(i-1, j, k, grid, b, C) * ∂zᶜᶜᶠ(i-1, j, k+1, grid, c) + - Sx⁺⁻(i-1, j, k, grid, b, C) * ∂zᶜᶜᶠ(i-1, j, k, grid, c) + - Sx⁻⁺(i, j, k, grid, b, C) * ∂zᶜᶜᶠ(i, j, k+1, grid, c) + - Sx⁻⁻(i, j, k, grid, b, C) * ∂zᶜᶜᶠ(i, j, k, grid, c)) / 4 + 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 - κ * (R₁₁_∂x_c + R₁₂_∂y_c + R₁₃_∂z_c) + return - (R₁₁_∂x_c + R₁₃_∂z_c) end # defined at cfc From 0ebdb499d613093ad2717dbbe2cf9e05dd6428c4 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 19 Nov 2024 17:15:10 -0800 Subject: [PATCH 06/15] continue formulation --- ..._skew_symmetric_diffusivity_with_triads.jl | 80 +++++++++---------- 1 file changed, 37 insertions(+), 43 deletions(-) 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 index fb3de07d90..5df5bfda42 100644 --- 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 @@ -178,22 +178,6 @@ end @inline ϵκy⁻⁺(i, j, k, grid, κ, clock, b, C) = triad_mask_y(i, j, j, k, k+1, grid) * κᶜᶜᶜ(i, j, k, grid, κ, clock) * tapering_factorᶜᶜᶜ(i, j, k, grid, b, C) @inline ϵκy⁻⁻(i, j, k, grid, κ, clock, b, C) = triad_mask_y(i, j, j, k, k, grid) * κᶜᶜᶜ(i, j, k, grid, κ, clock) * tapering_factorᶜᶜᶜ(i, j, k, grid, b, C) -# Tensor components -@inline Sxᶠᶜᶜ(i, j, k, grid, args...) = (Sx⁻⁺(i, j, k, grid, args...) + Sx⁻⁻(i, j, k, grid, args...) + # west cell, z-average - Sx⁻⁺(i-1, j, k, grid, args...) + Sx⁻⁻(i-1, j, k, grid, args...)) / 4 # east cell, z-average - -@inline Syᶜᶠᶜ(i, j, k, grid, args...) = (Sy⁻⁺(i, j, k, grid, args...) + Sy⁻⁻(i, j, k, grid, args...) + # south cell, z-average - Sy⁻⁺(i, j-1, k, grid, args...) + Sy⁻⁻(i, j-1, k, grid, args...)) / 4 # north cell, z-average - -@inline Syᶜᶜᶠ(i, j, k, grid, args...) = (Sy⁺⁻(i, j, k, grid, args...) + Sy⁻⁻(i, j, k, grid, args...) + # top cell, y-average - Sy⁺⁺(i, j, k-1, grid, args...) + Sy⁻⁺(i, j, k-1, grid, args...)) / 4 # bottom cell, y-average - -@inline Sxᶜᶜᶠ(i, j, k, grid, args...) = (Sx⁺⁻(i, j, k, grid, args...) + Sx⁻⁻(i, j, k, grid, args...) + # top cell, x-average - Sx⁺⁺(i, j, k-1, grid, args...) + Sx⁻⁺(i, j, k-1, grid, args...)) / 4 # bottom cell, x-average - -@inline Syᶜᶜᶠ(i, j, k, grid, args...) = (Sy⁺⁻(i, j, k, grid, args...) + Sy⁻⁻(i, j, k, grid, args...) + # top cell, y-average - Sy⁺⁺(i, j, k-1, grid, args...) + Sy⁻⁺(i, j, k-1, grid, args...)) / 4 # bottom cell, y-average - struct FluxTapering{FT} max_slope :: FT end @@ -274,15 +258,13 @@ end κ = get_tracer_κ(closure.κ_skew, id) - ϵκ⁺⁺ = ϵκx⁺⁺(i-1, j, k, grid, κ, clock) - ϵκ⁺⁻ = ϵκx⁺⁻(i-1, j, k, grid, κ, clock) - ϵκ⁻⁺ = ϵκx⁻⁺(i, j, k, grid, κ, clock) - ϵκ⁻⁻ = ϵκx⁻⁻(i, j, k, grid, κ, clock) + ϵκ⁺⁺ = ϵκx⁺⁺(i-1, j, k, grid, κ, clock, b, C) + ϵκ⁺⁻ = ϵκx⁺⁻(i-1, j, k, grid, κ, clock, b, C) + ϵκ⁻⁺ = ϵκx⁻⁺(i, j, k, grid, κ, clock, b, C) + ϵκ⁻⁻ = ϵκx⁻⁻(i, j, k, grid, κ, clock, b, C) # Small slope approximation ∂x_c = ∂xᶠᶜᶜ(i, j, k, grid, c) - - R₁₁_∂x_c = (ϵκ⁺ + ϵκ⁻) * ∂x_c / 2 # i-1 i # k+1 ------------- @@ -296,7 +278,7 @@ end ϵκ⁻⁺ * (∂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 - (R₁₁_∂x_c + R₁₃_∂z_c) + return - Fx end # defined at cfc @@ -306,18 +288,20 @@ end closure = getclosure(i, j, closure) κ = get_tracer_κ(closure.κ_skew, id) - κ = κᶜᶠᶜ(i, j, k, grid, issd_coefficient_loc, κ, clock) - # Small slope approximation - R₂₁_∂x_c = zero(grid) - R₂₂_∂y_c = ∂yᶜᶠᶜ(i, j, k, grid, c) + ∂y_c = ∂yᶜᶠᶜ(i, j, k, grid, c) - R₂₃_∂z_c = (Sy⁺⁺(i, j-1, k, grid, b, C) * ∂zᶜᶜᶠ(i, j-1, k+1, grid, c) + - Sy⁺⁻(i, j-1, k, grid, b, C) * ∂zᶜᶜᶠ(i, j-1, k, grid, c) + - Sy⁻⁺(i, j, k, grid, b, C) * ∂zᶜᶜᶠ(i, j, k+1, grid, c) + - Sy⁻⁻(i, j, k, grid, b, C) * ∂zᶜᶜᶠ(i, j, k, grid, c)) / 4 + ϵκ⁺⁺ = ϵκy⁺⁺(i, j-1, k, grid, κ, clock, b, C) + ϵκ⁺⁻ = ϵκy⁺⁻(i, j-1, k, grid, κ, clock, b, C) + ϵκ⁻⁺ = ϵκy⁻⁺(i, j, k, grid, κ, clock, b, C) + ϵκ⁻⁻ = ϵκy⁻⁻(i, j, k, grid, κ, 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 - κ * (R₂₁_∂x_c + R₂₂_∂y_c + R₂₃_∂z_c) + return - Fy end # defined at ccf @@ -327,7 +311,17 @@ end closure = getclosure(i, j, closure) κ = get_tracer_κ(closure.κ_skew, id) - κ = κᶜᶜᶠ(i, j, k, grid, issd_coefficient_loc, κ, clock) + + ϵκˣ⁻⁻ = ϵκx⁻⁻(i, j, k, grid, κ, clock, b, C) + ϵκˣ⁺⁻ = ϵκx⁺⁻(i, j, k, grid, κ, clock, b, C) + ϵκˣ⁻⁺ = ϵκx⁻⁺(i, j, k-1, grid, κ, clock, b, C) + ϵκˣ⁺⁺ = ϵκx⁺⁺(i, j, k-1, grid, κ, clock, b, C) + + ϵκʸ⁻⁻ = ϵκʸ⁻⁻(i, j, k, grid, κ, clock, b, C) + ϵκʸ⁺⁻ = ϵκʸ⁺⁻(i, j, k, grid, κ, clock, b, C) + ϵκʸ⁻⁺ = ϵκʸ⁻⁺(i, j, k-1, grid, κ, clock, b, C) + ϵκʸ⁺⁺ = ϵκʸ⁺⁺(i, j, k-1, grid, κ, clock, b, C) + # Triad diagram: # @@ -342,19 +336,19 @@ end # | | | | # -------------------- - 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₃₁_∂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₃₂_∂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) + κϵ_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) + return - κR₃₁_∂x_c - κR₃₂_∂y_c - κϵ_R₃₃_∂z_c end @inline function explicit_R₃₃_∂z_c(i, j, k, grid, ::ExplicitTimeDiscretization, c, closure, b, C) From 0db81610c1007e41b4ea8bbc78f8986865dd5203 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 19 Nov 2024 17:49:03 -0800 Subject: [PATCH 07/15] still need to complete the tapering --- ..._skew_symmetric_diffusivity_with_triads.jl | 278 ++++++------------ 1 file changed, 93 insertions(+), 185 deletions(-) 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 index 5df5bfda42..76174c9d72 100644 --- 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 @@ -60,11 +60,8 @@ end TriadIsopycnalSkewSymmetricDiffusivity(FT::DataType; kw...) = TriadIsopycnalSkewSymmetricDiffusivity(VerticallyImplicitTimeDiscretization(), FT; kw...) -function with_tracers(tracers, closure::TISSD{TD, N}) where {TD, N} - κ_skew = !isa(closure.κ_skew, NamedTuple) ? closure.κ_skew : tracer_diffusivities(tracers, closure.κ_skew) - κ_symmetric = !isa(closure.κ_symmetric, NamedTuple) ? closure.κ_symmetric : tracer_diffusivities(tracers, closure.κ_symmetric) - return TriadIsopycnalSkewSymmetricDiffusivity{TD, N}(κ_skew, κ_symmetric, closure.isopycnal_tensor, closure.slope_limiter) -end +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) @@ -81,37 +78,39 @@ function with_tracers(tracers, closure_vector::TISSDVector) end # Note: computing diffusivities at cell centers for now. -function DiffusivityFields(grid, tracer_names, bcs, closure::FlavorOfTISSD{TD}) where TD +function DiffusivityFields(grid, tracer_names, bcs, ::FlavorOfTISSD{TD}) where TD if TD() isa VerticallyImplicitTimeDiscretization # Precompute the _tapered_ 33 component of the isopycnal rotation tensor - return (; ϵ_R₃₃ = Field((Center, Center, Face), grid)) + K = (; ϵκR₃₃ = ZFaceField(grid)) else return nothing end + + return with_tracers(tracer_names, 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.ϵ_R₃₃, grid, closure, tracers, buoyancy) + diffusivities, grid, closure, clock, buoyancy, tracers) end return nothing end -@kernel function triad_compute_tapered_R₃₃!(ϵ_R₃₃, grid, closure, tracers, buoyancy) +@kernel function triad_compute_tapered_R₃₃!(K, grid, closure, clock, b, C) i, j, k, = @index(Global, NTuple) closure = getclosure(i, j, closure) - R₃₃ = isopycnal_rotation_tensor_zz_ccf(i, j, k, grid, buoyancy, tracers, closure.isopycnal_tensor) - ϵ = tapering_factorᶜᶜᶠ(i, j, k, grid, closure, tracers, buoyancy) - @inbounds ϵ_R₃₃[i, j, k] = ϵ * R₃₃ + κ = closure.κ_symmetric + @inbounds K.ϵκR₃₃[i, j, k] = ϵκR₃₃(i, j, k, grid, κ, clock, b, C) end ##### @@ -182,65 +181,11 @@ struct FluxTapering{FT} max_slope :: FT end -""" - 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 - # Diffusive fluxes @inline get_tracer_κ(κ::NamedTuple, tracer_index) = @inbounds κ[tracer_index] @inline get_tracer_κ(κ, tracer_index) = κ - # Triad diagram key # ================= # @@ -255,8 +200,7 @@ end c, clock, C, b) where id closure = getclosure(i, j, closure) - - κ = get_tracer_κ(closure.κ_skew, id) + κ = closure.κ_symmetric ϵκ⁺⁺ = ϵκx⁺⁺(i-1, j, k, grid, κ, clock, b, C) ϵκ⁺⁻ = ϵκx⁺⁻(i-1, j, k, grid, κ, clock, b, C) @@ -286,8 +230,7 @@ end c, clock, C, b) where id closure = getclosure(i, j, closure) - - κ = get_tracer_κ(closure.κ_skew, id) + κ = closure.κ_symmetric ∂y_c = ∂yᶜᶠᶜ(i, j, k, grid, c) @@ -309,8 +252,7 @@ end c, clock, C, b) where {TD, id} closure = getclosure(i, j, closure) - - κ = get_tracer_κ(closure.κ_skew, id) + κ = closure.κ_symmetric ϵκˣ⁻⁻ = ϵκx⁻⁻(i, j, k, grid, κ, clock, b, C) ϵκˣ⁺⁻ = ϵκx⁺⁻(i, j, k, grid, κ, clock, b, C) @@ -322,7 +264,6 @@ end ϵκʸ⁻⁺ = ϵκʸ⁻⁺(i, j, k-1, grid, κ, clock, b, C) ϵκʸ⁺⁺ = ϵκʸ⁺⁺(i, j, k-1, grid, κ, clock, b, C) - # Triad diagram: # # i-1 i i+1 @@ -337,136 +278,49 @@ end # -------------------- κ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 + ϵκˣ⁺⁻ * 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 + ϵκʸ⁺⁻ * 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) + κϵ_R₃₃_∂z_c = explicit_R₃₃_∂z_c(i, j, k, grid, TD(), c, ) return - κR₃₁_∂x_c - κR₃₂_∂y_c - κϵ_R₃₃_∂z_c end -@inline function explicit_R₃₃_∂z_c(i, j, k, grid, ::ExplicitTimeDiscretization, c, closure, b, C) - ϵ = tapering_factorᶜᶜᶠ(i, j, k, grid, closure, C, b) - ∂z_c = ∂zᶜᶜᶠ(i, j, k, grid, c) - Sx = Sxᶜᶜᶠ(i, j, k, grid, b, C) - Sy = Syᶜᶜᶠ(i, j, k, grid, b, C) - R₃₃ = Sx^2 + Sy^2 - return ϵ * R₃₃ * ∂z_c -end - -#= -# defined at fcc -@inline function diffusive_flux_x(i, j, k, grid, - closure::Union{TISSD, TISSDVector}, diffusivity_fields, ::Val{tracer_index}, - c, clock, fields, buoyancy) where tracer_index - - closure = getclosure(i, j, closure) - - κ_skew = get_tracer_κ(closure.κ_skew, tracer_index) - κ_symmetric = get_tracer_κ(closure.κ_symmetric, tracer_index) - - κ_skewᶠᶜᶜ = κᶠᶜᶜ(i, j, k, grid, issd_coefficient_loc, κ_skew, clock) - κ_symmetricᶠᶜᶜ = κᶠᶜᶜ(i, j, k, grid, issd_coefficient_loc, κ_symmetric, clock) - - ∂x_c = ∂xᶠᶜᶜ(i, j, k, grid, c) - - # Average... of... the gradient! - ∂y_c = ℑxyᶠᶜᵃ(i, j, k, grid, ∂yᶜᶠᶜ, c) - ∂z_c = ℑxzᶠᵃᶜ(i, j, k, grid, ∂zᶜᶜᶠ, c) - - R₁₁ = one(grid) - R₁₂ = zero(grid) - R₁₃ = isopycnal_rotation_tensor_xz_fcc(i, j, k, grid, buoyancy, fields, closure.isopycnal_tensor) - - ϵ = tapering_factorᶠᶜᶜ(i, j, k, grid, closure, fields, buoyancy) - - return - ϵ * ( κ_symmetricᶠᶜᶜ * R₁₁ * ∂x_c + - κ_symmetricᶠᶜᶜ * R₁₂ * ∂y_c + - (κ_symmetricᶠᶜᶜ - κ_skewᶠᶜᶜ) * R₁₃ * ∂z_c) -end - -# defined at cfc -@inline function diffusive_flux_y(i, j, k, grid, - closure::Union{TISSD, TISSDVector}, diffusivity_fields, ::Val{tracer_index}, - c, clock, fields, buoyancy) where tracer_index - - closure = getclosure(i, j, closure) - - κ_skew = get_tracer_κ(closure.κ_skew, tracer_index) - κ_symmetric = get_tracer_κ(closure.κ_symmetric, tracer_index) +@inline function ϵκR₃₃(i, j, k, grid, κ, clock, b, C) - κ_skewᶜᶠᶜ = κᶜᶠᶜ(i, j, k, grid, issd_coefficient_loc, κ_skew, clock) - κ_symmetricᶜᶠᶜ = κᶜᶠᶜ(i, j, k, grid, issd_coefficient_loc, κ_symmetric, clock) - - ∂y_c = ∂yᶜᶠᶜ(i, j, k, grid, c) - - # Average... of... the gradient! - ∂x_c = ℑxyᶜᶠᵃ(i, j, k, grid, ∂xᶠᶜᶜ, c) - ∂z_c = ℑyzᵃᶠᶜ(i, j, k, grid, ∂zᶜᶜᶠ, c) + ϵκˣ⁻⁻ = ϵκx⁻⁻(i, j, k, grid, κ, clock, b, C) + ϵκˣ⁺⁻ = ϵκx⁺⁻(i, j, k, grid, κ, clock, b, C) + ϵκˣ⁻⁺ = ϵκx⁻⁺(i, j, k-1, grid, κ, clock, b, C) + ϵκˣ⁺⁺ = ϵκx⁺⁺(i, j, k-1, grid, κ, clock, b, C) - R₂₁ = zero(grid) - R₂₂ = one(grid) - R₂₃ = isopycnal_rotation_tensor_yz_cfc(i, j, k, grid, buoyancy, fields, closure.isopycnal_tensor) + ϵκʸ⁻⁻ = ϵκʸ⁻⁻(i, j, k, grid, κ, clock, b, C) + ϵκʸ⁺⁻ = ϵκʸ⁺⁻(i, j, k, grid, κ, clock, b, C) + ϵκʸ⁻⁺ = ϵκʸ⁻⁺(i, j, k-1, grid, κ, clock, b, C) + ϵκʸ⁺⁺ = ϵκʸ⁺⁺(i, j, k-1, grid, κ, clock, b, C) - ϵ = tapering_factorᶜᶠᶜ(i, j, k, grid, closure, fields, buoyancy) + ϵκ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 - ϵ * (κ_symmetricᶜᶠᶜ * R₂₁ * ∂x_c + - κ_symmetricᶜᶠᶜ * R₂₂ * ∂y_c + - (κ_symmetricᶜᶠᶜ - κ_skewᶜᶠᶜ) * R₂₃ * ∂z_c) + return ϵκR₃₃ end -# defined at ccf -@inline function diffusive_flux_z(i, j, k, grid, - closure::FlavorOfTISSD{TD}, diffusivity_fields, ::Val{tracer_index}, - c, clock, fields, buoyancy) where {tracer_index, TD} - - closure = getclosure(i, j, closure) - - κ_skew = get_tracer_κ(closure.κ_skew, tracer_index) - κ_symmetric = get_tracer_κ(closure.κ_symmetric, tracer_index) - - κ_skewᶜᶜᶠ = κᶜᶜᶠ(i, j, k, grid, issd_coefficient_loc, κ_skew, clock) - κ_symmetricᶜᶜᶠ = κᶜᶜᶠ(i, j, k, grid, issd_coefficient_loc, κ_symmetric, clock) +@inline explicit_R₃₃_∂z_c(i, j, k, grid, ::ExplicitTimeDiscretization, c, closure, b, C) = ϵκR₃₃(i, j, k, grid, κ, clock, b, C) * ∂zᶜᶜᶠ(i, j, k, grid, c) - # Average... of... the gradient! - ∂x_c = ℑxzᶜᵃᶠ(i, j, k, grid, ∂xᶠᶜᶜ, c) - ∂y_c = ℑyzᵃᶜᶠ(i, j, k, grid, ∂yᶜᶠᶜ, c) - - R₃₁ = isopycnal_rotation_tensor_xz_ccf(i, j, k, grid, buoyancy, fields, closure.isopycnal_tensor) - R₃₂ = isopycnal_rotation_tensor_yz_ccf(i, j, k, grid, buoyancy, fields, closure.isopycnal_tensor) - - κ_symmetric_∂z_c = explicit_κ_∂z_c(i, j, k, grid, TD(), c, κ_symmetricᶜᶜᶠ, closure, buoyancy, fields) - - ϵ = tapering_factorᶜᶜᶠ(i, j, k, grid, closure, fields, buoyancy) - - return - ϵ * κ_symmetric_∂z_c - ϵ * ((κ_symmetricᶜᶜᶠ + κ_skewᶜᶜᶠ) * R₃₁ * ∂x_c + - (κ_symmetricᶜᶜᶠ + κ_skewᶜᶜᶠ) * R₃₂ * ∂y_c) -end -=# - -@inline function explicit_κ_∂z_c(i, j, k, grid, ::ExplicitTimeDiscretization, κ_symmetricᶜᶜᶠ, closure, buoyancy, tracers) - ∂z_c = ∂zᶜᶜᶠ(i, j, k, grid, c) - R₃₃ = isopycnal_rotation_tensor_zz_ccf(i, j, k, grid, buoyancy, tracers, closure.isopycnal_tensor) - ϵ = tapering_factorᶜᶜᶠ(i, j, k, grid, closure, tracers, buoyancy) - return ϵ * κ_symmetricᶜᶜᶠ * R₃₃ * ∂z_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, args...) = zero(grid) @inline explicit_R₃₃_∂z_u(i, j, k, grid, ::VerticallyImplicitTimeDiscretization, args...) = zero(grid) @inline explicit_R₃₃_∂z_v(i, j, k, grid, ::VerticallyImplicitTimeDiscretization, args...) = zero(grid) -@inline function κzᶜᶜᶠ(i, j, k, grid, closure::FlavorOfTISSD, K, ::Val{id}, clock) where id - closure = getclosure(i, j, closure) - κ_symmetric = get_tracer_κ(closure.κ_symmetric, id) - ϵ_R₃₃ = @inbounds K.ϵ_R₃₃[i, j, k] # tapered 33 component of rotation tensor - return ϵ_R₃₃ * κᶜᶜᶠ(i, j, k, grid, issd_coefficient_loc, κ_symmetric, clock) -end +@inline κzᶜᶜᶠ(i, j, k, grid, closure::FlavorOfTISSD, K, ::Val{id}, clock) where id = @inbounds K[id].ϵκ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) @@ -493,4 +347,58 @@ 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 From 6099b02ee4a80238e60fca86a94b3f14e0693f4d Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 19 Nov 2024 18:06:21 -0800 Subject: [PATCH 08/15] runs but probably wrong --- .../abstract_scalar_diffusivity_closure.jl | 6 +- ..._skew_symmetric_diffusivity_with_triads.jl | 82 ++++++++++--------- .../coarse_baroclinic_adjustment.jl | 3 +- 3 files changed, 49 insertions(+), 42 deletions(-) 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_with_triads.jl b/src/TurbulenceClosures/turbulence_closure_implementations/isopycnal_skew_symmetric_diffusivity_with_triads.jl index 76174c9d72..1b278bb658 100644 --- 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 @@ -16,7 +16,6 @@ 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 -const issd_coefficient_loc = (Center(), Center(), Center()) """ TriadIsopycnalSkewSymmetricDiffusivity([time_disc=VerticallyImplicitTimeDiscretization(), FT=Float64;] @@ -86,7 +85,7 @@ function DiffusivityFields(grid, tracer_names, bcs, ::FlavorOfTISSD{TD}) where T return nothing end - return with_tracers(tracer_names, K) + return K end function compute_diffusivities!(diffusivities, closure::FlavorOfTISSD{TD}, model; parameters = :xyz) where TD @@ -167,15 +166,15 @@ end @inline triad_mask_y(i, jy, jz, ky, kz, grid::ImmersedBoundaryGrid) = immersed_peripheral_node(i, jy, ky, grid) | immersed_peripheral_node(i, jz, kz, grid) -@inline ϵκx⁺⁺(i, j, k, grid, κ, clock, b, C) = triad_mask_x(i+1, i, j, k, k+1, grid) * κᶜᶜᶜ(i, j, k, grid, κ, clock) * tapering_factorᶜᶜᶜ(i, j, k, grid, b, C) -@inline ϵκx⁺⁻(i, j, k, grid, κ, clock, b, C) = triad_mask_x(i+1, i, j, k, k, grid) * κᶜᶜᶜ(i, j, k, grid, κ, clock) * tapering_factorᶜᶜᶜ(i, j, k, grid, b, C) -@inline ϵκx⁻⁺(i, j, k, grid, κ, clock, b, C) = triad_mask_x(i, i, j, k, k+1, grid) * κᶜᶜᶜ(i, j, k, grid, κ, clock) * tapering_factorᶜᶜᶜ(i, j, k, grid, b, C) -@inline ϵκx⁻⁻(i, j, k, grid, κ, clock, b, C) = triad_mask_x(i, i, j, k, k, grid) * κᶜᶜᶜ(i, j, k, grid, κ, 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+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, κ, clock, b, C) = triad_mask_y(i, j+1, j, k, k+1, grid) * κᶜᶜᶜ(i, j, k, grid, κ, clock) * tapering_factorᶜᶜᶜ(i, j, k, grid, b, C) -@inline ϵκy⁺⁻(i, j, k, grid, κ, clock, b, C) = triad_mask_y(i, j+1, j, k, k, grid) * κᶜᶜᶜ(i, j, k, grid, κ, clock) * tapering_factorᶜᶜᶜ(i, j, k, grid, b, C) -@inline ϵκy⁻⁺(i, j, k, grid, κ, clock, b, C) = triad_mask_y(i, j, j, k, k+1, grid) * κᶜᶜᶜ(i, j, k, grid, κ, clock) * tapering_factorᶜᶜᶜ(i, j, k, grid, b, C) -@inline ϵκy⁻⁻(i, j, k, grid, κ, clock, b, C) = triad_mask_y(i, j, j, k, k, grid) * κᶜᶜᶜ(i, j, k, grid, κ, 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) struct FluxTapering{FT} max_slope :: FT @@ -202,10 +201,12 @@ end closure = getclosure(i, j, closure) κ = closure.κ_symmetric - ϵκ⁺⁺ = ϵκx⁺⁺(i-1, j, k, grid, κ, clock, b, C) - ϵκ⁺⁻ = ϵκx⁺⁻(i-1, j, k, grid, κ, clock, b, C) - ϵκ⁻⁺ = ϵκx⁻⁺(i, j, k, grid, κ, clock, b, C) - ϵκ⁻⁻ = ϵκx⁻⁻(i, j, k, grid, κ, clock, b, C) + 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) @@ -232,12 +233,14 @@ end 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, κ, clock, b, C) - ϵκ⁺⁻ = ϵκy⁺⁻(i, j-1, k, grid, κ, clock, b, C) - ϵκ⁻⁺ = ϵκy⁻⁺(i, j, k, grid, κ, clock, b, C) - ϵκ⁻⁻ = ϵκy⁻⁻(i, j, k, grid, κ, clock, b, 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)) + @@ -254,15 +257,17 @@ end closure = getclosure(i, j, closure) κ = closure.κ_symmetric - ϵκˣ⁻⁻ = ϵκx⁻⁻(i, j, k, grid, κ, clock, b, C) - ϵκˣ⁺⁻ = ϵκx⁺⁻(i, j, k, grid, κ, clock, b, C) - ϵκˣ⁻⁺ = ϵκx⁻⁺(i, j, k-1, grid, κ, clock, b, C) - ϵκˣ⁺⁺ = ϵκx⁺⁺(i, j, k-1, grid, κ, clock, b, C) + loc = (Center(), Center(), Center()) - ϵκʸ⁻⁻ = ϵκʸ⁻⁻(i, j, k, grid, κ, clock, b, C) - ϵκʸ⁺⁻ = ϵκʸ⁺⁻(i, j, k, grid, κ, clock, b, C) - ϵκʸ⁻⁺ = ϵκʸ⁻⁺(i, j, k-1, grid, κ, clock, b, C) - ϵκʸ⁺⁺ = ϵκʸ⁺⁺(i, j, k-1, grid, κ, clock, b, C) + ϵκˣ⁻⁻ = ϵκ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: # @@ -287,22 +292,23 @@ 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, ) + κϵ_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, κ, clock, b, C) - ϵκˣ⁺⁻ = ϵκx⁺⁻(i, j, k, grid, κ, clock, b, C) - ϵκˣ⁻⁺ = ϵκx⁻⁺(i, j, k-1, grid, κ, clock, b, C) - ϵκˣ⁺⁺ = ϵκx⁺⁺(i, j, k-1, grid, κ, clock, b, C) + ϵκˣ⁻⁻ = ϵκ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) - ϵκʸ⁻⁻ = ϵκʸ⁻⁻(i, j, k, grid, κ, clock, b, C) - ϵκʸ⁺⁻ = ϵκʸ⁺⁻(i, j, k, grid, κ, clock, b, C) - ϵκʸ⁻⁺ = ϵκʸ⁻⁺(i, j, k-1, grid, κ, clock, b, C) - ϵκʸ⁺⁺ = ϵκʸ⁺⁺(i, j, k-1, grid, κ, 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 + @@ -313,14 +319,10 @@ end end @inline explicit_R₃₃_∂z_c(i, j, k, grid, ::ExplicitTimeDiscretization, c, closure, b, C) = ϵκR₃₃(i, j, k, grid, κ, clock, b, C) * ∂zᶜᶜᶠ(i, j, k, grid, c) - @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, args...) = zero(grid) -@inline explicit_R₃₃_∂z_u(i, j, k, grid, ::VerticallyImplicitTimeDiscretization, args...) = zero(grid) -@inline explicit_R₃₃_∂z_v(i, j, k, grid, ::VerticallyImplicitTimeDiscretization, args...) = zero(grid) -@inline κzᶜᶜᶠ(i, j, k, grid, closure::FlavorOfTISSD, K, ::Val{id}, clock) where id = @inbounds K[id].ϵκR₃₃[i, j, k] +@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) diff --git a/validation/isopycnal_skew_symmetric_diffusion/coarse_baroclinic_adjustment.jl b/validation/isopycnal_skew_symmetric_diffusion/coarse_baroclinic_adjustment.jl index b1eeec5693..96c5ee1d3a 100644 --- a/validation/isopycnal_skew_symmetric_diffusion/coarse_baroclinic_adjustment.jl +++ b/validation/isopycnal_skew_symmetric_diffusion/coarse_baroclinic_adjustment.jl @@ -30,7 +30,8 @@ grid = RectilinearGrid(topology = (Periodic, Bounded, Bounded), coriolis = FPlane(latitude = -45) gerdes_koberle_willebrand_tapering = FluxTapering(1e-2) -triad_closure = TriadIsopycnalSkewSymmetricDiffusivity(κ_skew = 1e3, +triad_closure = TriadIsopycnalSkewSymmetricDiffusivity(VerticallyImplicitTimeDiscretization(), + κ_skew = 1e3, κ_symmetric = 1e3, slope_limiter = gerdes_koberle_willebrand_tapering) From 12363a6972fbdd944184ebfbda88feaade3c5e59 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 19 Nov 2024 18:12:14 -0800 Subject: [PATCH 09/15] explicit runs --- .../isopycnal_skew_symmetric_diffusivity_with_triads.jl | 7 +++++-- .../coarse_baroclinic_adjustment.jl | 6 +++--- 2 files changed, 8 insertions(+), 5 deletions(-) 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 index 1b278bb658..9b2470d7f2 100644 --- 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 @@ -318,9 +318,12 @@ end return ϵκR₃₃ end -@inline explicit_R₃₃_∂z_c(i, j, k, grid, ::ExplicitTimeDiscretization, c, closure, b, C) = ϵκR₃₃(i, j, k, grid, κ, clock, b, C) * ∂zᶜᶜᶠ(i, j, k, grid, c) -@inline explicit_R₃₃_∂z_c(i, j, k, grid, ::VerticallyImplicitTimeDiscretization, c, closure, b, C) = zero(grid) +@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] diff --git a/validation/isopycnal_skew_symmetric_diffusion/coarse_baroclinic_adjustment.jl b/validation/isopycnal_skew_symmetric_diffusion/coarse_baroclinic_adjustment.jl index 96c5ee1d3a..e72aa85464 100644 --- a/validation/isopycnal_skew_symmetric_diffusion/coarse_baroclinic_adjustment.jl +++ b/validation/isopycnal_skew_symmetric_diffusion/coarse_baroclinic_adjustment.jl @@ -30,12 +30,12 @@ grid = RectilinearGrid(topology = (Periodic, Bounded, Bounded), coriolis = FPlane(latitude = -45) gerdes_koberle_willebrand_tapering = FluxTapering(1e-2) -triad_closure = TriadIsopycnalSkewSymmetricDiffusivity(VerticallyImplicitTimeDiscretization(), - κ_skew = 1e3, +triad_closure = TriadIsopycnalSkewSymmetricDiffusivity(κ_skew = 1e3, κ_symmetric = 1e3, slope_limiter = gerdes_koberle_willebrand_tapering) -cox_closure = IsopycnalSkewSymmetricDiffusivity(κ_skew = 1e3, +cox_closure = IsopycnalSkewSymmetricDiffusivity(VerticallyImplicitTimeDiscretization(), + κ_skew = 1e3, κ_symmetric = 1e3, slope_limiter = gerdes_koberle_willebrand_tapering) From f97cc4a943bf50dfa73239f7cac9f90dc0249059 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Thu, 21 Nov 2024 12:02:02 -0800 Subject: [PATCH 10/15] remove duplicate --- .../isopycnal_skew_symmetric_diffusivity_with_triads.jl | 4 ---- 1 file changed, 4 deletions(-) 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 index 9b2470d7f2..48d6f4e7c8 100644 --- 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 @@ -176,10 +176,6 @@ end @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) -struct FluxTapering{FT} - max_slope :: FT -end - # Diffusive fluxes @inline get_tracer_κ(κ::NamedTuple, tracer_index) = @inbounds κ[tracer_index] From 7680994d0b790d65066897e7a26cb4970f852bbb Mon Sep 17 00:00:00 2001 From: Jean-Michel Campin Date: Thu, 21 Nov 2024 19:04:43 -0500 Subject: [PATCH 11/15] Try a deeper test, isopycnal diffusion only --- .../new_front_relax.jl | 203 ++++++++++++++++++ 1 file changed, 203 insertions(+) create mode 100644 validation/isopycnal_skew_symmetric_diffusion/new_front_relax.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..404fa7e3d5 --- /dev/null +++ b/validation/isopycnal_skew_symmetric_diffusion/new_front_relax.jl @@ -0,0 +1,203 @@ +using Printf +using Statistics +using Random +using Oceananigans +using Oceananigans.Units +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 = 2days +save_fields_interval = 15minutes +stop_time = 1hours +Δt = 15minutes + +# viscosity +viscAh=300. +viscAz=2.e-4 + +grid = RectilinearGrid(topology = (Periodic, Bounded, Bounded), + size = (Nx, Ny, Nz), + x = (-Ly/2, Ly/2), + y = (-Ly/2, Ly/2), + z = ([0 cumsum(reverse(dz))']'.-Lz), + halo = (3, 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(κ_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." + +""" +Linear ramp from 0 to 1 between -Δy/2 and +Δy/2. + +For example: + +y < y₀ => ramp = 0 +y₀ < y < y₀ + Δy => ramp = y / Δy +y > y₀ + Δy => ramp = 1 +""" +function ramp(x, y, Δ) + gradient == "x" && return min(max(0, x / Δ + 1/2), 1) + gradient == "y" && return min(max(0, y / Δ + 1/2), 1) +end + +# 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=z[6] +Δy = Ly/8 +Δz = 500 + +#Δc = 2Δy +Δb = Ly/Ny * M² +ϵb = 1e-2 * Δb # noise amplitude + +# bᵢ(x, y, z) = N² * z + Δb * ramp(x, y, Δy) +# cᵢ(x, y, z) = exp(-y^2 / 2Δc^2) * exp(-(z + Lz/4)^2 / 2Δz^2) +bᵢ(x, y, z) = N² * z - M² * Ly*sin(pi*y/Ly) *exp(-(3*z/Lz)^2) +cᵢ(x, 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) + +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) + +if gradient == "y" # average in x + # un(n) = interior(ut[n], 1, :, :) + # bn(n) = interior(bt[n], 1, :, :) + # cn(n) = interior(ct[n], 1, :, :) + 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, :, :) +else # average in y + un(n) = interior(mean(ut[n], dims=2), :, 1, :) + bn(n) = interior(mean(bt[n], dims=2), :, 1, :) + cn(n) = interior(mean(ct[n], dims=2), :, 1, :) +end + +@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 + From f6b9cd6ad357f597f481802965d0284e99a53066 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Thu, 21 Nov 2024 16:08:17 -0800 Subject: [PATCH 12/15] remove extra definition --- .../isopycnal_skew_symmetric_diffusivity_with_triads.jl | 5 ----- 1 file changed, 5 deletions(-) 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 index 48d6f4e7c8..fa17b4ec2b 100644 --- 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 @@ -176,11 +176,6 @@ end @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) -# Diffusive fluxes - -@inline get_tracer_κ(κ::NamedTuple, tracer_index) = @inbounds κ[tracer_index] -@inline get_tracer_κ(κ, tracer_index) = κ - # Triad diagram key # ================= # From 07315c23726292d810b67024653d41fc87a84e6f Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Thu, 21 Nov 2024 16:38:27 -0800 Subject: [PATCH 13/15] index fixes --- .../isopycnal_skew_symmetric_diffusivity.jl | 1 - ..._skew_symmetric_diffusivity_with_triads.jl | 10 ++-- .../new_front_relax.jl | 47 +++++++------------ 3 files changed, 21 insertions(+), 37 deletions(-) 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 index fa17b4ec2b..bc8eab460e 100644 --- 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 @@ -151,10 +151,10 @@ end @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, k, j, k+1, grid, buoyancy, tracers) -@inline Sy⁺⁻(i, j, k, grid, buoyancy, tracers) = triad_Sy(i, j+1, k, j, k, grid, buoyancy, tracers) -@inline Sy⁻⁺(i, j, k, grid, buoyancy, tracers) = triad_Sy(i, j, k, j, k+1, grid, buoyancy, tracers) -@inline Sy⁻⁻(i, j, k, grid, buoyancy, tracers) = triad_Sy(i, j, k, j, 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 have at least one immersed interface @inline triad_mask_x(ix, iz, j, kx, kz, grid) = one(grid) @@ -237,7 +237,7 @@ end ϵκ⁺⁻ * (∂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 diff --git a/validation/isopycnal_skew_symmetric_diffusion/new_front_relax.jl b/validation/isopycnal_skew_symmetric_diffusion/new_front_relax.jl index 404fa7e3d5..eaf40c1220 100644 --- a/validation/isopycnal_skew_symmetric_diffusion/new_front_relax.jl +++ b/validation/isopycnal_skew_symmetric_diffusion/new_front_relax.jl @@ -3,6 +3,7 @@ using Statistics using Random using Oceananigans using Oceananigans.Units +using Oceananigans.Grids: znode using GLMakie using Oceananigans.TurbulenceClosures: IsopycnalSkewSymmetricDiffusivity using Oceananigans.TurbulenceClosures: TriadIsopycnalSkewSymmetricDiffusivity @@ -28,12 +29,11 @@ stop_time = 1hours viscAh=300. viscAz=2.e-4 -grid = RectilinearGrid(topology = (Periodic, Bounded, Bounded), - size = (Nx, Ny, Nz), - x = (-Ly/2, Ly/2), +grid = RectilinearGrid(topology = (Flat, Bounded, Bounded), + size = (Ny, Nz), y = (-Ly/2, Ly/2), z = ([0 cumsum(reverse(dz))']'.-Lz), - halo = (3, 3, 3)) + halo = (3, 3)) #coriolis = FPlane(latitude = -45) coriolis = FPlane( f = 1.e-4) @@ -42,7 +42,8 @@ h_visc= HorizontalScalarDiffusivity( ν=viscAh ) v_visc= VerticalScalarDiffusivity( ν=viscAz ) gerdes_koberle_willebrand_tapering = FluxTapering(1e-2) -triad_closure = TriadIsopycnalSkewSymmetricDiffusivity(κ_skew = 0e3, + +triad_closure = TriadIsopycnalSkewSymmetricDiffusivity(VerticallyImplicitTimeDiscretization(), κ_skew = 0e3, κ_symmetric = 1e3, slope_limiter = gerdes_koberle_willebrand_tapering) @@ -53,33 +54,19 @@ cox_closure = IsopycnalSkewSymmetricDiffusivity(κ_skew = 0e3, @info "Building a model..." model = HydrostaticFreeSurfaceModel(; grid, coriolis, - closure = (triad_closure, h_visc, v_visc), + closure = triad_closure, #, h_visc, v_visc), buoyancy = BuoyancyTracer(), tracers = (:b, :c)) @info "Built $model." -""" -Linear ramp from 0 to 1 between -Δy/2 and +Δy/2. - -For example: - -y < y₀ => ramp = 0 -y₀ < y < y₀ + Δy => ramp = y / Δy -y > y₀ + Δy => ramp = 1 -""" -function ramp(x, y, Δ) - gradient == "x" && return min(max(0, x / Δ + 1/2), 1) - gradient == "y" && return min(max(0, y / Δ + 1/2), 1) -end - # 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=z[6] +yC = 0 +zC = znode(6, grid, Center()) Δy = Ly/8 Δz = 500 @@ -87,10 +74,8 @@ zC=z[6] Δb = Ly/Ny * M² ϵb = 1e-2 * Δb # noise amplitude -# bᵢ(x, y, z) = N² * z + Δb * ramp(x, y, Δy) -# cᵢ(x, y, z) = exp(-y^2 / 2Δc^2) * exp(-(z + Lz/4)^2 / 2Δz^2) -bᵢ(x, y, z) = N² * z - M² * Ly*sin(pi*y/Ly) *exp(-(3*z/Lz)^2) -cᵢ(x, y, z) = exp( -((y-yC)/Δy)^2 )*exp( -((z-zC)/Δz).^2); +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ᵢ) @@ -98,7 +83,7 @@ set!(model, b=bᵢ, c=cᵢ) ##### Simulation building ##### -simulation = Simulation(model; Δt, stop_time) +simulation = Simulation(model; Δt, stop_time, stop_iteration = 10) wall_clock = Ref(time_ns()) @@ -120,10 +105,10 @@ 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) +# 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..." From 11f952ffdd502a8ddfb43875fbe5bedd0e519cca Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Thu, 21 Nov 2024 16:42:49 -0800 Subject: [PATCH 14/15] make plots --- .../new_front_relax.jl | 30 +++++++------------ 1 file changed, 10 insertions(+), 20 deletions(-) diff --git a/validation/isopycnal_skew_symmetric_diffusion/new_front_relax.jl b/validation/isopycnal_skew_symmetric_diffusion/new_front_relax.jl index eaf40c1220..20d6aa6f50 100644 --- a/validation/isopycnal_skew_symmetric_diffusion/new_front_relax.jl +++ b/validation/isopycnal_skew_symmetric_diffusion/new_front_relax.jl @@ -20,10 +20,9 @@ Nx = 3 Ny = 32 Nz = 15 #save_fields_interval = 1day -#stop_time = 2days +stop_time = 30days save_fields_interval = 15minutes -stop_time = 1hours -Δt = 15minutes +Δt = 30minutes # viscosity viscAh=300. @@ -83,7 +82,7 @@ set!(model, b=bᵢ, c=cᵢ) ##### Simulation building ##### -simulation = Simulation(model; Δt, stop_time, stop_iteration = 10) +simulation = Simulation(model; Δt, stop_time) #, stop_iteration = 10) wall_clock = Ref(time_ns()) @@ -105,10 +104,10 @@ 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) +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..." @@ -141,18 +140,9 @@ z = z .* zscale times = bt.times Nt = length(times) -if gradient == "y" # average in x - # un(n) = interior(ut[n], 1, :, :) - # bn(n) = interior(bt[n], 1, :, :) - # cn(n) = interior(ct[n], 1, :, :) - 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, :, :) -else # average in y - un(n) = interior(mean(ut[n], dims=2), :, 1, :) - bn(n) = interior(mean(bt[n], dims=2), :, 1, :) - cn(n) = interior(mean(ct[n], dims=2), :, 1, :) -end +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 From 00325640d3fd19f062d3fdd05afa9086baaccbe1 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Thu, 21 Nov 2024 16:48:02 -0800 Subject: [PATCH 15/15] not really correct but runs --- ...pycnal_skew_symmetric_diffusivity_with_triads.jl | 13 +++++-------- .../new_front_relax.jl | 4 ++-- 2 files changed, 7 insertions(+), 10 deletions(-) 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 index bc8eab460e..0aacb00910 100644 --- 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 @@ -156,15 +156,12 @@ end @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 have at least one immersed interface -@inline triad_mask_x(ix, iz, j, kx, kz, grid) = one(grid) -@inline triad_mask_y(ix, iz, j, kx, kz, grid) = one(grid) +# 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_x(ix, iz, j, kx, kz, grid::ImmersedBoundaryGrid) = - immersed_peripheral_node(ix, j, kx, grid) | immersed_peripheral_node(iz, j, kz, grid) - -@inline triad_mask_y(i, jy, jz, ky, kz, grid::ImmersedBoundaryGrid) = - immersed_peripheral_node(i, jy, ky, grid) | immersed_peripheral_node(i, jz, kz, grid) +@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) diff --git a/validation/isopycnal_skew_symmetric_diffusion/new_front_relax.jl b/validation/isopycnal_skew_symmetric_diffusion/new_front_relax.jl index 20d6aa6f50..5fbae78de7 100644 --- a/validation/isopycnal_skew_symmetric_diffusion/new_front_relax.jl +++ b/validation/isopycnal_skew_symmetric_diffusion/new_front_relax.jl @@ -21,8 +21,8 @@ Ny = 32 Nz = 15 #save_fields_interval = 1day stop_time = 30days -save_fields_interval = 15minutes -Δt = 30minutes +save_fields_interval = 1days +Δt = 15minutes # viscosity viscAh=300.