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 Pkgpkg"add Oceananigans, Oceanostics, CairoMakie"Model and simulation setup
using OceananigansWe 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 numberRe = 4e3 # Reynolds number (bounded by the grid; see the resolution note below)Pr = 1 # Prandtl numberν = U * h / Re # viscosityκ = ν / Pr # buoyancy diffusivity0.00025We 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 = 256k_max = 0.4446 / h # most unstable KH wavenumber (Michalke, 1964)Lx = 2π / k_max # one most-unstable wavelengthLz = 10grid = 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:
│ ├── momentum: Centered(order=4)
│ └── b: Centered(order=4)
├── tracers: b
├── closure: ScalarDiffusivity{ExplicitTimeDiscretization}(ν=0.00025, κ=(b=0.00025,))
├── buoyancy: BuoyancyTracer with ĝ = NegativeZDirection()
└── coriolis: NothingWe 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₀ / hperturbation_amplitude = 5e-2shear_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 RandomRandom.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 Oceanosticsprogress = 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 GradientRichardsonNumber
Ri = GradientRichardsonNumber(model)GradientRichardsonNumber 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 (∂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 ½(ΩᵢⱼΩᵢⱼ - 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 $e_k^l = \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 e_k^l\, \mathrm{d}V = \int \overline{w}\,\overline{b}\, \mathrm{d}V - \int \Pi_k\, \mathrm{d}V - \int \varepsilon_k^l\, \mathrm{d}V ,\]
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 the viscous dissipation of the filtered flow $\varepsilon_k^l$ (FilteredKineticEnergyDissipationRate). Note that $\varepsilon_k^l$ is not $\overline{\varepsilon_k}$: their difference is the subfilter dissipation of the subfilter kinetic energy budget.
using Oceananigans.AbstractOperations: @atA 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-widthu, w = model.velocities.u, model.velocities.wb = 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))eₖˡ = @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:
∫eₖˡ = Integral(eₖˡ)∫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 haloWe 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 ∫eₖˡ across that single step to estimate d/dt, exactly as in the Two-dimensional turbulence example.
using NCDatasetsfilename = "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, (; ∫eₖˡ, ∫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: (∫Πₖ, ∫εₖˡ, ∫w̄b̄, ∫eₖˡ)
├── array_type: Array{Float32}
├── file_splitting: NoFileSplitting
└── file size: (file not yet created)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 = 2.773 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 (57.761 seconds)
[ Info: Executing initial time step...
[ Info: ... initial time step complete (3.107 seconds).
┌ Info: iter = 200, [006.78%] time = 8.131 seconds, Δt = 43.611 ms, walltime = 3.380 minutes, walltime / timestep = 182.099 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 = 3.571 minutes, walltime / timestep = 57.222 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 = 3.766 minutes, walltime / timestep = 58.374 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 = 3.941 minutes, walltime / timestep = 52.478 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 = 4.097 minutes, walltime / timestep = 46.884 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 = 4.239 minutes, walltime / timestep = 42.623 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 = 4.397 minutes, walltime / timestep = 47.492 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 = 4.561 minutes, walltime / timestep = 49.068 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 = 4.704 minutes, walltime / timestep = 42.971 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 = 4.862 minutes, walltime / timestep = 47.442 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 = 5.004 minutes, walltime / timestep = 42.547 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 = 5.169 minutes, walltime / timestep = 49.460 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 = 5.312 minutes, walltime / timestep = 42.823 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 = 5.471 minutes, walltime / timestep = 47.757 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 = 5.630 minutes, walltime / timestep = 47.681 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 = 5.772 minutes, walltime / timestep = 42.608 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.937 minutes, walltime / timestep = 49.461 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 = 6.096 minutes, walltime / timestep = 47.755 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 = 6.254 minutes, walltime / timestep = 47.604 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 4.172 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].filepathRi_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 DatasetThe integrated budget scalars come in consecutive-iteration pairs (2k-1, 2k); a one-step finite difference inside each pair gives d(∫eₖˡ)/dt, and each source term is evaluated at the pair midpoint.
bud_filepath = simulation.output_writers[:budget].filepathds_bud = NCDataset(bud_filepath)times_bud = ds_bud["time"][:]∫eₖˡ_t = ds_bud["∫eₖˡ"][:]∫w̄b̄_t = ds_bud["∫w̄b̄"][:]∫Πₖ_t = ds_bud["∫Πₖ"][:]∫εₖˡ_t = ds_bud["∫εₖˡ"][:]close(ds_bud)i1 = 1:2:length(times_bud)-1 # primary snapshotsi2 = 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])deₖˡdt = (∫eₖˡ_t[i2] .- ∫eₖˡ_t[i1]) ./ Δt_pairw̄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]);Residual in sum-to-zero form: the negative tendency plus the three sources, so the plotted curves add to it
resid = @. -deₖˡdt + w̄b̄_pair - Πₖ_pair - εₖˡ_pairPlotting
We now use Makie to create the figure and its axes
using CairoMakieset_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...
529.8 ms ✓ FilePathsGlobExt (serial)
1 dependency successfully precompiled in 1 seconds
Precompiling packages...
428.6 ms ✓ DistancesChainRulesCoreExt (serial)
1 dependency successfully precompiled in 0 seconds
Precompiling packages...
1016.9 ms ✓ TaylorSeriesIAExt (serial)
1 dependency successfully precompiled in 1 seconds
Precompiling packages...
6272.7 ms ✓ OceananigansMakieExt (serial)
1 dependency successfully precompiled in 6 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(∫eₖˡ)/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 kinetic energy budget", height=140)lines!(ax_bud, t_pair, -deₖˡdt, label="−d(∫eₖˡ)/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[] = iend"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(∫eₖˡ)/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.