Library

AdvectionTerm(model; velocities, location)

Return a KernelFunctionOperation that computes the advection term, defined as

ADV = uᵢ∂ⱼ(uᵢuⱼ)

By default, the buoyancy production will be calculated using the resolved velocities and users cab use the keyword velocities to modify that behavior:

julia> using Oceananigans

julia> grid = RectilinearGrid(size = (1, 1, 4), extent = (1,1,1));

julia> model = NonhydrostaticModel(grid=grid);

julia> using Oceanostics.TKEBudgetTerms: BuoyancyProductionTerm

julia> using Oceanostics.TKEBudgetTerms: AdvectionTerm

julia> ADV = AdvectionTerm(model)
KernelFunctionOperation at (Center, Center, Center)
├── grid: 1×1×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 1×1×3 halo
├── kernel_function: uᵢ∂ⱼuⱼuᵢᶜᶜᶜ (generic function with 1 method)
└── arguments: ("(u=1×1×4 Field{Face, Center, Center} on RectilinearGrid on CPU, v=1×1×4 Field{Center, Face, Center} on RectilinearGrid on CPU, w=1×1×5 Field{Center, Center, Face} on RectilinearGrid on CPU)", "Centered reconstruction order 2")

BuoyancyProductionTerm(model; velocities, tracers, location)

Return a KernelFunctionOperation that computes the buoyancy production term, defined as

BP = uᵢbᵢ

where bᵢ is the component of the buoyancy acceleration in the i-th direction (which is zero for x and y, except when gravity_unit_vector isn't aligned with the grid's z-direction) and all three components of i=1,2,3 are added up.

By default, the buoyancy production will be calculated using the resolved velocities and tracers:

julia> using Oceananigans

julia> grid = RectilinearGrid(size = (1, 1, 4), extent = (1,1,1));

julia> model = NonhydrostaticModel(grid=grid, buoyancy=BuoyancyTracer(), tracers=:b);

julia> using Oceanostics.TKEBudgetTerms: BuoyancyProductionTerm

julia> wb = BuoyancyProductionTerm(model)
KernelFunctionOperation at (Center, Center, Center)
├── grid: 1×1×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 1×1×3 halo
├── kernel_function: uᵢbᵢᶜᶜᶜ (generic function with 1 method)
└── arguments: ("(u=1×1×4 Field{Face, Center, Center} on RectilinearGrid on CPU, v=1×1×4 Field{Center, Face, Center} on RectilinearGrid on CPU, w=1×1×5 Field{Center, Center, Face} on RectilinearGrid on CPU)", "BuoyancyTracer with ĝ = NegativeZDirection()", "(b=1×1×4 Field{Center, Center, Center} on RectilinearGrid on CPU,)")

If we want to calculate only the turbulent buoyancy production rate, we can do so by passing turbulent perturbations to the velocities and/or tracers options):

julia> w′ = Field(model.velocities.w - Field(Average(model.velocities.w)));

julia> b′ = Field(model.tracers.b - Field(Average(model.tracers.b)));

julia> w′b′ = BuoyancyProductionTerm(model, velocities=(u=model.velocities.u, v=model.velocities.v, w=w′), tracers=(b=b′,))
KernelFunctionOperation at (Center, Center, Center)
├── grid: 1×1×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 1×1×3 halo
├── kernel_function: uᵢbᵢᶜᶜᶜ (generic function with 1 method)
└── arguments: ("(u=1×1×4 Field{Face, Center, Center} on RectilinearGrid on CPU, v=1×1×4 Field{Center, Face, Center} on RectilinearGrid on CPU, w=1×1×5 Field{Center, Center, Face} on RectilinearGrid on CPU)", "BuoyancyTracer with ĝ = NegativeZDirection()", "(b=1×1×4 Field{Center, Center, Center} on RectilinearGrid on CPU,)")

IsotropicKineticEnergyDissipationRate(
    model;
    U,
    V,
    W,
    location
)

Calculate the Viscous Dissipation Rate, defined as

ε = 2 ν SᵢⱼSᵢⱼ,

where Sᵢⱼ is the strain rate tensor, for a fluid with an isotropic turbulence closure (i.e., a turbulence closure where ν (eddy or not) is the same for all directions.

KineticEnergy(model, u, v, w; location, kwargs...)

Calculate the kinetic energy of model manually specifying u, v and w.

KineticEnergy(model; kwargs...)

Calculate the kinetic energy of model.

KineticEnergyDissipationRate(model; U, V, W, location)

Calculate the Kinetic Energy Dissipation Rate, defined as

ε = ν (∂uᵢ/∂xⱼ) (∂uᵢ/∂xⱼ)
ε = ∂ⱼuᵢ ⋅ Fᵢⱼ

where ∂ⱼuᵢ is the velocity gradient tensor and Fᵢⱼ is the stress tensor.

KineticEnergyForcingTerm(model; location)

Return a KernelFunctionOperation that computes the forcing term of the KE prognostic equation:

    FORC = uᵢFᵤᵢ

where uᵢ are the velocity components and Fᵤᵢ is the forcing term(s) in the uᵢ prognostic equation (i.e. the forcing for uᵢ).

KineticEnergyStressTerm(model; location)

Return a KernelFunctionOperation that computes the diffusive term of the KE prognostic equation:

    DIFF = uᵢ∂ⱼτᵢⱼ

where uᵢ are the velocity components and τᵢⱼ is the diffusive flux of i momentum in the j-th direction.

KineticEnergyTendency(model; location)

Return a KernelFunctionOperation that computes the tendency uᵢGᵢ of the KE, excluding the nonhydrostatic pressure contribution:

KET = ½∂ₜuᵢ² = uᵢGᵢ - uᵢ∂ᵢpₙₕₛ

julia> using Oceananigans

julia> grid = RectilinearGrid(size = (1, 1, 4), extent = (1,1,1));

julia> model = NonhydrostaticModel(grid=grid);

julia> using Oceanostics.TKEBudgetTerms: KineticEnergyTendency

julia> ke_tendency = KineticEnergyTendency(model)
KernelFunctionOperation at (Center, Center, Center)
├── grid: 1×1×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 1×1×3 halo
├── kernel_function: uᵢ∂ₜuᵢᶜᶜᶜ (generic function with 1 method)
└── arguments: ("Centered reconstruction order 2", "Nothing", "Nothing", "Nothing", "FluxBoundaryCondition: Nothing", "FluxBoundaryCondition: Nothing", "FluxBoundaryCondition: Nothing", "Nothing", "(velocities=(u=ZeroField{Int64}, v=ZeroField{Int64}, w=ZeroField{Int64}), tracers=NamedTuple())", "(u=1×1×4 Field{Face, Center, Center} on RectilinearGrid on CPU, v=1×1×4 Field{Center, Face, Center} on RectilinearGrid on CPU, w=1×1×5 Field{Center, Center, Face} on RectilinearGrid on CPU)", "NamedTuple()", "NamedTuple()", "Nothing", "(u=zeroforcing (generic function with 1 method), v=zeroforcing (generic function with 1 method), w=zeroforcing (generic function with 1 method))", "Nothing", "Clock(time=0 seconds, iteration=0, last_Δt=Inf days)")

PressureRedistributionTerm(
    model;
    velocities,
    pressure,
    location
)

Return a KernelFunctionOperation that computes the pressure redistribution term:

PR = uᵢ∂ᵢp

where p is the pressure. By default p is taken to be the total pressure (nonhydrostatic + hydrostatic):

julia> using Oceananigans

julia> grid = RectilinearGrid(size = (1, 1, 4), extent = (1,1,1));

julia> model = NonhydrostaticModel(grid=grid);

julia> using Oceanostics.TKEBudgetTerms: PressureRedistributionTerm

julia> ∇u⃗p = PressureRedistributionTerm(model)
KernelFunctionOperation at (Center, Center, Center)
├── grid: 1×1×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 1×1×3 halo
├── kernel_function: uᵢ∂ᵢpᶜᶜᶜ (generic function with 1 method)
└── arguments: ("(u=1×1×4 Field{Face, Center, Center} on RectilinearGrid on CPU, v=1×1×4 Field{Center, Face, Center} on RectilinearGrid on CPU, w=1×1×5 Field{Center, Center, Face} on RectilinearGrid on CPU)", "1×1×4 Field{Center, Center, Center} on RectilinearGrid on CPU")

We can also pass velocities and pressure keywords to perform more specific calculations. The example below illustrates calculation of the nonhydrostatic contribution to the pressure redistrubution term:

julia> ∇u⃗pNHS = PressureRedistributionTerm(model, pressure=model.pressures.pNHS)
KernelFunctionOperation at (Center, Center, Center)
├── grid: 1×1×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 1×1×3 halo
├── kernel_function: uᵢ∂ᵢpᶜᶜᶜ (generic function with 1 method)
└── arguments: ("(u=1×1×4 Field{Face, Center, Center} on RectilinearGrid on CPU, v=1×1×4 Field{Center, Face, Center} on RectilinearGrid on CPU, w=1×1×5 Field{Center, Center, Face} on RectilinearGrid on CPU)", "1×1×4 Field{Center, Center, Center} on RectilinearGrid on CPU")

TurbulentKineticEnergy(model, u, v, w; U, V, W, location)

Calculate the turbulent kinetic energy of model manually specifying u, v, w and optionally background velocities U, V and W.

TurbulentKineticEnergy(model; kwargs...)

Calculate the turbulent kinetic energy of model.

XShearProductionRate(model, u, v, w, U, V, W; location)

Calculate the shear production rate in the model's x direction, considering velocities u, v, w and background (or average) velocities U, V and W.

XShearProductionRate(model; U, V, W, kwargs...)

Calculate the shear production rate in the model's x direction. At least one of the mean velocities U, V and W must be specified otherwise the output will be zero.

YShearProductionRate(model, u, v, w, U, V, W; location)

Calculate the shear production rate in the model's y direction, considering velocities u, v, w and background (or average) velocities U, V and W.

YShearProductionRate(model; U, V, W, kwargs...)

Calculate the shear production rate in the model's y direction. At least one of the mean velocities U, V and W must be specified otherwise the output will be zero.

ZShearProductionRate(model, u, v, w, U, V, W; location)

Calculate the shear production rate in the model's z direction, considering velocities u, v, w and background (or average) velocities U, V and W.

ZShearProductionRate(model; U, V, W, kwargs...)

Calculate the shear production rate in the model's z direction. At least one of the mean velocities U, V and W must be specified otherwise the output will be zero.

TracerVarianceDiffusiveTerm(model, tracer_name; location)

Return a KernelFunctionOperation that computes the diffusive term of the tracer variance prognostic equation using Oceananigans' diffusive tracer flux divergence kernel:

    DIFF = 2 c ∂ⱼFⱼ

where c is the tracer, and Fⱼ is the tracer's diffusive flux in the j-th direction.

grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1))
model = NonhydrostaticModel(grid=grid, tracers=:b, closure=SmagorinskyLilly())

DIFF = TracerVarianceDiffusiveTerm(model, :b)

TracerVarianceDissipationRate(
    model,
    tracer_name;
    tracer,
    location
)

Return a KernelFunctionOperation that computes the isotropic variance dissipation rate for tracer_name in model.tracers. The isotropic variance dissipation rate is defined as

    χ = 2 ∂ⱼc ⋅ Fⱼ

where Fⱼ is the diffusive flux of c in the j-th direction and ∂ⱼ is the gradient operator. χ is implemented in its conservative formulation based on the equation above.

Note that often χ is written as χ = 2κ (∇c ⋅ ∇c), which is the special case for Fickian diffusion (κ is the tracer diffusivity).

Here tracer_name is needed even when passing tracer in order to get the appropriate tracer_index. When passing tracer, this function should be used as

grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1))
model = NonhydrostaticModel(grid=grid, tracers=:b, closure=SmagorinskyLilly())

b̄ = Field(Average(model.tracers.b, dims=(1,2)))
b′ = model.tracers.b - b̄

χb = TracerVarianceDissipationRate(model, :b, tracer=b′)

TracerVarianceTendency(model, tracer_name; location)

Return a KernelFunctionOperation that computes the tracer variance tendency:

TEND = 2 c ∂ₜc

where c is the tracer and ∂ₜc is the tracer tendency (computed using Oceananigans' tracer tendency kernel).

grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1))
model = NonhydrostaticModel(grid=grid, tracers=:b, closure=SmagorinskyLilly())

DIFF = TracerVarianceTendency(model, :b)

DirectionalErtelPotentialVorticity(
    model,
    direction;
    location
)

Calculate the contribution from a given direction to the Ertel Potential Vorticity basde on a model and a direction. The Ertel Potential Vorticity is defined as

EPV = ωₜₒₜ ⋅ ∇b

where ωₜₒₜ is the total (relative + planetary) vorticity vector, b is the buoyancy and ∇ is the gradient operator.

ErtelPotentialVorticity(model; location, add_background)

Calculate the Ertel Potential Vorticty for model, where the characteristics of the Coriolis rotation are taken from model.coriolis. The Ertel Potential Vorticity is defined as

EPV = ωₜₒₜ ⋅ ∇b

where ωₜₒₜ is the total (relative + planetary) vorticity vector, b is the buoyancy and ∇ is the gradient operator.

QVelocityGradientTensorInvariant(model; location)

Calculate the value of the Q velocity gradient tensor invariant. This is usually just called Q and it is generally used for identifying and visualizing vortices in fluid flow.

The definition and nomenclature comes from the equation for the eigenvalues λ of the velocity gradient tensor ∂ⱼuᵢ:

    λ³ + P λ² + Q λ + T = 0

from where Q is defined as

    Q = ½ ( ΩᵢⱼΩᵢⱼ - SᵢⱼSᵢⱼ)

and where Sᵢⱼ= ½(∂ⱼuᵢ + ∂ᵢuⱼ) and Ωᵢⱼ= ½(∂ⱼuᵢ - ∂ᵢuⱼ). More info about it can be found in doi:10.1063/1.5124245.

RichardsonNumber(model; location, add_background)

Calculate the Richardson Number as

    Ri = (∂b/∂z) / (|∂u⃗ₕ/∂z|²)

where z is the true vertical direction (ie anti-parallel to gravity).

RossbyNumber(
    model;
    location,
    add_background,
    dWdy_bg,
    dVdz_bg,
    dUdz_bg,
    dWdx_bg,
    dUdy_bg,
    dVdx_bg
)

Calculate the Rossby number using the vorticity in the rotation axis direction according to model.coriolis. Rossby number is defined as

    Ro = ωᶻ / f

where ωᶻ is the vorticity in the Coriolis axis of rotation and f is the Coriolis rotation frequency.

StrainRateTensorModulus(model; location)

Calculate the modulus (absolute value) of the strain rate tensor S, which is defined as the symmetric part of the velocity gradient tensor:

    Sᵢⱼ = ½(∂ⱼuᵢ + ∂ᵢuⱼ)

Its modulus is then defined (using Einstein summation notation) as

    || Sᵢⱼ || = √( Sᵢⱼ Sᵢⱼ)

ThermalWindPotentialVorticity(model; f, location)

Calculate the Potential Vorticty assuming thermal wind balance for model, where the characteristics of the Coriolis rotation are taken from model.coriolis. The Potential Vorticity in this case is defined as

    TWPV = (f + ωᶻ) ∂b/∂z - f ((∂U/∂z)² + (∂V/∂z)²)

where f is the Coriolis frequency, ωᶻ is the relative vorticity in the z direction, b is the buoyancy, and ∂U/∂z and ∂V/∂z comprise the thermal wind shear.

VorticityTensorModulus(model; location)

Calculate the modulus (absolute value) of the vorticity tensor Ω, which is defined as the antisymmetric part of the velocity gradient tensor:

    Ωᵢⱼ = ½(∂ⱼuᵢ - ∂ᵢuⱼ)

Its modulus is then defined (using Einstein summation notation) as

    || Ωᵢⱼ || = √( Ωᵢⱼ Ωᵢⱼ)

PotentialEnergy(model; location, geopotential_height)

Return a KernelFunctionOperation to compute the PotentialEnergy per unit volume,

\[Eₚ = \frac{gρz}{ρ₀}\]

at each grid location in model. PotentialEnergy is defined for both BuoyancyTracer and SeawaterBuoyancy. See the relevant Oceananigans.jl documentation on buoyancy models for more information about available options.

The optional keyword argument geopotential_height is only used if ones wishes to calculate Eₚ with a potential density referenced to geopotential_height, rather than in-situ density, when using a BoussinesqEquationOfState.

Example

Usage with a BuoyancyTracer buoyacny model

julia> using Oceananigans

julia> using Oceanostics.PotentialEnergyEquationTerms: PotentialEnergy

julia> grid = RectilinearGrid(size=100, z=(-1000, 0), topology=(Flat, Flat, Bounded))
1×1×100 RectilinearGrid{Float64, Flat, Flat, Bounded} on CPU with 0×0×3 halo
├── Flat x
├── Flat y
└── Bounded  z ∈ [-1000.0, 0.0]   regularly spaced with Δz=10.0

julia> model = NonhydrostaticModel(; grid, buoyancy=BuoyancyTracer(), tracers=(:b,))
NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 1×1×100 RectilinearGrid{Float64, Flat, Flat, Bounded} on CPU with 0×0×3 halo
├── timestepper: QuasiAdamsBashforth2TimeStepper
├── advection scheme: Centered reconstruction order 2
├── tracers: b
├── closure: Nothing
├── buoyancy: BuoyancyTracer with ĝ = NegativeZDirection()
└── coriolis: Nothing

julia> PotentialEnergy(model)
KernelFunctionOperation at (Center, Center, Center)
├── grid: 1×1×100 RectilinearGrid{Float64, Flat, Flat, Bounded} on CPU with 0×0×3 halo
├── kernel_function: bz_ccc (generic function with 2 methods)
└── arguments: ("1×1×100 Field{Center, Center, Center} on RectilinearGrid on CPU",)

The default behaviour of PotentialEnergy uses the in-situ density in the calculation when the equation of state is a BoussinesqEquationOfState:

julia> using Oceananigans, SeawaterPolynomials.TEOS10

julia> using Oceanostics.PotentialEnergyEquationTerms: PotentialEnergy

julia> grid = RectilinearGrid(size=100, z=(-1000, 0), topology=(Flat, Flat, Bounded))
1×1×100 RectilinearGrid{Float64, Flat, Flat, Bounded} on CPU with 0×0×3 halo
├── Flat x
├── Flat y
└── Bounded  z ∈ [-1000.0, 0.0]   regularly spaced with Δz=10.0

julia> tracers = (:T, :S)
(:T, :S)

julia> eos = TEOS10EquationOfState()
BoussinesqEquationOfState{Float64}:
    ├── seawater_polynomial: TEOS10SeawaterPolynomial{Float64}
    └── reference_density: 1020.0

julia> buoyancy = SeawaterBuoyancy(equation_of_state=eos)
SeawaterBuoyancy{Float64}:
├── gravitational_acceleration: 9.80665
└── equation_of_state: BoussinesqEquationOfState{Float64}

julia> model = NonhydrostaticModel(; grid, buoyancy, tracers)
NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 1×1×100 RectilinearGrid{Float64, Flat, Flat, Bounded} on CPU with 0×0×3 halo
├── timestepper: QuasiAdamsBashforth2TimeStepper
├── advection scheme: Centered reconstruction order 2
├── tracers: (T, S)
├── closure: Nothing
├── buoyancy: SeawaterBuoyancy with g=9.80665 and BoussinesqEquationOfState{Float64} with ĝ = NegativeZDirection()
└── coriolis: Nothing

julia> PotentialEnergy(model)
KernelFunctionOperation at (Center, Center, Center)
├── grid: 1×1×100 RectilinearGrid{Float64, Flat, Flat, Bounded} on CPU with 0×0×3 halo
├── kernel_function: g′z_ccc (generic function with 1 method)
└── arguments: ("KernelFunctionOperation at (Center, Center, Center)", "(g=9.80665, ρ₀=1020.0)")

To use a reference density set a constant value for the keyword argument geopotential_height and pass this the function. For example,

julia> using Oceananigans, SeawaterPolynomials.TEOS10;

julia> using Oceanostics.PotentialEnergyEquationTerms: PotentialEnergy;

julia> grid = RectilinearGrid(size=100, z=(-1000, 0), topology=(Flat, Flat, Bounded));

julia> tracers = (:T, :S);

julia> eos = TEOS10EquationOfState();

julia> buoyancy = SeawaterBuoyancy(equation_of_state=eos);

julia> model = NonhydrostaticModel(; grid, buoyancy, tracers);

julia> geopotential_height = 0; # density variable will be σ₀

julia> PotentialEnergy(model)
KernelFunctionOperation at (Center, Center, Center)
├── grid: 1×1×100 RectilinearGrid{Float64, Flat, Flat, Bounded} on CPU with 0×0×3 halo
├── kernel_function: g′z_ccc (generic function with 1 method)
└── arguments: ("KernelFunctionOperation at (Center, Center, Center)", "(g=9.80665, ρ₀=1020.0)")