Kelvin-Helmholtz instability

In this example we simulate a simple 2D Kelvin-Helmholtz instability and then use Oceanostics to close the volume-integrated kinetic-energy budget of the filtered flow.

Before starting, make sure you have the required packages installed for this example, which can be done with

using Pkg
pkg"add Oceananigans, Oceanostics, CairoMakie"

Model and simulation setup

using Oceananigans

We work with nondimensional quantities, following the standard nondimensionalization of the stratified shear layer (Kaminski and Smyth, 2019). We nondimensionalize the Boussinesq equations using the shear-layer half-width h as the length scale and the velocity scale U (half the velocity difference across the layer), so that time is measured in units of h / U. The flow is then governed by three nondimensional numbers — the Richardson number Ri₀, the Reynolds number Re = U h / ν, and the Prandtl number Pr = ν / κ — from which the viscosity ν and the buoyancy diffusivity κ follow:

U   = 1     # velocity scale (half the velocity difference across the shear layer)
h   = 1     # length scale (shear-layer half-width)
Ri₀ = 0.1   # Richardson number
Re  = 4e3   # Reynolds number (bounded by the grid; see the resolution note below)
Pr  = 1     # Prandtl number

ν = U * h / Re   # viscosity
κ = ν / Pr       # buoyancy diffusivity

0.00025

We begin by creating a model with this isotropic diffusivity and centered advection on a xz grid, using a buoyancy b as the active scalar. We make the box one wavelength of the most unstable Kelvin-Helmholtz mode wide (k_max = 0.4446 / h; Michalke, 1964), so that the perturbation we seed below fits periodically:

N = 256
k_max = 0.4446 / h   # most unstable KH wavenumber (Michalke, 1964)
Lx = 2π / k_max      # one most-unstable wavelength
Lz = 10
grid = RectilinearGrid(size=(N, N), x=(-Lx/2, +Lx/2), z=(-Lz/2, +Lz/2), topology=(Periodic, Flat, Bounded))

model = NonhydrostaticModel(grid; timestepper = :RungeKutta3,
                            advection = Centered(order=4), # A centered scheme is used here to minimize numerical dissipation
                            closure = ScalarDiffusivity(; ν, κ),
                            buoyancy = BuoyancyTracer(), tracers = :b)

NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 256×1×256 RectilinearGrid{Float64, Periodic, Flat, Bounded} on CPU with 3×0×3 halo
├── timestepper: RungeKutta3TimeStepper
├── advection scheme: Centered(order=4)
├── tracers: b
├── closure: ScalarDiffusivity{ExplicitTimeDiscretization}(ν=0.00025, κ=(b=0.00025,))
├── buoyancy: BuoyancyTracer with ĝ = NegativeZDirection()
└── coriolis: Nothing

We use hyperbolic tangent profiles with the same length scale h for both the shear flow and the stratification. The buoyancy jump B₀ = U² Ri₀ / h is chosen so that the gradient Richardson number N² / (∂u/∂z)² reaches its minimum value Ri₀ = 0.1 — below the classical stability threshold of 1/4 — at the center of the shear layer (z = 0), where the flow is most unstable. To kick off the instability we perturb the vertical velocity w with the most unstable mode sin(k_max x), localized to the shear layer by a Gaussian envelope exp(-z²) and given a random amplitude. We seed the random number generator so the perturbation — and hence the movie — is reproducible:

B₀ = U^2 * Ri₀ / h
perturbation_amplitude = 5e-2

shear_flow(x, z) = U * tanh(z / h)
stratification(x, z) = B₀ * tanh(z / h)
perturbation(x, z) = perturbation_amplitude * abs(randn()) * exp(-z^2) * sin(x * k_max - π)

using Random
Random.seed!(43)
set!(model, u=shear_flow, b=stratification, w=perturbation)

Next create an adaptive-time-step simulation using the model above. The initial time step is set conservatively from the horizontal grid spacing and velocity scale; the TimeStepWizard below adapts it as the flow evolves:

Δx = minimum_xspacing(grid)
simulation = Simulation(model, Δt = 0.2 * Δx / U, stop_time=120)
conjure_time_step_wizard!(simulation, IterationInterval(2), cfl=0.8, max_Δt=1)

Model diagnostics

We set-up a progress messenger using the TimedMessenger, which displays, among other information, the time step duration

using Oceanostics

progress = ProgressMessengers.TimedMessenger()
simulation.callbacks[:progress] = Callback(progress, IterationInterval(200))

Callback of Oceanostics.ProgressMessengers.TimedMessenger{Oceanostics.ProgressMessengers.AbstractProgressMessenger} on IterationInterval(200)

We can also define some useful diagnostics of the flow, starting with the RichardsonNumber

Ri = RichardsonNumber(model)

RichardsonNumber KernelFunctionOperation at (Center, Center, Face)
├── grid: 256×1×256 RectilinearGrid{Float64, Periodic, Flat, Bounded} on CPU with 3×0×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|²

We also set-up the QVelocityGradientTensorInvariant, which is usually used for visualizing vortices in the flow:

Q = QVelocityGradientTensorInvariant(model)

QVelocityGradientTensorInvariant KernelFunctionOperation at (Center, Center, Center)
├── grid: 256×1×256 RectilinearGrid{Float64, Periodic, Flat, Bounded} on CPU with 3×0×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ᵢⱼ)

Q is one of the velocity gradient tensor invariants and it measures the amount of vorticity versus the strain in the flow and, when it's positive, indicates a vortex. This method of vortex visualization is called the Q-criterion.

Filtered kinetic energy budget

Kelvin-Helmholtz billows draw kinetic energy from the mean shear and pass it down to ever-smaller scales, so this is a natural flow in which to look at a filtered kinetic-energy budget in the spirit of Aluie et al. (2018). We define a box filter whose width is comparable to the shear-layer half-width h and use it to build every term in the budget of the filtered kinetic energy $\overline{K} = \tfrac{1}{2}\overline{u}_i\overline{u}_i$. Volume-integrated — advection and pressure work integrate to zero, since the flow is periodic in x and w = 0 with free slip at the z walls — that budget reads

\[\frac{d}{dt} \int \overline{K}\, dV = \int \overline{w}\,\overline{b}\, dV - \int \Pi_K\, dV - \int \overline{\varepsilon}\, dV ,\]

with a buoyancy production $\overline{w}\,\overline{b}$ (the conversion between filtered kinetic and potential energy), the cross-scale kinetic-energy flux $\Pi_K$ to subfilter scales (KineticEnergyCrossScaleFlux), and viscous dissipation due to the filtered flow $\overline{\varepsilon}$ (FilteredKineticEnergyDissipationRate).

using Oceananigans.AbstractOperations: @at

A box filter is specified by its stencil size N in grid points rather than by a physical width, so we pick the odd N whose stencil spans roughly the shear-layer half-width h. The grid is slightly anisotropic here (Δx ≈ 0.11 h, Δz ≈ 0.078 h), so no single N matches h exactly in both directions; N = 11 brackets it, spanning 11Δx ≈ 1.2 h in x and 11Δz ≈ 0.86 h in z.

bfilter = BoxFilter(; dims=(1, 3), N=11, boundary=:shrink)  # stencil ≈ h wide, the shear-layer half-width

u, w = model.velocities.u, model.velocities.w
b = model.tracers.b
# Materialize each filtered field so the multi-direction filter takes its fast staged (separable)
# path; composing the raw `bfilter(u)` into `ū^2 + w̄^2` below would instead run it fused (see the
# filter performance notes and `check_filter_staging`).
ū, w̄, b̄ = Field(bfilter(u)), Field(bfilter(w)), Field(bfilter(b))

Kˡ = @at (Center, Center, Center) (ū^2 + w̄^2) / 2   # filtered kinetic energy ½ūᵢūᵢ
w̄b̄ = @at (Center, Center, Center) (w̄ * b̄)           # buoyancy production of the filtered flow
Πₖ = KineticEnergyCrossScaleFlux(model, bfilter; dims=(1, 3))
εˡ = FilteredKineticEnergyDissipationRate(model, bfilter)

FilteredKineticEnergyDissipationRate KernelFunctionOperation at (Center, Center, Center)
├── grid: 256×1×256 RectilinearGrid{Float64, Periodic, Flat, Bounded} on CPU with 3×0×3 halo
├── kernel_function: filtered_dissipation_rate_ccc (generic function with 1 method)
└── arguments: ("NamedTuple", "NamedTuple")
└── computes: filtered kinetic energy dissipation rate  εˡ = ∂ⱼūᵢ·τ̄ᵢⱼ

The budget only needs the (cheap) volume integrals of these terms:

∫Kˡ = Integral(Kˡ)
∫w̄b̄ = Integral(w̄b̄)
∫Πₖ = Integral(Πₖ)
∫εˡ = Integral(εˡ)

Integral of BinaryOperation at (Center, Center, Center) over dims (1, 2, 3)
└── operand: BinaryOperation at (Center, Center, Center)
    └── grid: 256×1×256 RectilinearGrid{Float64, Periodic, Flat, Bounded} on CPU with 3×0×3 halo

We use two NetCDF writers. A snapshot writer stores the 2D fields on a plain TimeInterval(1), while a budget writer stores only the integrated scalars on ConsecutiveIterations(TimeInterval(1)) — a second sample one model step after each output time — which lets us finite-difference ∫Kˡ across that single step to estimate d/dt, exactly as in the Two-dimensional turbulence example.

using NCDatasets
filename = "kelvin_helmholtz"

simulation.output_writers[:nc] = NetCDFWriter(model, (; Ri, Q, b, w̄b̄, Πₖ, εˡ),
                                              filename=joinpath(@__DIR__, filename),
                                              schedule=TimeInterval(1),
                                              overwrite_existing=true)

simulation.output_writers[:budget] = NetCDFWriter(model, (; ∫Kˡ, ∫w̄b̄, ∫Πₖ, ∫εˡ),
                                                  filename=joinpath(@__DIR__, filename * "_budget"),
                                                  schedule=ConsecutiveIterations(TimeInterval(1)),
                                                  overwrite_existing=true)

NetCDFWriter scheduled on ConsecutiveIterations(TimeInterval(1 second), 1):
├── filepath: kelvin_helmholtz_budget.nc
├── dimensions: time(0), x_faa(256), x_caa(256), z_aaf(257), z_aac(256)
├── 4 outputs: (∫Πₖ, ∫Kˡ, ∫εˡ, ∫w̄b̄)
├── array_type: Array{Float32}
├── file_splitting: NoFileSplitting
└── file size: 33.5 KiB

Run the simulation and process results

To run the simulation:

run!(simulation)

[ Info: Initializing simulation...
┌ Info: iter =      0,  [000.00%] time = 0 seconds,  Δt = 12.145 ms,  walltime = 1.772 minutes,  walltime / timestep = 0 seconds
└       |u⃗|ₘₐₓ = [1.00e+00,  0.00e+00,  6.93e-02] m/s,  advective CFL = 0.22,  diffusive CFL = 0.002,  νₘₐₓ = 0.00025 m²/s
[ Info:     ... simulation initialization complete (23.712 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (3.346 seconds).
┌ Info: iter =    200,  [006.78%] time = 8.131 seconds,  Δt = 43.611 ms,  walltime = 2.384 minutes,  walltime / timestep = 183.633 ms
└       |u⃗|ₘₐₓ = [1.01e+00,  0.00e+00,  2.20e-02] m/s,  advective CFL = 0.8,  diffusive CFL = 0.0071,  νₘₐₓ = 0.00025 m²/s
┌ Info: iter =    400,  [013.79%] time = 16.549 seconds,  Δt = 42.133 ms,  walltime = 2.588 minutes,  walltime / timestep = 61.000 ms
└       |u⃗|ₘₐₓ = [1.03e+00,  0.00e+00,  6.11e-02] m/s,  advective CFL = 0.8,  diffusive CFL = 0.0069,  νₘₐₓ = 0.00025 m²/s
┌ Info: iter =    600,  [020.41%] time = 24.497 seconds,  Δt = 37.979 ms,  walltime = 2.786 minutes,  walltime / timestep = 59.554 ms
└       |u⃗|ₘₐₓ = [1.09e+00,  0.00e+00,  1.61e-01] m/s,  advective CFL = 0.8,  diffusive CFL = 0.0062,  νₘₐₓ = 0.00025 m²/s
┌ Info: iter =    800,  [026.16%] time = 31.386 seconds,  Δt = 31.991 ms,  walltime = 2.966 minutes,  walltime / timestep = 53.921 ms
└       |u⃗|ₘₐₓ = [1.19e+00,  0.00e+00,  3.30e-01] m/s,  advective CFL = 0.8,  diffusive CFL = 0.0052,  νₘₐₓ = 0.00025 m²/s
┌ Info: iter =   1000,  [031.07%] time = 37.285 seconds,  Δt = 28.438 ms,  walltime = 3.131 minutes,  walltime / timestep = 49.593 ms
└       |u⃗|ₘₐₓ = [1.28e+00,  0.00e+00,  4.63e-01] m/s,  advective CFL = 0.8,  diffusive CFL = 0.0047,  νₘₐₓ = 0.00025 m²/s
┌ Info: iter =   1200,  [035.67%] time = 42.806 seconds,  Δt = 27.829 ms,  walltime = 3.277 minutes,  walltime / timestep = 43.783 ms
└       |u⃗|ₘₐₓ = [1.31e+00,  0.00e+00,  4.94e-01] m/s,  advective CFL = 0.8,  diffusive CFL = 0.0046,  νₘₐₓ = 0.00025 m²/s
┌ Info: iter =   1400,  [040.31%] time = 48.370 seconds,  Δt = 28.505 ms,  walltime = 3.442 minutes,  walltime / timestep = 49.383 ms
└       |u⃗|ₘₐₓ = [1.31e+00,  0.00e+00,  5.25e-01] m/s,  advective CFL = 0.8,  diffusive CFL = 0.0047,  νₘₐₓ = 0.00025 m²/s
┌ Info: iter =   1600,  [045.05%] time = 54.058 seconds,  Δt = 28.963 ms,  walltime = 3.605 minutes,  walltime / timestep = 49.154 ms
└       |u⃗|ₘₐₓ = [1.31e+00,  0.00e+00,  5.21e-01] m/s,  advective CFL = 0.8,  diffusive CFL = 0.0047,  νₘₐₓ = 0.00025 m²/s
┌ Info: iter =   1800,  [049.86%] time = 59.831 seconds,  Δt = 29.763 ms,  walltime = 3.758 minutes,  walltime / timestep = 45.679 ms
└       |u⃗|ₘₐₓ = [1.29e+00,  0.00e+00,  4.95e-01] m/s,  advective CFL = 0.8,  diffusive CFL = 0.0049,  νₘₐₓ = 0.00025 m²/s
┌ Info: iter =   2000,  [054.73%] time = 1.095 minutes,  Δt = 27.107 ms,  walltime = 3.922 minutes,  walltime / timestep = 49.144 ms
└       |u⃗|ₘₐₓ = [1.27e+00,  0.00e+00,  5.23e-01] m/s,  advective CFL = 0.8,  diffusive CFL = 0.0044,  νₘₐₓ = 0.00025 m²/s
┌ Info: iter =   2200,  [058.87%] time = 1.177 minutes,  Δt = 28.360 ms,  walltime = 4.072 minutes,  walltime / timestep = 45.035 ms
└       |u⃗|ₘₐₓ = [1.31e+00,  0.00e+00,  5.37e-01] m/s,  advective CFL = 0.8,  diffusive CFL = 0.0046,  νₘₐₓ = 0.00025 m²/s
┌ Info: iter =   2400,  [063.45%] time = 1.269 minutes,  Δt = 28.691 ms,  walltime = 4.236 minutes,  walltime / timestep = 49.327 ms
└       |u⃗|ₘₐₓ = [1.30e+00,  0.00e+00,  6.96e-01] m/s,  advective CFL = 0.8,  diffusive CFL = 0.0047,  νₘₐₓ = 0.00025 m²/s
┌ Info: iter =   2600,  [068.04%] time = 1.361 minutes,  Δt = 28.688 ms,  walltime = 4.391 minutes,  walltime / timestep = 46.367 ms
└       |u⃗|ₘₐₓ = [1.25e+00,  0.00e+00,  5.16e-01] m/s,  advective CFL = 0.8,  diffusive CFL = 0.0047,  νₘₐₓ = 0.00025 m²/s
┌ Info: iter =   2800,  [072.66%] time = 1.453 minutes,  Δt = 27.602 ms,  walltime = 4.559 minutes,  walltime / timestep = 50.373 ms
└       |u⃗|ₘₐₓ = [1.23e+00,  0.00e+00,  5.51e-01] m/s,  advective CFL = 0.8,  diffusive CFL = 0.0045,  νₘₐₓ = 0.00025 m²/s
┌ Info: iter =   3000,  [077.63%] time = 1.553 minutes,  Δt = 30.459 ms,  walltime = 4.736 minutes,  walltime / timestep = 53.330 ms
└       |u⃗|ₘₐₓ = [1.28e+00,  0.00e+00,  4.24e-01] m/s,  advective CFL = 0.8,  diffusive CFL = 0.005,  νₘₐₓ = 0.00025 m²/s
┌ Info: iter =   3200,  [082.50%] time = 1.650 minutes,  Δt = 31.069 ms,  walltime = 4.895 minutes,  walltime / timestep = 47.674 ms
└       |u⃗|ₘₐₓ = [1.26e+00,  0.00e+00,  4.54e-01] m/s,  advective CFL = 0.8,  diffusive CFL = 0.0051,  νₘₐₓ = 0.00025 m²/s
┌ Info: iter =   3400,  [087.44%] time = 1.749 minutes,  Δt = 29.653 ms,  walltime = 5.070 minutes,  walltime / timestep = 52.277 ms
└       |u⃗|ₘₐₓ = [1.29e+00,  0.00e+00,  4.68e-01] m/s,  advective CFL = 0.8,  diffusive CFL = 0.0049,  νₘₐₓ = 0.00025 m²/s
┌ Info: iter =   3600,  [092.16%] time = 1.843 minutes,  Δt = 28.074 ms,  walltime = 5.237 minutes,  walltime / timestep = 50.221 ms
└       |u⃗|ₘₐₓ = [1.33e+00,  0.00e+00,  5.42e-01] m/s,  advective CFL = 0.8,  diffusive CFL = 0.0046,  νₘₐₓ = 0.00025 m²/s
┌ Info: iter =   3800,  [096.81%] time = 1.936 minutes,  Δt = 29.429 ms,  walltime = 5.409 minutes,  walltime / timestep = 51.577 ms
└       |u⃗|ₘₐₓ = [1.30e+00,  0.00e+00,  5.39e-01] m/s,  advective CFL = 0.8,  diffusive CFL = 0.0048,  νₘₐₓ = 0.00025 m²/s
[ Info: Simulation is stopping after running for 3.765 minutes.
[ Info: Simulation time 2 minutes equals or exceeds stop time 2 minutes.

Now we'll read the snapshot fields using FieldTimeSeries

filepath = simulation.output_writers[:nc].filepath
Ri_t = FieldTimeSeries(filepath, "Ri")
Q_t  = FieldTimeSeries(filepath, "Q")
b_t  = FieldTimeSeries(filepath, "b")
w̄b̄_t = FieldTimeSeries(filepath, "w̄b̄")
Πₖ_t = FieldTimeSeries(filepath, "Πₖ")
εˡ_t = FieldTimeSeries(filepath, "εˡ")

ds = NCDataset(filepath)
times = ds["time"][:]
close(ds)

closed Dataset

The integrated budget scalars come in consecutive-iteration pairs (2k-1, 2k); a one-step finite difference inside each pair gives d(∫Kˡ)/dt, and each source term is evaluated at the pair midpoint.

bud_filepath = simulation.output_writers[:budget].filepath
ds_bud = NCDataset(bud_filepath)
times_bud = ds_bud["time"][:]
∫Kˡ_t     = ds_bud["∫Kˡ"][:]
∫w̄b̄_t     = ds_bud["∫w̄b̄"][:]
∫Πₖ_t     = ds_bud["∫Πₖ"][:]
∫εˡ_t     = ds_bud["∫εˡ"][:]
close(ds_bud)

i1 = 1:2:length(times_bud)-1   # primary snapshots
i2 = 2:2:length(times_bud)       # consecutive-iteration snapshots
Δt_pair = times_bud[i2] .- times_bud[i1]
t_pair = @. 0.5 * (times_bud[i1] + times_bud[i2])

dKˡdt   = (∫Kˡ_t[i2] .- ∫Kˡ_t[i1]) ./ Δt_pair
w̄b̄_pair = @. 0.5 * (∫w̄b̄_t[i1] + ∫w̄b̄_t[i2])
Πₖ_pair = @. 0.5 * (∫Πₖ_t[i1] + ∫Πₖ_t[i2])
εˡ_pair = @. 0.5 * (∫εˡ_t[i1] + ∫εˡ_t[i2])

120-element Vector{Float64}:
 0.004655505996197462
 0.004653243348002434
 0.004651123192161322
 0.004649084992706776
 0.0046471161767840385
 0.004645168781280518
 0.004643262829631567
 0.00464139599353075
 0.004639582708477974
 0.004637842997908592
 0.0046361954882740974
 0.004634657874703407
 0.004633246921002865
 0.004631973337382078
 0.004630850628018379
 0.004629889503121376
 0.004629101604223251
 0.00462849996984005
 0.004628099501132965
 0.004627917893230915
 0.004627975169569254
 0.004628290422260761
 0.004628878086805344
 0.004629748873412609
 0.0046309055760502815
 0.0046323444694280624
 0.004634065553545952
 0.004636092111468315
 0.004638508427888155
 0.004641497507691383
 0.00464536901563406
 0.004650571383535862
 0.004657688084989786
 0.004667400382459164
 0.004680432379245758
 0.0046975016593933105
 0.0047193425707519054
 0.004746786784380674
 0.004780841059982777
 0.004822683520615101
 0.004873526282608509
 0.004934326745569706
 0.005005443003028631
 0.005086369812488556
 0.005175638012588024
 0.005270998924970627
 0.005369942635297775
 0.005470274947583675
 0.005570547189563513
 0.005670209415256977
 0.00576937198638916
 0.00586813734844327
 0.005965640768408775
 0.006059448700398207
 0.006146084517240524
 0.006222194526344538
 0.0062851072289049625
 0.006332657299935818
 0.006363268941640854
 0.006376385688781738
 0.006371200084686279
 0.0063433656468987465
 0.006285693030804396
 0.006199266295880079
 0.006100062280893326
 0.006003114394843578
 0.005901101976633072
 0.0057772016152739525
 0.005639564711600542
 0.005488804541528225
 0.0053433505818247795
 0.0052452003583312035
 0.005200877785682678
 0.005207362584769726
 0.005245379637926817
 0.005259605590254068
 0.005218904465436935
 0.005129063036292791
 0.005009748507291079
 0.004867195151746273
 0.004716275259852409
 0.004584795329719782
 0.004501937422901392
 0.004454781301319599
 0.004406527616083622
 0.004346755333244801
 0.004282193724066019
 0.004229191690683365
 0.004191863350570202
 0.004160469397902489
 0.004135490395128727
 0.004108055029064417
 0.0040740082040429115
 0.004034143872559071
 0.003996304236352444
 0.003963105380535126
 0.003933528903871775
 0.003911738283932209
 0.0038972648326307535
 0.0038884151726961136
 0.0038823112845420837
 0.003876729402691126
 0.0038737780414521694
 0.003878389485180378
 0.0038928999565541744
 0.003915505483746529
 0.003945005126297474
 0.003977849148213863
 0.0040076663717627525
 0.0040298812091350555
 0.004045859910547733
 0.004059404134750366
 0.0040726568549871445
 0.004089072346687317
 0.004113162402063608
 0.004149493761360645
 0.004201256204396486
 0.004267738200724125
 0.004344635643064976
 0.004425111226737499

Residual in sum-to-zero form: the negative tendency plus the three sources, so the plotted curves add to it

resid = @. -dKˡdt + w̄b̄_pair - Πₖ_pair - εˡ_pair

Plotting

We now use Makie to create the figure and its axes

using CairoMakie

set_theme!(Theme(fontsize=24))
fig = Figure()

kwargs = (xlabel="x", ylabel="z", height=150, width=250)
ax1 = Axis(fig[2, 1]; title="Ri", kwargs...)
ax2 = Axis(fig[2, 2]; title="Q", kwargs...)
ax3 = Axis(fig[2, 3]; title="b", kwargs...);

Precompiling packages...
    578.3 ms  ✓ FilePathsGlobExt (serial)
  1 dependency successfully precompiled in 1 seconds
Precompiling packages...
    474.4 ms  ✓ DistancesChainRulesCoreExt (serial)
  1 dependency successfully precompiled in 0 seconds
Precompiling packages...
   1073.2 ms  ✓ TaylorSeriesIAExt (serial)
  1 dependency successfully precompiled in 1 seconds
Precompiling packages...
   6945.9 ms  ✓ OceananigansMakieExt (serial)
  1 dependency successfully precompiled in 7 seconds

Next we use Observables to lift the values and plot heatmaps and their colorbars

n = Observable(1)

Riₙ = @lift Ri_t[$n]
hm1 = heatmap!(ax1, Riₙ; colormap=:bwr, colorrange=(-1, +1))
Colorbar(fig[3, 1], hm1, vertical=false, height=8)

Qₙ  = @lift Q_t[$n]
hm2 = heatmap!(ax2, Qₙ; colormap=:inferno, colorrange=(0, 0.2))
Colorbar(fig[3, 2], hm2, vertical=false, height=8)

bₙ = @lift b_t[$n]
hm3 = heatmap!(ax3, bₙ; colormap=:balance, colorrange=(-B₀, +B₀))
Colorbar(fig[3, 3], hm3, vertical=false, height=8);

The second row shows the (local) budget terms as 2D fields: the buoyancy production w̄b̄, the cross-scale kinetic-energy flux Πₖ, and the filtered dissipation εˡ. Each gets a symmetric (or, for the sign-definite εˡ, one-sided) color range set from its own peak magnitude over the run.

maxabs(fts) = maximum(maximum(abs, interior(fts[k])) for k in 1:length(times))
wb_lim = maxabs(w̄b̄_t)
Π_lim  = maxabs(Πₖ_t)
ε_lim  = maxabs(εˡ_t)

ax4 = Axis(fig[4, 1]; title="w̄b̄", kwargs...)
ax5 = Axis(fig[4, 2]; title="Πₖ", kwargs...)
ax6 = Axis(fig[4, 3]; title="εˡ", kwargs...)

w̄b̄ₙ = @lift w̄b̄_t[$n]
hm4 = heatmap!(ax4, w̄b̄ₙ; colormap=:balance, colorrange=(-wb_lim, wb_lim))
Colorbar(fig[5, 1], hm4, vertical=false, height=8)

Πₖₙ = @lift Πₖ_t[$n]
hm5 = heatmap!(ax5, Πₖₙ; colormap=:balance, colorrange=(-Π_lim, Π_lim))
Colorbar(fig[5, 2], hm5, vertical=false, height=8)

εˡₙ = @lift εˡ_t[$n]
hm6 = heatmap!(ax6, εˡₙ; colormap=:magma, colorrange=(0, ε_lim))
Colorbar(fig[5, 3], hm6, vertical=false, height=8);

The bottom panel shows the volume-integrated filtered kinetic-energy budget. We plot the negative tendency −d(∫Kˡ)/dt together with its three sources: buoyancy production ∫w̄b̄ dV, the cross-scale flux −∫Πₖ dV, and the filtered dissipation −∫εˡ dV. With the tendency negated, the four curves sum to the residual.

ax_bud = Axis(fig[6, 1:3]; xlabel="Time", title="Filtered KE budget", height=140)
lines!(ax_bud, t_pair, -dKˡdt, label="−d(∫Kˡ)/dt")
lines!(ax_bud, t_pair, w̄b̄_pair, label="∫w̄b̄ dV")
lines!(ax_bud, t_pair, -Πₖ_pair, label="−∫Πₖ dV")
lines!(ax_bud, t_pair, -εˡ_pair, label="−∫εˡ dV")
lines!(ax_bud, t_pair, resid, label="residual", color=:black, linestyle=:dash)
axislegend(ax_bud; position=:lb, labelsize=10)

Makie.Legend()

Now we mark the time by placing a vertical line in the bottom panel and adding a helpful title

tₙ = @lift times[$n]
vlines!(ax_bud, tₙ, color=:black, linestyle=:dash)

title = @lift "Time = " * string(round(times[$n], digits=2))
fig[1, 1:3] = Label(fig, title, fontsize=24, tellwidth=false);

Finally, we adjust the figure dimensions to fit all the panels and record a movie

resize_to_layout!(fig)

@info "Animating..."
record(fig, filename * ".mp4", 1:length(times), framerate=10) do i
    n[] = i
end

"kelvin_helmholtz.mp4"

The bottom panel shows the volume-integrated filtered kinetic-energy budget. As the billows grow and overturn, the filtered flow mostly loses kinetic energy to potential energy (∫w̄b̄ dV < 0) and feeds the subfilter scales through the cross-scale flux (−∫Πₖ dV), while the filtered viscous dissipation ∫εˡ dV stays comparatively small at this Reynolds number. The residual (dashed), the sum of the negative tendency −d(∫Kˡ)/dt and the three source terms, stays small. As in the Two-dimensional turbulence example, the centered scheme contributes no numerical dissipation of its own, so the budget closes against the explicit ∫εˡ dV alone with a negligible residual.

This page was generated using Literate.jl.