Kinetic energy equation

The KineticEnergyEquation module provides diagnostics for every term in the resolved kinetic energy (KE) budget. The kinetic energy per unit mass is defined as

\[K = \tfrac{1}{2} u_i u_i\]

and its prognostic equation is obtained by contracting the momentum equation with the velocity:

\[\partial_t K = \underbrace{-u_i \partial_j(u_i u_j)}_{\text{advection}} + \underbrace{u_i \partial_j \tau_{ij}}_{\text{stress}} - \underbrace{u_i \partial_i p}_{\text{pressure}} + \underbrace{u_i b_i}_{\text{buoyancy}} + \underbrace{u_i F_{u_i}}_{\text{forcing}}\]

where $\tau_{ij}$ is the viscous/subgrid stress tensor, $p$ is pressure, $b_i$ is the buoyancy acceleration component in the $i$-th direction, and $F_{u_i}$ is the forcing on the $i$-th momentum equation.

This decomposition is essential for understanding how kinetic energy is generated (e.g. by buoyancy production or forcing), redistributed (by advection or pressure work), and removed (by viscous dissipation). The module also provides two formulations of the dissipation rate: a general one based on the full stress tensor ($\varepsilon = \partial_j u_i \cdot F_{ij}$), and an isotropic version ($\varepsilon = 2\nu S_{ij} S_{ij}$) valid when the turbulence closure uses a single scalar viscosity.

All diagnostics are computed at (Center, Center, Center).

Example

julia> using Oceananigans, Oceanostics

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

julia> model = NonhydrostaticModel(grid; buoyancy=BuoyancyTracer(), tracers=:b, closure=ScalarDiffusivity(ν=1e-4));

julia> KE = KineticEnergyEquation.KineticEnergy(model)
KernelFunctionOperation at (Center, Center, Center)
├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── kernel_function: kinetic_energy_ccc (generic function with 1 method)
└── arguments: ("Field", "Field", "Field")

julia> ε = KineticEnergyEquation.KineticEnergyIsotropicDissipationRate(model)
KernelFunctionOperation at (Center, Center, Center)
├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── kernel_function: isotropic_viscous_dissipation_rate_ccc (generic function with 1 method)
└── arguments: ("Field", "Field", "Field", "NamedTuple")

julia> wb = KineticEnergyEquation.BuoyancyProduction(model)
KernelFunctionOperation at (Center, Center, Center)
├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── kernel_function: uᵢbᵢᶜᶜᶜ (generic function with 1 method)
└── arguments: ("NamedTuple", "BuoyancyForce", "NamedTuple")

Kinetic energy

Oceanostics.KineticEnergyEquation.KineticEnergy — Type

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

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

source

KineticEnergy(model; kwargs...)

Calculate the kinetic energy of model.

julia> using Oceananigans, Oceanostics

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

julia> model = NonhydrostaticModel(grid);

julia> KE = KineticEnergyEquation.KineticEnergy(model)
KernelFunctionOperation at (Center, Center, Center)
├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── kernel_function: kinetic_energy_ccc (generic function with 1 method)
└── arguments: ("Field", "Field", "Field")

source

Tendency

Oceanostics.KineticEnergyEquation.KineticEnergyTendency — Type

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);

julia> using Oceanostics.KineticEnergyEquation: 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ᵢGᵢᶜᶜᶜ (generic function with 1 method)
└── arguments: ("Centered", "Nothing", "Nothing", "Nothing", "Nothing", "Nothing", "Nothing", "Nothing", "Oceananigans.Models.NonhydrostaticModels.BackgroundFields", "NamedTuple", "NamedTuple", "NamedTuple", "Nothing", "Nothing", "Clock", "NamedTuple")

source

Advection

Oceanostics.KineticEnergyEquation.KineticEnergyAdvection — Type

KineticEnergyAdvection(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);

julia> using Oceanostics.KineticEnergyEquation: KineticEnergyAdvection

julia> ADV = KineticEnergyAdvection(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: ("NamedTuple", "Centered")

source

Stress (diffusive) term

Oceanostics.KineticEnergyEquation.KineticEnergyStress — Type

KineticEnergyStress(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.

julia> using Oceananigans, Oceanostics

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

julia> model = NonhydrostaticModel(grid; closure=ScalarDiffusivity(ν=1e-4));

julia> DIFF = KineticEnergyEquation.KineticEnergyStress(model)
KernelFunctionOperation at (Center, Center, Center)
├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── kernel_function: uᵢ∂ⱼ_τᵢⱼᶜᶜᶜ (generic function with 1 method)
└── arguments: ("ScalarDiffusivity", "Nothing", "Clock", "NamedTuple", "Nothing")

source

Forcing

Oceanostics.KineticEnergyEquation.KineticEnergyForcing — Type

KineticEnergyForcing(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ᵢ).

julia> using Oceananigans, Oceanostics

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

julia> model = NonhydrostaticModel(grid);

julia> FORC = KineticEnergyEquation.KineticEnergyForcing(model)
KernelFunctionOperation at (Center, Center, Center)
├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── kernel_function: uᵢFᵤᵢᶜᶜᶜ (generic function with 1 method)
└── arguments: ("NamedTuple", "Clock", "NamedTuple")

source

Pressure redistribution

Oceanostics.KineticEnergyEquation.KineticEnergyPressureRedistribution — Type

KineticEnergyPressureRedistribution(
    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);

julia> using Oceanostics.KineticEnergyEquation: KineticEnergyPressureRedistribution

julia> ∇u⃗p = KineticEnergyPressureRedistribution(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: ("NamedTuple", "Field")

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 = KineticEnergyPressureRedistribution(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: ("NamedTuple", "Field")

source

Buoyancy production

Oceanostics.KineticEnergyEquation.BuoyancyProduction — Type

BuoyancyProduction(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; buoyancy=BuoyancyTracer(), tracers=:b);

julia> using Oceanostics.KineticEnergyEquation: BuoyancyProduction

julia> wb = BuoyancyProduction(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: ("NamedTuple", "BuoyancyForce", "NamedTuple")

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′ = BuoyancyProduction(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: ("NamedTuple", "BuoyancyForce", "NamedTuple")

source

Dissipation rate

Oceanostics.KineticEnergyEquation.DissipationRate — Type

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

Calculate the Kinetic Energy Dissipation Rate, defined as

ε = ∂ⱼuᵢ ⋅ Fᵢⱼ

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

julia> using Oceananigans, Oceanostics

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

julia> model = NonhydrostaticModel(grid; closure=ScalarDiffusivity(ν=1e-4));

julia> ε = KineticEnergyEquation.DissipationRate(model)
KernelFunctionOperation at (Center, Center, Center)
├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── kernel_function: viscous_dissipation_rate_ccc (generic function with 1 method)
└── arguments: ("Nothing", "NamedTuple", "NamedTuple")

source

Oceanostics.KineticEnergyEquation.KineticEnergyIsotropicDissipationRate — Type

KineticEnergyIsotropicDissipationRate(
    u,
    v,
    w,
    closure,
    closure_fields,
    model_fields,
    clock;
    location
)

Calculate the Viscous Dissipation Rate 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).

julia> using Oceananigans, Oceanostics

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

julia> model = NonhydrostaticModel(grid; closure=ScalarDiffusivity(ν=1e-4));

julia> ε = KineticEnergyEquation.KineticEnergyIsotropicDissipationRate(model)
KernelFunctionOperation at (Center, Center, Center)
├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── kernel_function: isotropic_viscous_dissipation_rate_ccc (generic function with 1 method)
└── arguments: ("Field", "Field", "Field", "NamedTuple")

source