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(εₖˡ)∂ₜ∫eₖˡ = TimeDerivative(∫eₖˡ)TimeDerivative of 1×1×1 Field{Nothing, Nothing, Nothing} reduced over dims = (1, 2, 3) on RectilinearGrid on CPUWe use two NetCDF writers. A snapshot writer stores the 2D fields and a budget writer only the integrated scalars, both on TimeInterval(1).
using NCDatasetsfilename = "kelvin_helmholtz"simulation.output_writers[:nc] = NetCDFWriter(model, (; Ri, Q, b, w̄b̄, Πₖ, εₖˡ), filename=joinpath(@__DIR__, filename), schedule=TimeInterval(1), overwrite_files=true)simulation.output_writers[:budget] = NetCDFWriter(model, (; ∂ₜ∫eₖˡ, ∫w̄b̄, ∫Πₖ, ∫εₖˡ), filename=joinpath(@__DIR__, filename * "_budget"), schedule=TimeInterval(1), overwrite_files=true)NetCDFWriter scheduled on TimeInterval(1 second):
├── filepath: kelvin_helmholtz_budget.nc
├── dimensions: time(0), x_faa(256), x_caa(256), z_aaf(257), z_aac(256)
├── 4 outputs: (∂ₜ∫eₖˡ, ∫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.115 minutes, walltime / timestep = 0 seconds
└ |u⃗|ₘₐₓ = [1.00e+00, 0.00e+00, 6.68e-02] m/s, advective CFL = 0.22, diffusive CFL = 0.002, νₘₐₓ = 0.00025 m²/s
[ Info: ... simulation initialization complete (33.687 seconds)
[ Info: Executing initial time step...
[ Info: ... initial time step complete (1.079 seconds).
┌ Info: iter = 200, [006.78%] time = 8.131 seconds, Δt = 43.607 ms, walltime = 1.491 minutes, walltime / timestep = 112.775 ms
└ |u⃗|ₘₐₓ = [1.01e+00, 0.00e+00, 2.16e-02] m/s, advective CFL = 0.8, diffusive CFL = 0.0071, νₘₐₓ = 0.00025 m²/s
┌ Info: iter = 400, [013.76%] time = 16.507 seconds, Δt = 42.152 ms, walltime = 1.599 minutes, walltime / timestep = 32.466 ms
└ |u⃗|ₘₐₓ = [1.03e+00, 0.00e+00, 6.08e-02] m/s, advective CFL = 0.8, diffusive CFL = 0.0069, νₘₐₓ = 0.00025 m²/s
┌ Info: iter = 600, [020.38%] time = 24.459 seconds, Δt = 38.009 ms, walltime = 1.708 minutes, walltime / timestep = 32.605 ms
└ |u⃗|ₘₐₓ = [1.09e+00, 0.00e+00, 1.60e-01] m/s, advective CFL = 0.8, diffusive CFL = 0.0062, νₘₐₓ = 0.00025 m²/s
┌ Info: iter = 800, [026.13%] time = 31.354 seconds, Δt = 32.017 ms, walltime = 1.812 minutes, walltime / timestep = 31.195 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.05%] time = 37.256 seconds, Δt = 28.440 ms, walltime = 1.911 minutes, walltime / timestep = 29.876 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.65%] time = 42.778 seconds, Δt = 27.830 ms, walltime = 2.002 minutes, walltime / timestep = 27.196 ms
└ |u⃗|ₘₐₓ = [1.30e+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.28%] time = 48.342 seconds, Δt = 28.500 ms, walltime = 2.103 minutes, walltime / timestep = 30.286 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.02%] time = 54.029 seconds, Δt = 28.948 ms, walltime = 2.198 minutes, walltime / timestep = 28.593 ms
└ |u⃗|ₘₐₓ = [1.31e+00, 0.00e+00, 5.22e-01] m/s, advective CFL = 0.8, diffusive CFL = 0.0047, νₘₐₓ = 0.00025 m²/s
┌ Info: iter = 1800, [049.83%] time = 59.800 seconds, Δt = 29.709 ms, walltime = 2.288 minutes, walltime / timestep = 26.936 ms
└ |u⃗|ₘₐₓ = [1.29e+00, 0.00e+00, 4.97e-01] m/s, advective CFL = 0.8, diffusive CFL = 0.0049, νₘₐₓ = 0.00025 m²/s
┌ Info: iter = 2000, [054.70%] time = 1.094 minutes, Δt = 27.277 ms, walltime = 2.384 minutes, walltime / timestep = 28.860 ms
└ |u⃗|ₘₐₓ = [1.28e+00, 0.00e+00, 5.38e-01] m/s, advective CFL = 0.8, diffusive CFL = 0.0045, νₘₐₓ = 0.00025 m²/s
┌ Info: iter = 2200, [058.71%] time = 1.174 minutes, Δt = 26.928 ms, walltime = 2.473 minutes, walltime / timestep = 26.568 ms
└ |u⃗|ₘₐₓ = [1.40e+00, 0.00e+00, 5.30e-01] m/s, advective CFL = 0.8, diffusive CFL = 0.0044, νₘₐₓ = 0.00025 m²/s
┌ Info: iter = 2400, [063.04%] time = 1.261 minutes, Δt = 28.647 ms, walltime = 2.559 minutes, walltime / timestep = 25.760 ms
└ |u⃗|ₘₐₓ = [1.40e+00, 0.00e+00, 6.55e-01] m/s, advective CFL = 0.8, diffusive CFL = 0.0047, νₘₐₓ = 0.00025 m²/s
┌ Info: iter = 2600, [067.48%] time = 1.350 minutes, Δt = 26.805 ms, walltime = 2.648 minutes, walltime / timestep = 26.861 ms
└ |u⃗|ₘₐₓ = [1.26e+00, 0.00e+00, 4.94e-01] m/s, advective CFL = 0.8, diffusive CFL = 0.0044, νₘₐₓ = 0.00025 m²/s
┌ Info: iter = 2800, [072.06%] time = 1.441 minutes, Δt = 29.925 ms, walltime = 2.745 minutes, walltime / timestep = 29.127 ms
└ |u⃗|ₘₐₓ = [1.26e+00, 0.00e+00, 5.52e-01] m/s, advective CFL = 0.8, diffusive CFL = 0.0049, νₘₐₓ = 0.00025 m²/s
┌ Info: iter = 3000, [077.18%] time = 1.544 minutes, Δt = 31.088 ms, walltime = 2.839 minutes, walltime / timestep = 28.033 ms
└ |u⃗|ₘₐₓ = [1.21e+00, 0.00e+00, 4.05e-01] m/s, advective CFL = 0.8, diffusive CFL = 0.0051, νₘₐₓ = 0.00025 m²/s
┌ Info: iter = 3200, [082.15%] time = 1.643 minutes, Δt = 28.092 ms, walltime = 2.932 minutes, walltime / timestep = 27.977 ms
└ |u⃗|ₘₐₓ = [1.34e+00, 0.00e+00, 4.76e-01] m/s, advective CFL = 0.8, diffusive CFL = 0.0046, νₘₐₓ = 0.00025 m²/s
┌ Info: iter = 3400, [086.63%] time = 1.733 minutes, Δt = 28.687 ms, walltime = 3.019 minutes, walltime / timestep = 26.115 ms
└ |u⃗|ₘₐₓ = [1.34e+00, 0.00e+00, 5.44e-01] m/s, advective CFL = 0.8, diffusive CFL = 0.0047, νₘₐₓ = 0.00025 m²/s
┌ Info: iter = 3600, [091.27%] time = 1.825 minutes, Δt = 29.301 ms, walltime = 3.113 minutes, walltime / timestep = 28.260 ms
└ |u⃗|ₘₐₓ = [1.31e+00, 0.00e+00, 5.41e-01] m/s, advective CFL = 0.8, diffusive CFL = 0.0048, νₘₐₓ = 0.00025 m²/s
┌ Info: iter = 3800, [096.11%] time = 1.922 minutes, Δt = 29.745 ms, walltime = 3.211 minutes, walltime / timestep = 29.390 ms
└ |u⃗|ₘₐₓ = [1.31e+00, 0.00e+00, 4.95e-01] m/s, advective CFL = 0.8, diffusive CFL = 0.0049, νₘₐₓ = 0.00025 m²/s
[ Info: Simulation is stopping after running for 2.489 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 DatasetEvery budget record carries the tendency and the three source terms at the same time. The first record has no earlier state to difference against and is written as zero, so the budget starts from the second.
bud_filepath = simulation.output_writers[:budget].filepathds_bud = NCDataset(bud_filepath)nb = 2:length(ds_bud["time"])t_bud = ds_bud["time"][nb]deₖˡdt = ds_bud["∂ₜ∫eₖˡ"][nb]w̄b̄_bud = ds_bud["∫w̄b̄"][nb]Πₖ_bud = ds_bud["∫Πₖ"][nb]εₖˡ_bud = ds_bud["∫εₖˡ"][nb]close(ds_bud)closed DatasetResidual in sum-to-zero form: the negative tendency plus the three sources, so the plotted curves add to it
resid = @. -deₖˡdt + w̄b̄_bud - Πₖ_bud - εₖˡ_budPlotting
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...);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_bud, -deₖˡdt, label="−d(∫eₖˡ)/dt")lines!(ax_bud, t_bud, w̄b̄_bud, label="∫w̄b̄ dV")lines!(ax_bud, t_bud, -Πₖ_bud, label="−∫Πₖ dV")lines!(ax_bud, t_bud, -εₖˡ_bud, label="−∫εₖˡ dV")lines!(ax_bud, t_bud, 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.