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
\[e_k = \tfrac{1}{2} u_i u_i\]
and its prognostic equation is obtained by contracting the momentum equation with the velocity:
\[\partial_t e_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 momentum flux (negative of the 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. As throughout Oceanostics, the lower-case $e_k$ is the pointwise energy density and the upper-case $E_k = \int e_k \, \mathrm{d}V$ its volume integral; the potential energy page states the convention in full.
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 momentum flux ($\varepsilon_k = -\partial_j u_i \, \tau_{ij}$), and an isotropic version ($\varepsilon_k = 2\nu S_{ij} S_{ij}$) valid when the turbulence closure uses a single scalar viscosity. The two agree for a constant viscosity, and both are non-negative for a down-gradient closure. As in the momentum equation, $\tau_{ij}$ is the momentum flux Oceananigans' kernels carry (i.e. , the negative of the stress tensor; $-2\nu S_{ij}$ for a constant viscosity).
All diagnostics are computed at (Center, Center, Center).
Example
julia> using Oceananigans, Oceanosticsjulia> grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1));julia> model = NonhydrostaticModel(grid; buoyancy=BuoyancyTracer(), tracers=:b, closure=ScalarDiffusivity(ν=1e-4));julia> eₖ = KineticEnergyEquation.KineticEnergy(model)KineticEnergy 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")└── computes: kinetic energy ½uᵢuᵢjulia> εₖ = KineticEnergyEquation.KineticEnergyIsotropicDissipationRate(model)KineticEnergyIsotropicDissipationRate 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")└── computes: isotropic kinetic energy dissipation rate 2νSᵢⱼSᵢⱼjulia> wb = KineticEnergyEquation.PotentialEnergyConversion(model)PotentialToKineticEnergyConversion 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")└── computes: potential to kinetic energy conversion uᵢbᵢKinetic energy
Oceanostics.KineticEnergyEquation.KineticEnergy — Type
KineticEnergy(model; kwargs...)
Calculate the kinetic energy per unit mass eₖ = ½uᵢuᵢ of model.
julia> using Oceananigans, Oceanosticsjulia> grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1));julia> model = NonhydrostaticModel(grid);julia> KE = KineticEnergyEquation.KineticEnergy(model)KineticEnergy 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")└── computes: kinetic energy ½uᵢuᵢTendency
Oceanostics.KineticEnergyEquation.KineticEnergyTendency — Type
KineticEnergyTendency(model; location)
Return a KernelFunctionOperation that computes the tendency uᵢGᵢ of the kinetic energy eₖ, excluding the nonhydrostatic pressure contribution:
∂ₜeₖ = ½∂ₜuᵢ² = uᵢGᵢ - uᵢ∂ᵢpₙₕₛjulia> using Oceananigansjulia> grid = RectilinearGrid(size = (1, 1, 4), extent = (1, 1, 1));julia> model = NonhydrostaticModel(grid);julia> using Oceanostics.KineticEnergyEquation: KineticEnergyTendencyjulia> ke_tendency = KineticEnergyTendency(model)KineticEnergyTendency 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")└── computes: kinetic energy tendency uᵢGᵢ (excl. nonhydrostatic pressure)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 Oceananigansjulia> grid = RectilinearGrid(size = (1, 1, 4), extent = (1,1,1));julia> model = NonhydrostaticModel(grid);julia> using Oceanostics.KineticEnergyEquation: KineticEnergyAdvectionjulia> ADV = KineticEnergyAdvection(model)KineticEnergyAdvection 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")└── computes: kinetic energy advection uᵢ∂ⱼ(uᵢuⱼ)Stress (diffusive) term
Oceanostics.KineticEnergyEquation.KineticEnergyStress — Type
KineticEnergyStress(model; location)
Return a KernelFunctionOperation that computes the diffusive term of the eₖ prognostic equation:
DIFF = uᵢ∂ⱼτᵢⱼwhere uᵢ are the velocity components and τᵢⱼ is the viscous momentum flux (minus the stress tensor). It enters the kinetic energy budget as -DIFF, and splits into a transport divergence plus the dissipation rate εₖ (DissipationRate).
julia> using Oceananigans, Oceanosticsjulia> grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1));julia> model = NonhydrostaticModel(grid; closure=ScalarDiffusivity(ν=1e-4));julia> DIFF = KineticEnergyEquation.KineticEnergyStress(model)KineticEnergyStress 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")└── computes: kinetic energy stress/diffusion uᵢ∂ⱼτᵢⱼ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, Oceanosticsjulia> grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1));julia> model = NonhydrostaticModel(grid);julia> FORC = KineticEnergyEquation.KineticEnergyForcing(model)KineticEnergyForcing 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")└── computes: kinetic energy forcing uᵢFᵤᵢPressure redistribution
Oceanostics.KineticEnergyEquation.KineticEnergyPressureRedistribution — Type
KineticEnergyPressureRedistribution(
model;
velocities,
pressure,
location
)
Return a KernelFunctionOperation that computes the pressure redistribution term:
PR = uᵢ∂ᵢpwhere p is the pressure. By default p is taken to be the total pressure (nonhydrostatic + hydrostatic):
julia> using Oceananigansjulia> grid = RectilinearGrid(size = (1, 1, 4), extent = (1,1,1));julia> model = NonhydrostaticModel(grid);julia> using Oceanostics.KineticEnergyEquation: KineticEnergyPressureRedistributionjulia> ∇u⃗p = KineticEnergyPressureRedistribution(model)KineticEnergyPressureRedistribution 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")└── computes: kinetic energy pressure redistribution uᵢ∂ᵢpWe 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)KineticEnergyPressureRedistribution 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")└── computes: kinetic energy pressure redistribution uᵢ∂ᵢpPotential energy conversion
Oceanostics.KineticEnergyEquation.PotentialEnergyConversion — Type
PotentialEnergyConversion(
model;
velocities,
tracers,
location
)
Return a KernelFunctionOperation that computes the conversion of potential energy into kinetic energy, defined as
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 conversion will be calculated using the resolved velocities and tracers:
julia> using Oceananigansjulia> grid = RectilinearGrid(size = (1, 1, 4), extent = (1,1,1));julia> model = NonhydrostaticModel(grid; buoyancy=BuoyancyTracer(), tracers=:b);julia> using Oceanostics.KineticEnergyEquation: PotentialEnergyConversionjulia> wb = PotentialEnergyConversion(model)PotentialToKineticEnergyConversion 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")└── computes: potential to kinetic energy conversion uᵢbᵢIf we want to calculate only the turbulent conversion 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′ = PotentialEnergyConversion(model, velocities=(u=model.velocities.u, v=model.velocities.v, w=w′), tracers=(b=b′,))PotentialToKineticEnergyConversion 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")└── computes: potential to kinetic energy conversion uᵢbᵢDissipation rate
Oceanostics.KineticEnergyEquation.DissipationRate — Type
DissipationRate(model; U, V, W, location)
Calculate the kinetic energy dissipation rate, defined as
εₖ = -∂ⱼuᵢ ⋅ τᵢⱼwhere ∂ⱼuᵢ is the velocity gradient tensor and τᵢⱼ is the closure's viscous momentum flux (minus the stress tensor), in the same convention Oceananigans uses. That definition holds for any closure. For a Newtonian flux with a scalar viscosity, τᵢⱼ = -2νSᵢⱼ, it reduces to εₖ = 2νSᵢⱼSᵢⱼ ≥ 0, which KineticEnergyIsotropicDissipationRate computes directly; neither that reduction nor, for a flux that is not down-gradient, the sign carries over to other closures.
julia> using Oceananigans, Oceanosticsjulia> grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1));julia> model = NonhydrostaticModel(grid; closure=ScalarDiffusivity(ν=1e-4));julia> ε = KineticEnergyEquation.DissipationRate(model)KineticEnergyDissipationRate 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")└── computes: kinetic energy dissipation rate -∂ⱼuᵢ·τᵢⱼOceanostics.KineticEnergyEquation.KineticEnergyIsotropicDissipationRate — Type
KineticEnergyIsotropicDissipationRate(
u,
v,
w,
closure,
closure_fields,
model_fields,
clock;
location
)
Calculate the isotropic kinetic energy 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, Oceanosticsjulia> grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1));julia> model = NonhydrostaticModel(grid; closure=ScalarDiffusivity(ν=1e-4));julia> ε = KineticEnergyEquation.KineticEnergyIsotropicDissipationRate(model)KineticEnergyIsotropicDissipationRate 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")└── computes: isotropic kinetic energy dissipation rate 2νSᵢⱼSᵢⱼ