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"

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: Nothing

We 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(TracerVarianceEquation.DissipationRate(model, :b))
∫χ = Integral(TracerVarianceEquation.Diffusion(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 halo

Now 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: 30.8 KiB

Run the simulation and process results

To run the simulation:

run!(simulation)

[ Info: Initializing simulation...
┌ Info: iter =      0,  [000.00%] time = 0 seconds,  Δt = 59.055 ms,  walltime = 51.321 seconds,  walltime / timestep = 0 seconds
└           |u⃗|ₘₐₓ = [1.04e+00,  0.00e+00,  3.64e-02] m/s,  advective CFL = 0.8,  diffusive CFL = 0.00019,  νₘₐₓ = 2e-05 m²/s
[ Info:     ... simulation initialization complete (17.552 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (3.894 seconds).
┌ Info: iter =    200,  [011.79%] time = 11.793 seconds,  Δt = 61.031 ms,  walltime = 1.217 minutes,  walltime / timestep = 108.353 ms
└           |u⃗|ₘₐₓ = [1.02e+00,  0.00e+00,  3.16e-02] m/s,  advective CFL = 0.8,  diffusive CFL = 0.0002,  νₘₐₓ = 2e-05 m²/s
┌ Info: iter =    400,  [023.55%] time = 23.551 seconds,  Δt = 61.237 ms,  walltime = 1.252 minutes,  walltime / timestep = 10.552 ms
└           |u⃗|ₘₐₓ = [1.02e+00,  0.00e+00,  3.79e-02] m/s,  advective CFL = 0.8,  diffusive CFL = 0.0002,  νₘₐₓ = 2e-05 m²/s
┌ Info: iter =    600,  [035.30%] time = 35.298 seconds,  Δt = 59.691 ms,  walltime = 1.286 minutes,  walltime / timestep = 10.211 ms
└           |u⃗|ₘₐₓ = [1.03e+00,  0.00e+00,  8.91e-02] m/s,  advective CFL = 0.8,  diffusive CFL = 0.0002,  νₘₐₓ = 2e-05 m²/s
┌ Info: iter =    800,  [046.43%] time = 46.431 seconds,  Δt = 53.626 ms,  walltime = 1.318 minutes,  walltime / timestep = 9.717 ms
└           |u⃗|ₘₐₓ = [1.10e+00,  0.00e+00,  2.62e-01] m/s,  advective CFL = 0.8,  diffusive CFL = 0.00018,  νₘₐₓ = 2e-05 m²/s
┌ Info: iter =   1000,  [056.25%] time = 56.246 seconds,  Δt = 49.159 ms,  walltime = 1.349 minutes,  walltime / timestep = 9.252 ms
└           |u⃗|ₘₐₓ = [1.16e+00,  0.00e+00,  3.85e-01] m/s,  advective CFL = 0.8,  diffusive CFL = 0.00016,  νₘₐₓ = 2e-05 m²/s
┌ Info: iter =   1200,  [065.81%] time = 1.097 minutes,  Δt = 50.574 ms,  walltime = 1.380 minutes,  walltime / timestep = 9.260 ms
└           |u⃗|ₘₐₓ = [1.15e+00,  0.00e+00,  3.43e-01] m/s,  advective CFL = 0.8,  diffusive CFL = 0.00017,  νₘₐₓ = 2e-05 m²/s
┌ Info: iter =   1400,  [076.00%] time = 1.267 minutes,  Δt = 53.550 ms,  walltime = 1.413 minutes,  walltime / timestep = 9.918 ms
└           |u⃗|ₘₐₓ = [1.10e+00,  0.00e+00,  2.75e-01] m/s,  advective CFL = 0.8,  diffusive CFL = 0.00018,  νₘₐₓ = 2e-05 m²/s
┌ Info: iter =   1600,  [086.30%] time = 1.438 minutes,  Δt = 50.285 ms,  walltime = 1.446 minutes,  walltime / timestep = 9.861 ms
└           |u⃗|ₘₐₓ = [1.15e+00,  0.00e+00,  3.36e-01] m/s,  advective CFL = 0.8,  diffusive CFL = 0.00016,  νₘₐₓ = 2e-05 m²/s
┌ Info: iter =   1800,  [095.89%] time = 1.598 minutes,  Δt = 49.798 ms,  walltime = 1.475 minutes,  walltime / timestep = 8.875 ms
└           |u⃗|ₘₐₓ = [1.16e+00,  0.00e+00,  3.46e-01] m/s,  advective CFL = 0.8,  diffusive CFL = 0.00016,  νₘₐₓ = 2e-05 m²/s
[ Info: Simulation is stopping after running for 40.267 seconds.
[ Info: Simulation time 1.667 minutes equals or exceeds stop time 1.667 minutes.

Now we'll read the results using FieldTimeSeries

filepath = simulation.output_writers[:nc].filepath
Ri_t = FieldTimeSeries(filepath, "Ri")
Q_t  = FieldTimeSeries(filepath, "Q")
b_t  = FieldTimeSeries(filepath, "b")

128×1×128×101 FieldTimeSeries{InMemory} located at (Center, Center, Center) of b at /home/runner/work/Oceanostics.jl/Oceanostics.jl/docs/build/generated/kelvin_helmholtz.nc
├── grid: 128×1×128 RectilinearGrid{Float64, Periodic, Flat, Bounded} on CPU with 3×0×3 halo
├── indices: (1:128, 1:1, 1:128)
├── time_indexing: Linear()
├── backend: InMemory()
├── path: /home/runner/work/Oceanostics.jl/Oceanostics.jl/docs/build/generated/kelvin_helmholtz.nc
├── name: b
└── data: 128×1×128×101 OffsetArray(::Array{Float64, 4}, 1:128, 1:1, 1:128, 1:101) with eltype Float64 with indices 1:128×1:1×1:128×1:101
    └── max=0.0498958, min=-0.0515952, mean=1.93555e-13

Volume-integrated quantities are scalar time series, so we read them directly with NCDatasets:

ds = NCDataset(filepath)
∫χ  = ds["∫χ"][:]
∫χᴰ = ds["∫χᴰ"][:]
close(ds)

closed Dataset

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...);

Precompiling packages...
    533.3 ms  ✓ FilePathsGlobExt (serial)
  1 dependency successfully precompiled in 1 seconds
Precompiling packages...
    428.9 ms  ✓ DistancesChainRulesCoreExt (serial)
  1 dependency successfully precompiled in 0 seconds
Precompiling packages...
   1020.3 ms  ✓ TaylorSeriesIAExt (serial)
  1 dependency successfully precompiled in 1 seconds
Precompiling packages...
   6032.1 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 = (-2.5e-2, +2.5e-2))
Colorbar(fig[3, 3], hm3, vertical=false, height=8);

We now plot the time evolution of our integrated quantities

axb = Axis(fig[4, 1:3]; xlabel="Time", height=100)
times = b_t.times
lines!(axb, times, ∫χ,  label = "∫χdV")
lines!(axb, times, ∫χᴰ, 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 TracerVarianceDiffusion 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.