Kelvin-Helmholtz instability
In this example we simulate a simple 2D Kelvin-Helmholtz instability and then use Oceanostics to close the volume-integrated kinetic and available potential energy budgets 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 horizontal 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: N = 19 spans 19Δx ≈ 1.05 h. The filter acts along x only, since the filtered APE diagnostics below accept only horizontal filters.
bfilter = BoxFilter(; dims=1, N=19) # stencil ≈ h wide, the shear-layer half-widthu, w = model.velocities.u, model.velocities.wb = model.tracers.b# Materialize each filtered field, so the filter runs once per compute rather than inside every term# built from it below.ū, 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 CPUFiltered available potential energy budget
The same filter gives the available potential energy (APE) of the filtered buoyancy, $e_a^l = e_a(\bar b, z)$ (FilteredAvailablePotentialEnergy). Its transport terms integrate to zero as well, since the walls carry no buoyancy flux either, so its volume-integrated budget reads
\[\frac{d}{dt} \int e_a^l\, \mathrm{d}V = -\int \overline{w}\,\overline{b_r}\, \mathrm{d}V - \int \Pi_a\, \mathrm{d}V - \int \varepsilon_a^l\, \mathrm{d}V + \int R^l\, \mathrm{d}V ,\]
with the conversion to filtered kinetic energy $\overline{w}\,\overline{b_r}$ (FilteredAvailablePotentialToKineticEnergyConversion), the cross-scale APE flux $\Pi_a$ (AvailablePotentialEnergyCrossScaleFlux), the filtered dissipation $\varepsilon_a^l$ (FilteredAvailablePotentialEnergyDissipationRate) and the reference tendency $R^l$ (FilteredAvailablePotentialEnergyReferenceTendency), derived in the Filtered available potential energy equation. The conversion has the same volume integral as the $\overline{w}\,\overline{b}$ of the kinetic-energy budget, since $\overline{w}$ has zero horizontal mean at every height. $R^l$ is the rate at which the sorted profile $b^\star$, which mixing reshapes, changes the filtered APE. Unlike its full-field counterpart it does not integrate to zero.
These diagnostics share one sorted column of the buoyancy through ProfileLookup, and $R^l$ integrates the time derivative of its profile. Reading a TimeDerivative does not update it, so a callback does, on the iteration before each output and at the output:
z✶_column = reference_height(model, method=VerticalSort())lookup = ProfileLookup(z✶_column)∂ₜb✶ = TimeDerivative(reference_buoyancy(z✶_column))simulation.callbacks[:∂ₜb✶] = Callback(∂ₜb✶, PrecedingIterations(TimeInterval(1)))eₐˡ = FilteredAvailablePotentialEnergy(model, bfilter; method=lookup)w̄b̄ᵣ = FilteredAvailablePotentialToKineticEnergyConversion(model, bfilter; method=lookup)Πₐ = AvailablePotentialEnergyCrossScaleFlux(model, bfilter; dims=(1, 3), method=lookup)εₐˡ = FilteredAvailablePotentialEnergyDissipationRate(model, bfilter; method=lookup)Rˡ = FilteredAvailablePotentialEnergyReferenceTendency(model, bfilter, ∂ₜb✶; method=lookup)∫w̄b̄ᵣ = Integral(w̄b̄ᵣ)∫Πₐ = Integral(Πₐ)∫εₐˡ = Integral(εₐˡ)∫Rˡ = Integral(Rˡ)∂ₜ∫eₐˡ = TimeDerivative(Integral(eₐˡ))TimeDerivative of 1×1×1 Field{Nothing, Nothing, Nothing} reduced over dims = (1, 2, 3) on RectilinearGrid on CPUA single NetCDF writer stores, every time unit, the six fields the movie shows and the tendencies and volume integrals of the two budgets:
using NCDatasetsfilename = "kelvin_helmholtz"outputs = (; Q, Πₖ, εₖˡ, Ri, Πₐ, εₐˡ, ∂ₜ∫eₖˡ, ∫w̄b̄, ∫Πₖ, ∫εₖˡ, ∂ₜ∫eₐˡ, ∫w̄b̄ᵣ, ∫Πₐ, ∫εₐˡ, ∫Rˡ)simulation.output_writers[:nc] = NetCDFWriter(model, outputs, filename=joinpath(@__DIR__, filename), schedule=TimeInterval(1), overwrite_files=true)NetCDFWriter scheduled on TimeInterval(1 second):
├── filepath: kelvin_helmholtz.nc
├── dimensions: time(0), x_faa(256), x_caa(256), z_aaf(257), z_aac(256)
├── 15 outputs: (Q, Πₖ, εₖˡ, Ri, Πₐ, εₐˡ, ∂ₜ∫eₖˡ, ∫w̄b̄, ∫Πₖ, ∫εₖˡ, ∂ₜ∫eₐˡ, ∫w̄b̄ᵣ, ∫Πₐ, ∫εₐˡ, ∫Rˡ)
├── 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.226 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 (32.324 seconds)
[ Info: Executing initial time step...
[ Info: ... initial time step complete (816.610 ms).
┌ Info: iter = 200, [006.78%] time = 8.131 seconds, Δt = 43.607 ms, walltime = 2.561 minutes, walltime / timestep = 100.350 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 = 2.647 minutes, walltime / timestep = 26.048 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 = 2.732 minutes, walltime / timestep = 25.509 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 = 2.818 minutes, walltime / timestep = 25.681 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 = 2.895 minutes, walltime / timestep = 23.158 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.969 minutes, walltime / timestep = 22.119 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 = 3.046 minutes, walltime / timestep = 23.127 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 = 3.128 minutes, walltime / timestep = 24.607 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 = 3.204 minutes, walltime / timestep = 22.872 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 = 3.284 minutes, walltime / timestep = 23.981 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 = 3.358 minutes, walltime / timestep = 22.087 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 = 3.436 minutes, walltime / timestep = 23.521 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 = 3.512 minutes, walltime / timestep = 22.747 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 = 3.596 minutes, walltime / timestep = 25.143 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 = 3.674 minutes, walltime / timestep = 23.400 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 = 3.754 minutes, walltime / timestep = 24.057 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.831 minutes, walltime / timestep = 22.958 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.913 minutes, walltime / timestep = 24.834 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.996 minutes, walltime / timestep = 24.690 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.152 minutes.
[ Info: Simulation time 2 minutes equals or exceeds stop time 2 minutes.
Now we'll read the output back as FieldTimeSeries
filepath = simulation.output_writers[:nc].filepathQ_t = FieldTimeSeries(filepath, "Q")Πₖ_t = FieldTimeSeries(filepath, "Πₖ")εₖˡ_t = FieldTimeSeries(filepath, "εₖˡ")Ri_t = FieldTimeSeries(filepath, "Ri")Πₐ_t = FieldTimeSeries(filepath, "Πₐ")εₐˡ_t = FieldTimeSeries(filepath, "εₐˡ")times = Q_t.times0.0:1.0:120.0Every record carries both tendencies and their source terms at the same time, and each volume integral comes back as a FieldTimeSeries of a single point, which interior turns into a vector. The first record of a tendency has no earlier state to difference against and is written as zero, so the budget starts from the second.
nb = 2:length(times)t_bud = times[nb]budget_term(name) = interior(FieldTimeSeries(filepath, name), 1, 1, 1, nb)deₖˡdt = budget_term("∂ₜ∫eₖˡ")w̄b̄_bud = budget_term("∫w̄b̄")Πₖ_bud = budget_term("∫Πₖ")εₖˡ_bud = budget_term("∫εₖˡ")deₐˡdt = budget_term("∂ₜ∫eₐˡ")w̄b̄ᵣ_bud = budget_term("∫w̄b̄ᵣ")Πₐ_bud = budget_term("∫Πₐ")εₐˡ_bud = budget_term("∫εₐˡ")Rˡ_bud = budget_term("∫Rˡ");┌ Warning: Reading boundary conditions from NetCDF files is not supported for FieldTimeSeries. Using default FieldBoundaryConditions for `grid` and `location`.
└ @ OceananigansNCDatasetsExt ~/.julia/packages/Oceananigans/EhPjs/ext/OceananigansNCDatasetsExt/output_readers.jl:44
┌ Warning: Reading boundary conditions from NetCDF files is not supported for FieldTimeSeries. Using default FieldBoundaryConditions for `grid` and `location`.
└ @ OceananigansNCDatasetsExt ~/.julia/packages/Oceananigans/EhPjs/ext/OceananigansNCDatasetsExt/output_readers.jl:44
┌ Warning: Reading boundary conditions from NetCDF files is not supported for FieldTimeSeries. Using default FieldBoundaryConditions for `grid` and `location`.
└ @ OceananigansNCDatasetsExt ~/.julia/packages/Oceananigans/EhPjs/ext/OceananigansNCDatasetsExt/output_readers.jl:44
┌ Warning: Reading boundary conditions from NetCDF files is not supported for FieldTimeSeries. Using default FieldBoundaryConditions for `grid` and `location`.
└ @ OceananigansNCDatasetsExt ~/.julia/packages/Oceananigans/EhPjs/ext/OceananigansNCDatasetsExt/output_readers.jl:44
┌ Warning: Reading boundary conditions from NetCDF files is not supported for FieldTimeSeries. Using default FieldBoundaryConditions for `grid` and `location`.
└ @ OceananigansNCDatasetsExt ~/.julia/packages/Oceananigans/EhPjs/ext/OceananigansNCDatasetsExt/output_readers.jl:44
┌ Warning: Reading boundary conditions from NetCDF files is not supported for FieldTimeSeries. Using default FieldBoundaryConditions for `grid` and `location`.
└ @ OceananigansNCDatasetsExt ~/.julia/packages/Oceananigans/EhPjs/ext/OceananigansNCDatasetsExt/output_readers.jl:44
┌ Warning: Reading boundary conditions from NetCDF files is not supported for FieldTimeSeries. Using default FieldBoundaryConditions for `grid` and `location`.
└ @ OceananigansNCDatasetsExt ~/.julia/packages/Oceananigans/EhPjs/ext/OceananigansNCDatasetsExt/output_readers.jl:44
┌ Warning: Reading boundary conditions from NetCDF files is not supported for FieldTimeSeries. Using default FieldBoundaryConditions for `grid` and `location`.
└ @ OceananigansNCDatasetsExt ~/.julia/packages/Oceananigans/EhPjs/ext/OceananigansNCDatasetsExt/output_readers.jl:44
┌ Warning: Reading boundary conditions from NetCDF files is not supported for FieldTimeSeries. Using default FieldBoundaryConditions for `grid` and `location`.
└ @ OceananigansNCDatasetsExt ~/.julia/packages/Oceananigans/EhPjs/ext/OceananigansNCDatasetsExt/output_readers.jl:44
Residuals in sum-to-zero form: each negative tendency plus its sources, so the plotted curves add to it
eₖˡ_resid = @. -deₖˡdt + w̄b̄_bud - Πₖ_bud - εₖˡ_budeₐˡ_resid = @. -deₐˡdt - w̄b̄ᵣ_bud - Πₐ_bud - εₐˡ_bud + Rˡ_budPlotting
We build the figure in four rows, as in the Rayleigh-Taylor instability example. The top two hold snapshots, three panels each: Q and the kinetic-energy terms Πₖ and εₖˡ in the first row, and Ri and the APE terms Πₐ and εₐˡ in the second. The bottom two rows hold the two volume-integrated filtered budgets.
using CairoMakieset_theme!(Theme(fontsize=18))fig = Figure(size=(1500, 1300))n = Observable(1);Q and Ri keep colormaps of their own. The budget terms get a diverging colormap centered on zero, and the two cross-scale fluxes share one color range, as do the two dissipations, so the kinetic-energy and APE terms stacked in each column compare directly.
Πlim = 0.5 * max(maximum(abs, Πₖ_t), maximum(abs, Πₐ_t))εlim = 0.5 * max(maximum(abs, εₖˡ_t), maximum(abs, εₐˡ_t))panels = [(1, 1, Q_t, "velocity gradient invariant, Q", :inferno, (0, 0.2)), (1, 2, Πₖ_t, "cross-scale KE flux, Πₖ", :balance, (-Πlim, Πlim)), (1, 3, εₖˡ_t, "filtered KE dissipation, εₖˡ", :balance, (-εlim, εlim)), (2, 1, Ri_t, "gradient Richardson number, Ri", :bwr, (-1, 1)), (2, 2, Πₐ_t, "cross-scale APE flux, Πₐ", :balance, (-Πlim, Πlim)), (2, 3, εₐˡ_t, "filtered APE dissipation, εₐˡ", :balance, (-εlim, εlim))]for (row, column, series, label, colormap, colorrange) in panels ax = Axis(fig[row + 1, 2column - 1]; title=label, xlabel="x", ylabel="z", aspect=DataAspect()) hm = heatmap!(ax, @lift(series[$n]); colormap, colorrange) Colorbar(fig[row + 1, 2column], hm)endThe third row 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. The bottom row does the same for the filtered APE budget, whose reference tendency ∫Rˡ dV is a fourth source.
ax_bud = Axis(fig[4, 1:6]; xlabel="Time", title="Filtered kinetic energy budget")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, eₖˡ_resid, label="residual", color=:black, linestyle=:dash)axislegend(ax_bud; position=:lb, labelsize=12)ax_ape = Axis(fig[5, 1:6]; xlabel="Time", title="Filtered available potential energy budget")lines!(ax_ape, t_bud, -deₐˡdt, label="−d(∫eₐˡ)/dt")lines!(ax_ape, t_bud, -w̄b̄ᵣ_bud, label="−∫w̄b̄ᵣ dV")lines!(ax_ape, t_bud, -Πₐ_bud, label="−∫Πₐ dV")lines!(ax_ape, t_bud, -εₐˡ_bud, label="−∫εₐˡ dV")lines!(ax_ape, t_bud, Rˡ_bud, label="∫Rˡ dV")lines!(ax_ape, t_bud, eₐˡ_resid, label="residual", color=:black, linestyle=:dash)axislegend(ax_ape; position=:lb, labelsize=12, nbanks=2)Makie.Legend()Finally, we mark the time with a vertical line in both budget panels, add a title and record a movie
for ax in (ax_bud, ax_ape) vlines!(ax, @lift(times[$n]), color=:black, linestyle=:dash)endtitle = @lift "Kelvin-Helmholtz instability, t = " * string(round(times[$n], digits=2))fig[1, 1:6] = Label(fig, title, fontsize=22, tellwidth=false)@info "Animating..."record(fig, filename * ".mp4", 1:length(times), framerate=10) do i n[] = iend"kelvin_helmholtz.mp4"The first budget 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). The cross-scale flux sends kinetic energy to the subfilter scales while the billows grow (∫Πₖ dV > 0) and turns mostly inverse once they have overturned, after t ≈ 40. The filtered viscous dissipation ∫εₖˡ dV stays comparatively small at this Reynolds number.
The bottom panel shows the filtered APE budget, whose conversion is that of the kinetic-energy budget with the opposite sign. The filtered APE passes energy on to the subfilter scales throughout the run (∫Πₐ dV ≥ 0), so after t ≈ 40 kinetic energy moves up across the filter scale while APE keeps moving down. The diffusive sink ∫εₐˡ dV is mostly smaller than ∫Πₐ dV and peaks later. As the billows mix, the reference profile changes and ∫Rˡ dV removes filtered APE at a rate comparable to ∫εₐˡ dV, with both peaking around t ≈ 70. Since the full-field ∫R dV vanishes, the subfilter scales gain what ∫Rˡ dV removes.
In both panels the residual (dashed), the sum of the negative tendency and the source terms, stays small. As in the Two-dimensional turbulence example, the centered scheme contributes no numerical dissipation of its own, so both budgets close against the explicit dissipation rates alone.
This page was generated using Literate.jl.