Kinetic energy equation

KineticEnergy(model; kwargs...)

Calculate the kinetic energy per unit mass eₖ = ½uᵢuᵢ of model. On a HydrostaticFreeSurfaceModel, which steps u and v and diagnoses w from continuity, this is ½(u² + v²).

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ᵢ

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ₙₕₛ

On a HydrostaticFreeSurfaceModel the sum runs over u and v only, and Gᵢ is Oceananigans' hydrostatic tendency, which leaves out the part of the barotropic pressure gradient the free surface treats implicitly (all of it for an ImplicitFreeSurface or a SplitExplicitFreeSurface), so the change of eₖ over a step is uₕGₕ - uₕ⋅∇ₕ(gη) there.

On either model Gᵢ is the explicit tendency. A closure with a VerticallyImplicitTimeDiscretization (CATKE's default, and an option of ScalarDiffusivity and VerticalScalarDiffusivity) takes the ∂zu and ∂zv parts of its vertical fluxes out of Gᵢ and applies them in the implicit step, so this term, like KineticEnergyStress, leaves out the vertical viscous work the model does. The change of eₖ over a time step then differs from uᵢGᵢ - uᵢ∂ᵢpₙₕₛ by that work, which DissipationRate still accounts for. On an ImmersedBoundaryGrid Gᵢ also carries the stress through the immersed faces, which KineticEnergyStress leaves out.

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", "BackgroundFields", "NamedTuple", "NamedTuple", "NamedTuple", "Nothing", "Nothing", "Clock", "NamedTuple")└── computes: kinetic energy tendency  uᵢGᵢ (excl. nonhydrostatic pressure)

KineticEnergyAdvection(model; velocities, location)

Return a KernelFunctionOperation that computes the advection term, defined as

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

On a HydrostaticFreeSurfaceModel it is uₕ⋅(u⃗⋅∇)uₕ over u and v only, built on the kernels (U_dot_∇u, U_dot_∇v) the hydrostatic momentum equations advect with, so it follows the model's scheme, VectorInvariant included. The hydrostatic metric terms that a flux-form scheme picks up on a curvilinear grid are left out, as they are by UMomentumEquation.Advection; they vanish on a RectilinearGrid and for a VectorInvariant scheme.

By default, the advection will be calculated using the resolved velocities and users can 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ⱼ)

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

A HydrostaticFreeSurfaceModel has no w-momentum equation, so on that model the sum runs over i = 1, 2 only and DIFF is the stress term of the budget of ½(u² + v²). The fluxes τ₁ⱼ and τ₂ⱼ are still computed with the diagnostic w, as in the model's own u and v tendencies.

∂ⱼτᵢⱼ is Oceananigans' own explicit viscous term, with the closure's time discretization applied. A closure with a VerticallyImplicitTimeDiscretization (CATKE's default, and an option of ScalarDiffusivity and VerticalScalarDiffusivity) takes the ∂zu and ∂zv parts of its vertical fluxes out of that term and applies them in the implicit step, so DIFF then leaves out the vertical viscous work the model does, as KineticEnergyTendency does, while DissipationRate keeps the full explicit flux. Under such a closure ∫DIFF dV = ∫εₖ dV no longer holds, and an integrated budget should be written with εₖ rather than with DIFF.

On an ImmersedBoundaryGrid the tendency kernels also subtract the stress through the immersed faces, immersed_∂ⱼτᵢⱼ, which DIFF does not include. It vanishes under the default no-flux immersed boundary conditions, but with an immersed ValueBoundaryCondition or a drag FluxBoundaryCondition on the velocities the wall-stress work is left out of the budget; the momentum modules expose it per component as ImmersedViscousDissipation.

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: ("NamedTuple", "ScalarDiffusivity", "Nothing", "Clock", "NamedTuple", "Nothing")└── computes: kinetic energy stress/diffusion  uᵢ∂ⱼτᵢⱼ

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ᵢ). On a HydrostaticFreeSurfaceModel the sum runs over u and v only, since w has no momentum equation, and so no forcing, there.

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ᵤᵢ

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 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ᵢ∂ᵢp

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)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ᵢ∂ᵢp

On a HydrostaticFreeSurfaceModel the total pressure is p′ + gη: pressure defaults to the hydrostatic pressure anomaly model.pressure.pHY′, the pressure the hydrostatic momentum tendencies carry, and the work of the surface pressure gη of free_surface (model.free_surface by default) is added, so that PR is the work of the total pressure on both model types; on a NonhydrostaticModel with a free surface gη is the surface value of pNHS, and the default there already contains it. free_surface = nothing leaves the work of the anomaly alone. The sum runs over all three components with the diagnosed w, although the hydrostatic tendencies have no w row. The z row is w∂zp′ = wb by hydrostatic balance, and exactly so on Oceananigans' grid away from the top face, since p′ is integrated from the same face buoyancy that PotentialEnergyConversion multiplies by w. So -PR + wb is the work of the horizontal gradients, the one the model applies, and the budget reads ∂ₜeₖ = -uᵢ∂ⱼ(uⱼuᵢ) - uᵢ∂ⱼτᵢⱼ - PR + wb + uᵢFᵤᵢ on both model types; the two-component work of the anomaly alone is PR - uₕ⋅∇ₕ(gη) - wb.

On a NonhydrostaticModel that has taken a time step, Oceananigans (0.106 to 0.116 at least) leaves the halos of pNHS at the pressure times Δt: it corrects the velocities with p Δt and then divides the interior alone by Δt. The default p is a computed Field, which fills its own halos when the term is computed, so it is not affected. A pressure passed by hand, such as model.pressures.pNHS above, is differentiated across its own halos, so in the cells next to a periodic boundary the term is wrong until they are refilled with fill_halo_regions!(model.pressures.pNHS) right before it is evaluated. The error is a redistribution and barely changes a volume integral.

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.

On a HydrostaticFreeSurfaceModel, w is the vertical velocity diagnosed from continuity. Buoyancy enters the hydrostatic momentum equations through the hydrostatic pressure, so wb is not one of the terms of uₕGₕ there; it is the conversion that the horizontal pressure work -uₕ⋅∇ₕp′ contains, since ∂p′/∂z = b. KineticEnergyPressureRedistribution contracts all three components on that model, so that its z row is this wb and the budget reads -PR + wb on both model types.

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ᵢ

DissipationRate(
    model;
    U,
    V,
    W,
    perturbation_boundary_conditions,
    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 on a NonhydrostaticModel 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.

A HydrostaticFreeSurfaceModel has no w-momentum equation, so on that model the sum runs over i = 1, 2 only and εₖ is the dissipation rate of ½(u² + v²). The fluxes τ₁ⱼ and τ₂ⱼ are still computed with the diagnostic w, as in the model's own u and v tendencies. The reduction to 2νSᵢⱼSᵢⱼ needs the w row, so there εₖ is not that, and the isotropic rate, which keeps the row, is not a term of the hydrostatic budget. ∫εₖ dV is still the viscous sink of ½(u² + v²). Pointwise, εₖ ≥ 0 with a horizontal or vertical formulation (HorizontalScalarDiffusivity, VerticalScalarDiffusivity, CATKE), whose τ is built from u and v alone. With a three-dimensional ScalarDiffusivity, whose τ₁₃ carries ∂ₓw, the cross term ν ∂zu ∂ₓw has no w row to complete its square, so εₖ is indefinite cell by cell, by a term of order aspect ratio squared relative to ν(∂zu)², which is the order the hydrostatic approximation neglects elsewhere.

The fluxes are the closure's explicit viscous_flux_*, with no time discretization applied, so under a closure with a VerticallyImplicitTimeDiscretization (CATKE's default, and an option of ScalarDiffusivity and VerticalScalarDiffusivity) εₖ still includes the vertical shear term that the model applies in the implicit step. KineticEnergyStress and KineticEnergyTendency leave that term out, so ∫DIFF dV = ∫εₖ dV fails there and an integrated budget should be written with εₖ.

On an ImmersedBoundaryGrid the same bare fluxes are not masked across immersed faces, as the _viscous_flux_* behind ∂ⱼτᵢⱼ are, so εₖ there includes the shear across faces the model treats as solid and ∫DIFF dV = ∫εₖ dV does not hold either.

U, V and W are means removed from the velocities before εₖ is computed, which makes it the dissipation rate of the turbulent kinetic energy. Each perturbation is a computed Field whose boundary conditions are zero flux at every wall unless perturbation_boundary_conditions sets them, keyed by velocity. At a wall where a velocity has a value or gradient condition, the zero-flux default drops the wall shear from the cells next to it, so pass the condition the perturbation satisfies there. For a no-slip bottom and a horizontal mean, that is perturbation_boundary_conditions = (u = FieldBoundaryConditions(bottom = ValueBoundaryCondition(0)),). See Oceanostics.perturbation_fields for the details.

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ᵢ·τᵢⱼ

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). It keeps the w row on both model types, so on a HydrostaticFreeSurfaceModel it is not the dissipation term of the budget of ½(u² + v²), which DissipationRate computes: it includes the dissipation of ½w², which that budget does not carry.

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