Skip to content

Commit

Permalink
Fix RiBasedVerticalDiffusivity (#3510)
Browse files Browse the repository at this point in the history 
* Fix max/min bug with RiBasedVerticalDiffusivity

* Add optional horizontal filtering

* Add a minimum entrainment buoyancy gradient

* Fix Ri calculation

* add missing comma

* Update src/TurbulenceClosures/turbulence_closure_implementations/ri_based_vertical_diffusivity.jl

Co-authored-by: Navid C. Constantinou <navidcy@users.noreply.github.com>

* More bugfix

* Update CATKEVerticalDiffusivities.jl

---------

Co-authored-by: Navid C. Constantinou <navidcy@users.noreply.github.com>
  • Loading branch information
@glwagner @navidcy
glwagner and navidcy committed Apr 10, 2024
1 parent b3d7a8a commit 00f028b
Show file tree
Hide file tree
Showing 2 changed files with 60 additions and 24 deletions.
2 changes: 1 addition & 1 deletion ...rbulence_closure_implementations/CATKEVerticalDiffusivities/CATKEVerticalDiffusivities.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -145,7 +145,7 @@ function CATKEVerticalDiffusivity(time_discretization::TD = VerticallyImplicitTi
maximum_tke_diffusivity = Inf,
maximum_viscosity = Inf,
minimum_turbulent_kinetic_energy = 1e-6,
minimum_convective_buoyancy_flux = 1e-8,
minimum_convective_buoyancy_flux = 1e-11,
negative_turbulent_kinetic_energy_damping_time_scale = 1minute) where TD

mixing_length = convert_eltype(FT, mixing_length)
Expand Down
82 changes: 59 additions & 23 deletions src/TurbulenceClosures/turbulence_closure_implementations/ri_based_vertical_diffusivity.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@ using Oceananigans.Operators
using Oceananigans.Grids: inactive_node
using Oceananigans.Operators: ℑzᵃᵃᶜ

struct RiBasedVerticalDiffusivity{TD, FT, R} <: AbstractScalarDiffusivity{TD, VerticalFormulation, 1}
struct RiBasedVerticalDiffusivity{TD, FT, R, HR} <: AbstractScalarDiffusivity{TD, VerticalFormulation, 1}
ν₀ :: FT
κ₀ :: FT
κᶜᵃ :: FT
Expand All @@ -13,6 +13,8 @@ struct RiBasedVerticalDiffusivity{TD, FT, R} <: AbstractScalarDiffusivity{TD, Ve
Ri₀ :: FT
Riᵟ :: FT
Ri_dependent_tapering :: R
horizontal_Ri_filter :: HR
minimum_entrainment_buoyancy_gradient :: FT
maximum_diffusivity :: FT
maximum_viscosity :: FT
end
Expand All @@ -25,14 +27,18 @@ function RiBasedVerticalDiffusivity{TD}(ν₀::FT,
Ri₀::FT,
Riᵟ::FT,
Ri_dependent_tapering::R,
horizontal_Ri_filter::HR,
minimum_entrainment_buoyancy_gradient::FT,
maximum_diffusivity::FT,
maximum_viscosity::FT) where {TD, FT, R}
maximum_viscosity::FT) where {TD, FT, R, HR}


return RiBasedVerticalDiffusivity{TD, FT, R}(ν₀, κ₀, κᶜᵃ, Cᵉⁿ, Cᵃᵛ, Ri₀, Riᵟ,
Ri_dependent_tapering,
maximum_diffusivity,
maximum_viscosity)
return RiBasedVerticalDiffusivity{TD, FT, R, HR}(ν₀, κ₀, κᶜᵃ, Cᵉⁿ, Cᵃᵛ, Ri₀, Riᵟ,
Ri_dependent_tapering,
horizontal_Ri_filter,
minimum_entrainment_buoyancy_gradient,
maximum_diffusivity,
maximum_viscosity)
end

# Ri-dependent tapering flavor
Expand All @@ -44,10 +50,19 @@ Base.summary(::HyperbolicTangentRiDependentTapering) = "HyperbolicTangentRiDepen
Base.summary(::ExponentialRiDependentTapering) = "ExponentialRiDependentTapering"
Base.summary(::PiecewiseLinearRiDependentTapering) = "PiecewiseLinearRiDependentTapering"

# Horizontal filtering for the Richardson number
struct FivePointHorizontalFilter end
@inline filter_horizontally(i, j, k, grid, ::Nothing, ϕ) = @inbounds ϕ[i, j, k]
@inline filter_horizontally(i, j, k, grid, ::FivePointHorizontalFilter, ϕ) = ℑxyᶜᶜᵃ(i, j, k, grid, ℑxyᶠᶠᵃ, ϕ)

"""
RiBasedVerticalDiffusivity([time_discretization = VerticallyImplicitTimeDiscretization(),
FT = Float64;]
Ri_dependent_tapering = HyperbolicTangentRiDependentTapering(),
horizontal_Ri_filter = nothing,
minimum_entrainment_buoyancy_gradient = 1e-10,
maximum_diffusivity = Inf,
maximum_viscosity = Inf,
ν₀ = 0.7,
κ₀ = 0.5,
κᶜᵃ = 1.7,
Expand Down Expand Up @@ -79,28 +94,39 @@ Keyword arguments
`HyperbolicTangentRiDependentTapering()` (default), and
`ExponentialRiDependentTapering()`.
* `ν₀`: Non-convective viscosity.
* `ν₀`: Non-convective viscosity (units of kinematic viscosity, typically m² s⁻¹).
* `κ₀`: Non-convective diffusivity for tracers (units of diffusivity, typically m² s⁻¹).
* `κ₀`: Non-convective diffusivity for tracers.
* `κᶜᵃ`: Convective adjustment diffusivity for tracers (units of diffusivity, typically m² s⁻¹).
* `κᶜᵃ`: Convective adjustment diffusivity for tracers.
* `Cᵉⁿ`: Entrainment coefficient for tracers (non-dimensional).
Set `Cᵉⁿ = 0` to turn off the penetrative entrainment diffusivity.
* `Cᵉⁿ`: Entrainment coefficient for tracers.
* `Cᵃᵛ`: Time-averaging coefficient for viscosity and diffusivity (non-dimensional).
* `Cᵃᵛ`: Time-averaging coefficient for viscosity and diffusivity.
* `Ri₀`: ``Ri`` threshold for decreasing viscosity and diffusivity (non-dimensional).
* `Ri₀`: ``Ri`` threshold for decreasing viscosity and diffusivity.
* `Riᵟ`: ``Ri``-width over which viscosity and diffusivity decreases to 0 (non-dimensional).
* `Riᵟ`: ``Ri``-width over which viscosity and diffusivity decreases to 0.
* `minimum_entrainment_buoyancy_gradient`: Minimum buoyancy gradient for application of the entrainment
diffusvity. If the entrainment buoyancy gradient is less than the
minimum value, the entrainment diffusivity is 0. Units of
buoyancy gradient (typically s⁻²).
* `maximum_diffusivity`: A cap on the diffusivity
* `maximum_diffusivity`: A limiting maximum tracer diffusivity (units of diffusivity, typically m² s⁻¹).
* `maximum_viscosity`: A cap on viscosity
* `maximum_viscosity`: A limiting maximum viscosity (units of kinematic viscosity, typically m² s⁻¹).
* `horizontal_Ri_filter`: Horizontal filter to apply to Ri, which can help alleviate noise for
some simulations. The default is `nothing`, or no filtering. The other
option is `horizontal_Ri_filter = FivePointHorizontalFilter()`.
"""
function RiBasedVerticalDiffusivity(time_discretization = VerticallyImplicitTimeDiscretization(),
FT = Float64;
Ri_dependent_tapering = HyperbolicTangentRiDependentTapering(),
horizontal_Ri_filter = nothing,
minimum_entrainment_buoyancy_gradient = 1e-10,
maximum_diffusivity = Inf,
maximum_viscosity = Inf,
ν₀ = 0.7,
Expand Down Expand Up @@ -129,6 +155,8 @@ function RiBasedVerticalDiffusivity(time_discretization = VerticallyImplicitTime
convert(FT, Ri₀),
convert(FT, Riᵟ),
Ri_dependent_tapering,
horizontal_Ri_filter,
convert(FT, minimum_entrainment_buoyancy_gradient),
convert(FT, maximum_diffusivity),
convert(FT, maximum_viscosity))
end
Expand Down Expand Up @@ -182,6 +210,10 @@ function compute_diffusivities!(diffusivities, closure::FlavorOfRBVD, model; par
top_tracer_bcs,
clock)

# Use `only_local_halos` to ensure that no communication occurs during
# this call to fill_halo_regions!
fill_halo_regions!(diffusivities.Ri; only_local_halos=true)

launch!(arch, grid, parameters,
compute_ri_based_diffusivities!,
diffusivities,
Expand Down Expand Up @@ -217,9 +249,7 @@ const Tanh = HyperbolicTangentRiDependentTapering
end

@inline function Riᶜᶜᶠ(i, j, k, grid, velocities, buoyancy, tracers)
∂z_u² = ℑxᶜᵃᵃ(i, j, k, grid, ϕ², ∂zᶠᶜᶠ, velocities.u)
∂z_v² = ℑyᵃᶜᵃ(i, j, k, grid, ϕ², ∂zᶜᶠᶠ, velocities.v)
S² = ∂z_u² + ∂z_v²
S² = shear_squaredᶜᶜᶠ(i, j, k, grid, velocities)
N² = ∂z_b(i, j, k, grid, buoyancy, tracers)
Ri = N² / S²

Expand All @@ -243,6 +273,7 @@ end
velocities, tracers, buoyancy, tracer_bcs, clock)
end


@inline function _compute_ri_based_diffusivities!(i, j, k, diffusivities, grid, closure,
velocities, tracers, buoyancy, tracer_bcs, clock)

Expand All @@ -257,6 +288,8 @@ end
Ri₀ = closure_ij.Ri₀
Riᵟ = closure_ij.Riᵟ
tapering = closure_ij.Ri_dependent_tapering
Ri_filter = closure_ij.horizontal_Ri_filter
N²ᵉⁿ = closure_ij.minimum_entrainment_buoyancy_gradient
Qᵇ = top_buoyancy_flux(i, j, grid, buoyancy, tracer_bcs, clock, merge(velocities, tracers))

# Convection and entrainment
Expand All @@ -266,16 +299,18 @@ end
# Conditions
# TODO: apply a minimum entrainment buoyancy gradient?
convecting = N² < 0 # applies regardless of Qᵇ
entraining = (N² > 0) & (N²_above < 0) & (Qᵇ > 0)
entraining = (N² > N²ᵉⁿ) & (N²_above < 0) & (Qᵇ > 0)

# Convective adjustment diffusivity
κᶜᵃ = ifelse(convecting, κᶜᵃ, zero(grid))

# Entrainment diffusivity
κᵉⁿ = ifelse(entraining, Cᵉⁿ * Qᵇ / N², zero(grid))

# (Potentially) apply a horizontal filter to the Richardson number
Ri = filter_horizontally(i, j, k, grid, Ri_filter, diffusivities.Ri)

# Shear mixing diffusivity and viscosity
Ri = ℑxyᶜᶜᵃ(i, j, k, grid, ℑxyᶠᶠᵃ, diffusivities.Ri)
τ = taper(tapering, Ri, Ri₀, Riᵟ)
κᶜ★ = κ₀ * τ
κᵘ★ = ν₀ * τ
Expand All @@ -289,8 +324,8 @@ end
κᵘ⁺ = κᵘ★

# Limit by specified maximum
κᶜ⁺ = max(κᶜ⁺, closure_ij.maximum_diffusivity)
κᵘ⁺ = max(κᵘ⁺, closure_ij.maximum_viscosity)
κᶜ⁺ = min(κᶜ⁺, closure_ij.maximum_diffusivity)
κᵘ⁺ = min(κᵘ⁺, closure_ij.maximum_viscosity)

# Set to zero on periphery and NaN within inactive region
on_periphery = peripheral_node(i, j, k, grid, c, c, f)
Expand All @@ -315,11 +350,12 @@ function Base.show(io::IO, closure::RiBasedVerticalDiffusivity)
print(io, summary(closure), '\n')
print(io, "├── Ri_dependent_tapering: ", prettysummary(closure.Ri_dependent_tapering), '\n')
print(io, "├── κ₀: ", prettysummary(closure.κ₀), '\n')
print(io, "├── ν₀: ", prettysummary(closure.ν₀), '\n')
print(io, "├── κᶜᵃ: ", prettysummary(closure.κᶜᵃ), '\n')
print(io, "├── Cᵉⁿ: ", prettysummary(closure.Cᵉⁿ), '\n')
print(io, "├── Cᵃᵛ: ", prettysummary(closure.Cᵃᵛ), '\n')
print(io, "├── Ri₀: ", prettysummary(closure.Ri₀), '\n')
print(io, "├── Riᵟ: ", prettysummary(closure.Riᵟ), '\n')
print(io, "├── maximum_diffusivity: ", prettysummary(closure.maximum_diffusivity), '\n')
print(io, "└── maximum_viscosity: ", prettysummary(closure.maximum_viscosity))
end

0 comments on commit 00f028b

Please sign in to comment.