Kelvin-Helmholtz instability
This example simulates a simple 2D Kelvin-Helmholtz instability and is based on the similar Oceananigans example.
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, Rasters"Model and simulation setup
We begin by creating a model with an isotropic diffusivity and fifth-order advection on a xz 128² grid using a buoyancy b as the active scalar. We'll work here with nondimensional quantities.
using Oceananigans
N = 128
L = 10
grid = RectilinearGrid(size=(N, N), x=(-L/2, +L/2), z=(-L/2, +L/2), topology=(Periodic, Flat, Bounded))
model = NonhydrostaticModel(; grid, timestepper = :RungeKutta3,
advection = UpwindBiased(order=5),
closure = ScalarDiffusivity(ν=2e-5, κ=2e-5),
buoyancy = BuoyancyTracer(), tracers = :b)NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 128×1×128 RectilinearGrid{Float64, Periodic, Flat, Bounded} on CPU with 3×0×3 halo
├── timestepper: RungeKutta3TimeStepper
├── advection scheme: UpwindBiased(order=5)
├── tracers: b
├── closure: ScalarDiffusivity{ExplicitTimeDiscretization}(ν=2.0e-5, κ=(b=2.0e-5,))
├── buoyancy: BuoyancyTracer with ĝ = NegativeZDirection()
└── coriolis: NothingWe use hyperbolic tangent functions for the initial conditions and set the maximum Richardson number below the threshold of 1/4. We also add some grid-scale small-amplitude noise to u to kick the instability off:
noise(x, z) = 2e-2 * randn()
shear_flow(x, z) = tanh(z) + noise(x, z)
Ri₀ = 0.1; h = 1/4
stratification(x, z) = h * Ri₀ * tanh(z / h)
set!(model, u=shear_flow, b=stratification)Next create an adaptive-time-step simulation using the model above:
simulation = Simulation(model, Δt=0.1, stop_time=100)
wizard = TimeStepWizard(cfl=0.8, max_Δt=1)
simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(2))Callback of TimeStepWizard(cfl=0.8, max_Δt=1.0, min_Δt=0.0) on IterationInterval(2)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 for of the flow, starting with the RichardsonNumber
Ri = RichardsonNumber(model)KernelFunctionOperation at (Center, Center, Face)
├── grid: 128×1×128 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")We also set-up the QVelocityGradientTensorInvariant, which is usually used for visualizing vortices in the flow:
Q = QVelocityGradientTensorInvariant(model)KernelFunctionOperation at (Center, Center, Center)
├── grid: 128×1×128 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")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.
Let's also keep track of the amount of buoyancy mixing by measuring the buoyancy variance dissipation rate and diffusive term. When volume-integrated, these two quantities should be equal.
∫χᴰ = Integral(TracerVarianceDissipationRate(model, :b))
∫χ = Integral(TracerVarianceDiffusiveTerm(model, :b))Integral of BinaryOperation at (Center, Center, Center) over dims (1, 2, 3)
└── operand: BinaryOperation at (Center, Center, Center)
└── grid: 128×1×128 RectilinearGrid{Float64, Periodic, Flat, Bounded} on CPU with 3×0×3 haloNow we write these quantities, along with b, to a NetCDF:
output_fields = (; Ri, Q, model.tracers.b, ∫χ, ∫χᴰ)
using NCDatasets
filename = "kelvin_helmholtz"
simulation.output_writers[:nc] = NetCDFWriter(model, output_fields,
filename = joinpath(@__DIR__, filename),
schedule = TimeInterval(1),
overwrite_existing = true)NetCDFWriter scheduled on TimeInterval(1 second):
├── filepath: build/generated/kelvin_helmholtz.nc
├── dimensions: time(0), x_faa(128), x_caa(128), z_aaf(129), z_aac(128)
├── 5 outputs: (Q, ∫χᴰ, Ri, b, ∫χ)
└── array type: Array{Float32}
├── file_splitting: NoFileSplitting
└── file size: 24.7 KiBRun the simulation and process results
To run the simulation:
run!(simulation)[ Info: Initializing simulation...
┌ Info: iter = 0, [000.00%] time = 0 seconds, Δt = 58.673 ms, walltime = 29.150 seconds, walltime / timestep = 0 seconds
└ |u⃗|ₘₐₓ = [1.04e+00, 0.00e+00, 3.67e-02] m/s, advective CFL = 0.8, diffusive CFL = 0.00019, νₘₐₓ = 2e-05 m²/s
[ Info: ... simulation initialization complete (12.934 seconds)
[ Info: Executing initial time step...
[ Info: ... initial time step complete (6.144 seconds).
┌ Info: iter = 200, [011.86%] time = 11.858 seconds, Δt = 61.293 ms, walltime = 45.859 seconds, walltime / timestep = 83.543 ms
└ |u⃗|ₘₐₓ = [1.02e+00, 0.00e+00, 3.11e-02] m/s, advective CFL = 0.8, diffusive CFL = 0.0002, νₘₐₓ = 2e-05 m²/s
┌ Info: iter = 400, [023.55%] time = 23.552 seconds, Δt = 61.246 ms, walltime = 48.349 seconds, walltime / timestep = 12.452 ms
└ |u⃗|ₘₐₓ = [1.02e+00, 0.00e+00, 3.90e-02] m/s, advective CFL = 0.8, diffusive CFL = 0.0002, νₘₐₓ = 2e-05 m²/s
┌ Info: iter = 600, [035.30%] time = 35.301 seconds, Δt = 60.102 ms, walltime = 50.790 seconds, walltime / timestep = 12.202 ms
└ |u⃗|ₘₐₓ = [1.03e+00, 0.00e+00, 8.71e-02] m/s, advective CFL = 0.8, diffusive CFL = 0.0002, νₘₐₓ = 2e-05 m²/s
┌ Info: iter = 800, [046.38%] time = 46.378 seconds, Δt = 53.892 ms, walltime = 53.105 seconds, walltime / timestep = 11.578 ms
└ |u⃗|ₘₐₓ = [1.10e+00, 0.00e+00, 2.60e-01] m/s, advective CFL = 0.8, diffusive CFL = 0.00018, νₘₐₓ = 2e-05 m²/s
┌ Info: iter = 1000, [056.15%] time = 56.147 seconds, Δt = 48.842 ms, walltime = 55.307 seconds, walltime / timestep = 11.006 ms
└ |u⃗|ₘₐₓ = [1.17e+00, 0.00e+00, 3.84e-01] m/s, advective CFL = 0.8, diffusive CFL = 0.00016, νₘₐₓ = 2e-05 m²/s
┌ Info: iter = 1200, [065.71%] time = 1.095 minutes, Δt = 50.477 ms, walltime = 57.436 seconds, walltime / timestep = 10.645 ms
└ |u⃗|ₘₐₓ = [1.15e+00, 0.00e+00, 3.53e-01] m/s, advective CFL = 0.8, diffusive CFL = 0.00017, νₘₐₓ = 2e-05 m²/s
┌ Info: iter = 1400, [075.75%] time = 1.262 minutes, Δt = 53.295 ms, walltime = 59.660 seconds, walltime / timestep = 11.120 ms
└ |u⃗|ₘₐₓ = [1.11e+00, 0.00e+00, 2.86e-01] m/s, advective CFL = 0.8, diffusive CFL = 0.00017, νₘₐₓ = 2e-05 m²/s
┌ Info: iter = 1600, [085.95%] time = 1.433 minutes, Δt = 49.961 ms, walltime = 1.031 minutes, walltime / timestep = 11.133 ms
└ |u⃗|ₘₐₓ = [1.15e+00, 0.00e+00, 3.42e-01] m/s, advective CFL = 0.8, diffusive CFL = 0.00016, νₘₐₓ = 2e-05 m²/s
┌ Info: iter = 1800, [095.40%] time = 1.590 minutes, Δt = 49.491 ms, walltime = 1.068 minutes, walltime / timestep = 11.078 ms
└ |u⃗|ₘₐₓ = [1.16e+00, 0.00e+00, 3.55e-01] m/s, advective CFL = 0.8, diffusive CFL = 0.00016, νₘₐₓ = 2e-05 m²/s
[ Info: Simulation is stopping after running for 41.317 seconds.
[ Info: Simulation time 1.667 minutes equals or exceeds stop time 1.667 minutes.Now we'll read the results using Rasters.jl, which works somewhat similarly to Python's Xarray and can speed-up the work the workflow
using Rasters
ds = RasterStack(simulation.output_writers[:nc].filepath)┌ 128×128×129×101×128 RasterStack ┐
├─────────────────────────────────┴────────────────────────────────────── dims ┐
↓ x_caa Sampled{Float32} [-4.9609375f0, …, 4.9609375f0] ForwardOrdered Regular Points,
→ z_aac Sampled{Float32} [-4.9609375f0, …, 4.9609375f0] ForwardOrdered Regular Points,
↗ z_aaf Sampled{Float32} [-5.0f0, …, 5.0f0] ForwardOrdered Regular Points,
⬔ Ti Sampled{Float64} [0.0, …, 100.0] ForwardOrdered Regular Points,
◩ x_faa Sampled{Float32} [-5.0f0, …, 4.921875f0] ForwardOrdered Regular Points
├────────────────────────────────────────────────────────────────────── layers ┤
:Δx_caa eltype: Union{Missing, Float32} dims: x_caa size: 128
:Δx_faa eltype: Union{Missing, Float32} dims: x_faa size: 128
:Δz_aac eltype: Union{Missing, Float32} dims: z_aac size: 128
:Δz_aaf eltype: Union{Missing, Float32} dims: z_aaf size: 129
:Q eltype: Union{Missing, Float32} dims: x_caa, z_aac, Ti size: 128×128×101
:Ri eltype: Union{Missing, Float32} dims: x_caa, z_aaf, Ti size: 128×129×101
:b eltype: Union{Missing, Float32} dims: x_caa, z_aac, Ti size: 128×128×101
:∫χ eltype: Union{Missing, Float32} dims: Ti size: 101
:∫χᴰ eltype: Union{Missing, Float32} dims: Ti size: 101
├──────────────────────────────────────────────────────────────────── metadata ┤
Metadata{Rasters.NCDsource} of Dict{String, Any} with 6 entries:
"interval" => 1.0
"Oceananigans" => "This file was generated using "
"Julia" => "This file was generated using Julia Version 1.10.9…
"output time interval" => "Output was saved every 1 second."
"date" => "This file was generated on 2025-06-03T15:00:00.286…
"schedule" => "TimeInterval"
├────────────────────────────────────────────────────────────────────── raster ┤
missingval: missing
extent: Extent(x_caa = (-4.9609375f0, 4.9609375f0), z_aac = (-4.9609375f0, 4.9609375f0), z_aaf = (-5.0f0, 5.0f0), Ti = (0.0, 100.0), x_faa = (-5.0f0, 4.921875f0))
└──────────────────────────────────────────────────────────────────────────────┘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...);Next we use Observables to lift the values and plot heatmaps and their colorbars
n = Observable(1)
Riₙ = @lift set(ds.Ri[Ti=$n, y_aca=Near(0)], :x_caa => X, :z_aaf => Z)
hm1 = heatmap!(ax1, Riₙ; colormap = :bwr, colorrange = (-1, +1))
Colorbar(fig[3, 1], hm1, vertical=false, height=8)
Qₙ = @lift set(ds.Q[Ti=$n, y_aca=Near(0)], :x_caa => X, :z_aac => Z)
hm2 = heatmap!(ax2, Qₙ; colormap = :inferno, colorrange = (0, 0.2))
Colorbar(fig[3, 2], hm2, vertical=false, height=8)
bₙ = @lift set(ds.b[Ti=$n, y_aca=Near(0)], :x_caa => X, :z_aac => Z)
hm3 = heatmap!(ax3, bₙ; colormap = :balance, colorrange = (-2.5e-2, +2.5e-2))
Colorbar(fig[3, 3], hm3, vertical=false, height=8);┌ Warning: (DimensionalData.Dimensions.Dim{:y_aca},) dims were not found in object.
└ @ DimensionalData.Dimensions ~/.julia/packages/DimensionalData/KlTtO/src/Dimensions/primitives.jl:852
┌ Warning: (DimensionalData.Dimensions.Dim{:y_aca},) dims were not found in object.
└ @ DimensionalData.Dimensions ~/.julia/packages/DimensionalData/KlTtO/src/Dimensions/primitives.jl:852
┌ Warning: (DimensionalData.Dimensions.Dim{:y_aca},) dims were not found in object.
└ @ DimensionalData.Dimensions ~/.julia/packages/DimensionalData/KlTtO/src/Dimensions/primitives.jl:852We now plot the time evolution of our integrated quantities
axb = Axis(fig[4, 1:3]; xlabel="Time", height=100)
times = dims(ds, :Ti)
lines!(axb, Array(times), Array(ds.∫χ), label = "∫χdV")
lines!(axb, Array(times), Array(ds.∫χᴰ), label = "∫χᴰdV", linestyle=:dash)
axislegend(position=:lb, labelsize=14)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!(axb, 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"Similarly to the kinetic energy dissipation rate (see the Two-dimensional turbulence example), TracerVarianceDissipationRate and TracerVarianceDiffusiveTerm are implemented with a energy-conserving formulation, which means that (for NoFlux boundary conditions) their volume-integral should be exactly (up to machine precision) the same.
This page was generated using Literate.jl.