Flow diagnostics

RichardsonNumber(model; loc)

Calculate the Richardson Number as

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

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

julia> using Oceananigans, Oceanostics

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

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

julia> Ri = RichardsonNumber(model)
RichardsonNumber KernelFunctionOperation at (Center, Center, Face)
├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── kernel_function: richardson_number_ccf (generic function with 1 method)
└── arguments: ("Field", "Field", "Field", "Field", "Tuple")
└── computes: Richardson number  Ri = (∂b/∂z) / |∂u⃗ₕ/∂z|²

RossbyNumber(
    model;
    loc,
    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.

julia> using Oceananigans, Oceanostics

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

julia> model = NonhydrostaticModel(grid; coriolis=FPlane(f=1e-4));

julia> Ro = RossbyNumber(model)
RossbyNumber KernelFunctionOperation at (Face, Face, Face)
├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── kernel_function: rossby_number_fff (generic function with 1 method)
└── arguments: ("Field", "Field", "Field", "NamedTuple")
└── computes: Rossby number  Ro = ωᶻ/f

ErtelPotentialVorticity(
    model;
    tracer_name,
    thermal_wind,
    loc
)

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.

If thermal_wind = true, the thermal-wind approximation is used instead, giving

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

where f is the (vertical component of the) Coriolis frequency, ωᶻ is the vertical relative vorticity, and ∂U/∂z, ∂V/∂z comprise the thermal wind shear. The returned object is an instance of both ErtelPotentialVorticity and ThermalWindPotentialVorticity, so the thermal-wind variant can be identified separately.

julia> using Oceananigans

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

julia> N² = 1e-6;

julia> b_bcs = FieldBoundaryConditions(top=GradientBoundaryCondition(N²));

julia> model = NonhydrostaticModel(grid; coriolis=FPlane(f=1e-4), buoyancy=BuoyancyTracer(), tracers=:b, boundary_conditions=(; b=b_bcs));

julia> stratification(z) = N² * z;

julia> set!(model, b=stratification)

julia> using Oceanostics: ErtelPotentialVorticity

julia> EPV = ErtelPotentialVorticity(model)
ErtelPotentialVorticity KernelFunctionOperation at (Face, Face, Face)
├── grid: 1×1×4 RectilinearGrid{Float64, Flat, Flat, Bounded} on CPU with 0×0×3 halo
├── kernel_function: ertel_potential_vorticity_fff (generic function with 1 method)
└── arguments: ("Field", "Field", "Field", "Field", "Int64", "Int64", "Float64")
└── computes: Ertel potential vorticity  q = ω⃗ₜₒₜ · ∇b

julia> interior(Field(EPV))
1×1×5 view(::Array{Float64, 3}, 1:1, 1:1, 4:8) with eltype Float64:
[:, :, 1] =
 0.0

[:, :, 2] =
 1.0000000000000002e-10

[:, :, 3] =
 9.999999999999998e-11

[:, :, 4] =
 1.0000000000000002e-10

[:, :, 5] =
 1.0e-10

Note that EPV values are correctly calculated both in the interior and the boundaries. In the interior and top boundary, EPV = f×N² = 10⁻¹⁰, while EPV = 0 at the bottom boundary since ∂b/∂z is zero there.

ThermalWindPotentialVorticity

Narrower type alias matching only the thermal-wind variant of ErtelPotentialVorticity. Useful for identifying or dispatching on the thermal-wind variant via isa. Construct via ErtelPotentialVorticity(model; thermal_wind = true).

julia> using Oceananigans, Oceanostics

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

julia> model = NonhydrostaticModel(grid; coriolis=FPlane(f=1e-4), buoyancy=BuoyancyTracer(), tracers=:b);

julia> PV = ErtelPotentialVorticity(model; thermal_wind = true);

julia> PV isa ThermalWindPotentialVorticity
true

julia> PV isa ErtelPotentialVorticity
true

DirectionalErtelPotentialVorticity(
    model,
    direction;
    tracer_name,
    loc
)

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.

julia> using Oceananigans, Oceanostics

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

julia> model = NonhydrostaticModel(grid; coriolis=FPlane(f=1e-4), buoyancy=BuoyancyTracer(), tracers=:b);

julia> DEPV = DirectionalErtelPotentialVorticity(model, (0, 0, 1))
DirectionalErtelPotentialVorticity KernelFunctionOperation at (Face, Face, Face)
├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── kernel_function: directional_ertel_potential_vorticity_fff (generic function with 1 method)
└── arguments: ("Field", "Field", "Field", "Field", "NamedTuple")
└── computes: directional contribution to Ertel PV  (f̂ + ω̂)·∇b along a direction

StrainRateTensor(model; dims)

Return the components of the strain rate tensor S, defined as the symmetric part of the velocity gradient tensor:

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

The result is a NamedTuple with the independent components, each a KernelFunctionOperation living at its natural location on the staggered grid:

The tensor is symmetric, so the remaining components follow from Sⱼᵢ = Sᵢⱼ (i.e. S₂₁ = S₁₂, S₃₁ = S₁₃, S₃₂ = S₂₃).

dims selects which spatial directions (1 → x, 2 → y, 3 → z) enter the tensor: component Sᵢⱼ is included only when both i and j are in dims. The default dims = (1, 2, 3) returns the full tensor, while e.g. dims = (1, 3) returns the 2D strain rate tensor in the x–z plane (S₁₁, S₃₃, S₁₃). Components are always ordered diagonals-first, independently of the order of dims.

Each component can be wrapped in a Field and used with output writers, time-averaging, etc. Can also be called as StrainRateTensor(grid, u, v, w; dims) to build the components from individual velocity fields. See also StrainRateTensorModulus for the scalar modulus √(SᵢⱼSᵢⱼ).

julia> using Oceananigans, Oceanostics

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

julia> model = NonhydrostaticModel(grid);

julia> S = StrainRateTensor(model);

julia> keys(S)
(:S₁₁, :S₂₂, :S₃₃, :S₁₂, :S₁₃, :S₂₃)

julia> S.S₁₃
StrainRateTensor KernelFunctionOperation at (Face, Center, Face)
├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── kernel_function: Oceanostics.FlowDiagnostics.StrainRateTensorKernel{1, 3}
└── arguments: ("Field", "Field")
└── computes: strain-rate tensor component  Sᵢⱼ = ½(∂ⱼuᵢ + ∂ᵢuⱼ)

StrainRateTensorModulus(model; loc)

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ᵢⱼ)

julia> using Oceananigans, Oceanostics

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

julia> model = NonhydrostaticModel(grid);

julia> S = StrainRateTensorModulus(model)
StrainRateTensorModulus KernelFunctionOperation at (Center, Center, Center)
├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── kernel_function: strain_rate_tensor_modulus_ccc (generic function with 1 method)
└── arguments: ("Field", "Field", "Field")
└── computes: strain-rate tensor modulus  √(SᵢⱼSᵢⱼ)

VorticityTensor(model; dims)

Return the components of the vorticity tensor Ω, defined as the antisymmetric part of the velocity gradient tensor:

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

The tensor is antisymmetric, so its diagonal vanishes (Ω₁₁ = Ω₂₂ = Ω₃₃ = 0) and only the independent off-diagonal components are returned, as a NamedTuple of KernelFunctionOperations, each living at its natural location on the staggered grid:

The remaining off-diagonal components follow from antisymmetry, Ωⱼᵢ = -Ωᵢⱼ (i.e. Ω₂₁ = -Ω₁₂, Ω₃₁ = -Ω₁₃, Ω₃₂ = -Ω₂₃).

dims selects which spatial directions (1 → x, 2 → y, 3 → z) enter the tensor: component Ωᵢⱼ is included only when both i and j are in dims. The default dims = (1, 2, 3) returns all three off-diagonal components, while e.g. dims = (1, 3) returns the single component in the x–z plane (Ω₁₃). Because every component couples two distinct directions, a single-direction dims (e.g. dims = (1,)) yields an empty tensor. Components are always ordered Ω₁₂, Ω₁₃, Ω₂₃, independently of the order of dims.

Each component can be wrapped in a Field and used with output writers, time-averaging, etc. Can also be called as VorticityTensor(grid, u, v, w; dims) to build the components from individual velocity fields. See also VorticityTensorModulus for the scalar modulus √(ΩᵢⱼΩᵢⱼ).

julia> using Oceananigans, Oceanostics

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

julia> model = NonhydrostaticModel(grid);

julia> Ω = VorticityTensor(model);

julia> keys(Ω)
(:Ω₁₂, :Ω₁₃, :Ω₂₃)

julia> Ω.Ω₁₃
VorticityTensor KernelFunctionOperation at (Face, Center, Face)
├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── kernel_function: Oceanostics.FlowDiagnostics.VorticityTensorKernel{1, 3}
└── arguments: ("Field", "Field")
└── computes: vorticity tensor component  Ωᵢⱼ = ½(∂ⱼuᵢ - ∂ᵢuⱼ)

VorticityTensorModulus(model; loc)

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

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

julia> using Oceananigans, Oceanostics

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

julia> model = NonhydrostaticModel(grid);

julia> Ω = VorticityTensorModulus(model)
VorticityTensorModulus KernelFunctionOperation at (Center, Center, Center)
├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── kernel_function: vorticity_tensor_modulus_ccc (generic function with 1 method)
└── arguments: ("Field", "Field", "Field")
└── computes: vorticity tensor modulus  √(ΩᵢⱼΩᵢⱼ)

QVelocityGradientTensorInvariant(model; loc)

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.

julia> using Oceananigans, Oceanostics

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

julia> model = NonhydrostaticModel(grid);

julia> Qcrit = QVelocityGradientTensorInvariant(model)
QVelocityGradientTensorInvariant KernelFunctionOperation at (Center, Center, Center)
├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── kernel_function: Q_velocity_gradient_tensor_invariant_ccc (generic function with 1 method)
└── arguments: ("Field", "Field", "Field")
└── computes: Q velocity-gradient invariant  Q = ½(ΩᵢⱼΩᵢⱼ - SᵢⱼSᵢⱼ)

StressTensor(model; dims, collocate_diagonals)

Return the components of the (kinematic) stress tensor τ, defined as the outer product of the velocity field with itself:

    τᵢⱼ = uᵢ uⱼ

The result is a NamedTuple of the independent components, each a KernelFunctionOperation living at a location on the staggered grid.

The off-diagonal components are always evaluated at their natural edge location. Because each couples two different, mutually-staggered velocities, the two factors must be interpolated to a common point before multiplying — this is unavoidable, and the edge locations below are the interpolation-minimal choice (one interpolation per factor):

The collocate_diagonals keyword controls only where the diagonal components τ₁₁, τ₂₂, τ₃₃ live. It trades interpolation against collocation, and defaults to false (minimal interpolation):

• collocate_diagonals = false (default): each diagonal is computed as τᵢᵢ = uᵢ² read directly at uᵢ's own location, performing no interpolation at all. This is the cheapest and most accurate option. The trade-off is that the three diagonals end up at three different locations, so they are not collocated with each other or with the off-diagonals:

  Component	Definition	Location
  τ₁₁	u u	fcc
  τ₂₂	v v	cfc
  τ₃₃	w w	ccf

• collocate_diagonals = true: each velocity is first interpolated to ccc and then squared, τᵢᵢ = (ℑ uᵢ)², so all three diagonals share the single location ccc. Choose this when you need the diagonals collocated — e.g. to form the trace τ₁₁ + τ₂₂ + τ₃₃ (twice the kinetic energy) or to treat τ as a single collocated tensor. The cost is one interpolation per diagonal:

  Component	Definition	Location
  τ₁₁	u u	ccc
  τ₂₂	v v	ccc
  τ₃₃	w w	ccc

The tensor is symmetric, so the remaining components follow from τⱼᵢ = τᵢⱼ (i.e. τ₂₁ = τ₁₂, τ₃₁ = τ₁₃, τ₃₂ = τ₂₃).

dims selects which spatial directions (1 → x, 2 → y, 3 → z) enter the tensor: component τᵢⱼ is included only when both i and j are in dims. The default dims = (1, 2, 3) returns the full tensor, while e.g. dims = (1, 3) returns the 2D stress tensor in the x–z plane (τ₁₁, τ₃₃, τ₁₃). Components are always ordered diagonals-first, independently of the order of dims.

Each component can be wrapped in a Field and used with output writers, time-averaging, etc. Can also be called as StressTensor(grid, u, v, w; dims, collocate_diagonals) to build the components from individual velocity fields. Building the tensor from perturbation velocities (e.g. via perturbation_fields) yields the kinematic Reynolds stress tensor.

julia> using Oceananigans, Oceanostics

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

julia> model = NonhydrostaticModel(grid);

julia> τ = StressTensor(model);

julia> keys(τ)
(:τ₁₁, :τ₂₂, :τ₃₃, :τ₁₂, :τ₁₃, :τ₂₃)

julia> τ.τ₁₁
StressTensor KernelFunctionOperation at (Face, Center, Center)
├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── kernel_function: Oceanostics.FlowDiagnostics.StressTensorKernel{1, 1, false}
└── arguments: ("Field",)
└── computes: stress tensor component  τᵢⱼ = uᵢuⱼ

MixedLayerDepth(grid, args; criterion)

Returns the mixed layer depth defined as the depth at which criterion is true.

Defaults to DensityAnomalyCriterion where the depth is that at which the density is some threshold (defaults to 0.125kg/m³) higher than the surface density.

When DensityAnomalyCriterion is used, the arguments buoyancy_formulation and C should be supplied where buoyancy_formulation should be the buoyancy model, and C should be a named tuple of (; T, S), (; T) or (; S) (the latter two if the buoyancy model specifies a constant salinity or temperature).

julia> using Oceananigans, Oceanostics

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

julia> C = (; b = CenterField(grid));

julia> h = MixedLayerDepth(grid, BuoyancyTracer(), C);

julia> h isa KernelFunctionOperation
true

struct BuoyancyAnomalyCriterion{FT} <: Oceanostics.FlowDiagnostics.AbstractAnomalyCriterion

Defines the mixed layer to be the depth at which the buoyancy is more than threshold greater than the surface buoyancy (but the pertubaton is usually negative).

When this model is used, the arguments buoyancy_formulation and C should be supplied where C should be the named tuple (; b), with b the buoyancy tracer.

julia> using Oceanostics

julia> BuoyancyAnomalyCriterion(threshold = -1e-4)
BuoyancyAnomalyCriterion{Float64}(-0.0001)

struct DensityAnomalyCriterion{FT} <: Oceanostics.FlowDiagnostics.AbstractAnomalyCriterion

Defines the mixed layer to be the depth at which the density is more than threshold greater than the surface density.

When this model is used, the arguments buoyancy_formulation and C should be supplied where buoyancy_formulation should be the buoyancy model, and C should be a named tuple of (; T, S), (; T) or (; S) (the latter two if the buoyancy model specifies a constant salinity or temperature).

julia> using Oceanostics

julia> DensityAnomalyCriterion()
DensityAnomalyCriterion{Float64}(1020.0, 9.80665, 0.125)

BottomCellValue(diagnostic)

Returns the value of the given diagnostic at the bottom, which can be either the bottom of the domain (lowest vertical level) or an immersed bottom.

julia> using Oceananigans, Oceanostics

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

julia> c = CenterField(grid);

julia> set!(c, (x, y, z) -> z);

julia> cᵇ = BottomCellValue(c);

julia> Field(cᵇ) isa Field
true