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: 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 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: @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))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 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 ∫Kˡ 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, (; ∫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: (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 = 1.331 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 (59.934 seconds)
[ Info: Executing initial time step...
[ Info: ... initial time step complete (3.777 seconds).
┌ Info: iter = 200, [006.78%] time = 8.131 seconds, Δt = 43.611 ms, walltime = 2.592 minutes, walltime / timestep = 378.194 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.821 minutes, walltime / timestep = 68.649 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.054 minutes, walltime / timestep = 70.002 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.268 minutes, walltime / timestep = 64.004 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.455 minutes, walltime / timestep = 56.261 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.610 minutes, walltime / timestep = 46.609 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.784 minutes, walltime / timestep = 52.207 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.962 minutes, walltime / timestep = 53.220 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.114 minutes, walltime / timestep = 45.578 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.288 minutes, walltime / timestep = 52.394 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.439 minutes, walltime / timestep = 45.028 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.609 minutes, walltime / timestep = 51.129 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.766 minutes, walltime / timestep = 47.205 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.938 minutes, walltime / timestep = 51.434 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.109 minutes, walltime / timestep = 51.425 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.264 minutes, walltime / timestep = 46.344 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.452 minutes, walltime / timestep = 56.595 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.641 minutes, walltime / timestep = 56.697 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.826 minutes, walltime / timestep = 55.321 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.633 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(∫Kˡ)/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"][:]∫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 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])dKˡdt = (∫Kˡ_t[i2] .- ∫Kˡ_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 = @. -dKˡ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...
7315.2 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[] = 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(∫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.