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 = UpwindBiasedFifthOrder(),
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: Upwind Biased reconstruction 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: ("128×1×128 Field{Face, Center, Center} on RectilinearGrid on CPU", "128×1×128 Field{Center, Face, Center} on RectilinearGrid on CPU", "128×1×129 Field{Center, Center, Face} on RectilinearGrid on CPU", "128×1×128 Field{Center, Center, Center} on RectilinearGrid on CPU", "Tuple{Int64, Int64, Int64}")
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: ("128×1×128 Field{Face, Center, Center} on RectilinearGrid on CPU", "128×1×128 Field{Center, Face, Center} on RectilinearGrid on CPU", "128×1×129 Field{Center, Center, Face} on RectilinearGrid on CPU")
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))
sum! over dims (1, 2, 3) of BinaryOperation at (Center, Center, Center)
└── 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, ∫χ, ∫χᴰ)
filename = "kelvin_helmholtz"
simulation.output_writers[:nc] = NetCDFOutputWriter(model, output_fields,
filename = joinpath(@__DIR__, filename),
schedule = TimeInterval(1),
overwrite_existing = true)
NetCDFOutputWriter scheduled on TimeInterval(1 second):
├── filepath: /home/runner/work/Oceanostics.jl/Oceanostics.jl/docs/build/generated/kelvin_helmholtz.nc
├── dimensions: zC(128), zF(129), xC(128), yF(1), xF(128), yC(1), time(0)
├── 5 outputs: (Q, ∫χᴰ, Ri, b, ∫χ)
└── array type: Array{Float64}
├── file_splitting: NoFileSplitting
└── file size: 20.0 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 = 58.690 ms, walltime = 19.912 seconds, walltime / timestep = 0 seconds
└ |u⃗|ₘₐₓ = [1.05e+00, 0.00e+00, 3.27e-02] m/s, advective CFL = 0.8, diffusive CFL = 0.00019, νₘₐₓ = 2e-05 m²/s
[ Info: ... simulation initialization complete (30.869 minutes)
[ Info: Executing initial time step...
[ Info: ... initial time step complete (5.413 seconds).
┌ Info: iter = 200, [011.79%] time = 11.792 seconds, Δt = 60.934 ms, walltime = 31.330 minutes, walltime / timestep = 9.299 seconds
└ |u⃗|ₘₐₓ = [1.02e+00, 0.00e+00, 3.39e-02] m/s, advective CFL = 0.8, diffusive CFL = 0.0002, νₘₐₓ = 2e-05 m²/s
┌ Info: iter = 400, [023.37%] time = 23.367 seconds, Δt = 61.184 ms, walltime = 31.393 minutes, walltime / timestep = 19.083 ms
└ |u⃗|ₘₐₓ = [1.02e+00, 0.00e+00, 3.07e-02] m/s, advective CFL = 0.8, diffusive CFL = 0.0002, νₘₐₓ = 2e-05 m²/s
┌ Info: iter = 600, [035.12%] time = 35.122 seconds, Δt = 60.825 ms, walltime = 31.456 minutes, walltime / timestep = 18.796 ms
└ |u⃗|ₘₐₓ = [1.02e+00, 0.00e+00, 5.21e-02] m/s, advective CFL = 0.8, diffusive CFL = 0.0002, νₘₐₓ = 2e-05 m²/s
┌ Info: iter = 800, [046.82%] time = 46.825 seconds, Δt = 58.520 ms, walltime = 31.516 minutes, walltime / timestep = 18.081 ms
└ |u⃗|ₘₐₓ = [1.04e+00, 0.00e+00, 1.23e-01] m/s, advective CFL = 0.8, diffusive CFL = 0.00019, νₘₐₓ = 2e-05 m²/s
┌ Info: iter = 1000, [057.46%] time = 57.463 seconds, Δt = 51.185 ms, walltime = 31.576 minutes, walltime / timestep = 17.848 ms
└ |u⃗|ₘₐₓ = [1.14e+00, 0.00e+00, 3.22e-01] m/s, advective CFL = 0.8, diffusive CFL = 0.00017, νₘₐₓ = 2e-05 m²/s
┌ Info: iter = 1200, [067.00%] time = 1.117 minutes, Δt = 48.852 ms, walltime = 31.631 minutes, walltime / timestep = 16.467 ms
└ |u⃗|ₘₐₓ = [1.17e+00, 0.00e+00, 3.85e-01] m/s, advective CFL = 0.8, diffusive CFL = 0.00016, νₘₐₓ = 2e-05 m²/s
┌ Info: iter = 1400, [076.72%] time = 1.279 minutes, Δt = 51.306 ms, walltime = 31.688 minutes, walltime / timestep = 17.121 ms
└ |u⃗|ₘₐₓ = [1.13e+00, 0.00e+00, 3.20e-01] m/s, advective CFL = 0.8, diffusive CFL = 0.00017, νₘₐₓ = 2e-05 m²/s
┌ Info: iter = 1600, [087.00%] time = 1.450 minutes, Δt = 52.915 ms, walltime = 31.744 minutes, walltime / timestep = 17.003 ms
└ |u⃗|ₘₐₓ = [1.11e+00, 0.00e+00, 2.90e-01] m/s, advective CFL = 0.8, diffusive CFL = 0.00017, νₘₐₓ = 2e-05 m²/s
┌ Info: iter = 1800, [096.94%] time = 1.616 minutes, Δt = 49.190 ms, walltime = 31.801 minutes, walltime / timestep = 16.865 ms
└ |u⃗|ₘₐₓ = [1.17e+00, 0.00e+00, 3.77e-01] m/s, advective CFL = 0.8, diffusive CFL = 0.00016, νₘₐₓ = 2e-05 m²/s
[ Info: Simulation is stopping after running for 31.510 minutes.
[ 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)
RasterStack with dimensions:
Dim{:xC} Sampled{Float64} Float64[-4.9609375, -4.8828125, …, 4.8828125, 4.9609375] ForwardOrdered Regular Points,
Dim{:yC} Sampled{Float64} Float64[1.0] ForwardOrdered Regular Points,
Dim{:zC} Sampled{Float64} Float64[-4.9609375, -4.8828125, …, 4.8828125, 4.9609375] ForwardOrdered Regular Points,
Ti Sampled{Float64} Float64[0.0, 1.0, …, 99.0, 100.0] ForwardOrdered Regular Points,
Dim{:zF} Sampled{Float64} Float64[-5.0, -4.921875, …, 4.921875, 5.0] ForwardOrdered Regular Points
and 5 layers:
:Q Float64 dims: Dim{:xC}, Dim{:yC}, Dim{:zC}, Ti (128×1×128×101)
:∫χᴰ Float64 dims: Ti (101)
:Ri Float64 dims: Dim{:xC}, Dim{:yC}, Dim{:zF}, Ti (128×1×129×101)
:b Float64 dims: Dim{:xC}, Dim{:yC}, Dim{:zC}, Ti (128×1×128×101)
:∫χ Float64 dims: Ti (101)
with 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.9.4\…
"output time interval" => "Output was saved every 1 second."
"date" => "This file was generated on 2024-05-23T10:49:34.890…
"schedule" => "TimeInterval"
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...);
ERROR: LoadError: UndefVarError: `IntervalBox` not defined
Stacktrace:
[1] top-level scope
@ ~/.julia/packages/TaylorSeries/2qRvJ/ext/TaylorSeriesIAExt.jl:182
[2] include
@ ./Base.jl:457 [inlined]
[3] include_package_for_output(pkg::Base.PkgId, input::String, depot_path::Vector{String}, dl_load_path::Vector{String}, load_path::Vector{String}, concrete_deps::Vector{Pair{Base.PkgId, UInt128}}, source::String)
@ Base ./loading.jl:2049
[4] top-level scope
@ stdin:3
in expression starting at /home/runner/.julia/packages/TaylorSeries/2qRvJ/ext/TaylorSeriesIAExt.jl:1
in expression starting at stdin:3
┌ Error: Error during loading of extension TaylorSeriesIAExt of TaylorSeries, use `Base.retry_load_extensions()` to retry.
│ exception =
│ 1-element ExceptionStack:
│ Failed to precompile TaylorSeriesIAExt [ed7ef945-33a4-511e-97fe-2b89c7a130ca] to "/home/runner/.julia/compiled/v1.9/TaylorSeriesIAExt/jl_6YN9nT".
│ Stacktrace:
│ [1] error(s::String)
│ @ Base ./error.jl:35
│ [2] compilecache(pkg::Base.PkgId, path::String, internal_stderr::IO, internal_stdout::IO, keep_loaded_modules::Bool)
│ @ Base ./loading.jl:2294
│ [3] compilecache
│ @ ./loading.jl:2167 [inlined]
│ [4] _require(pkg::Base.PkgId, env::Nothing)
│ @ Base ./loading.jl:1805
│ [5] _require_prelocked(uuidkey::Base.PkgId, env::Nothing)
│ @ Base ./loading.jl:1660
│ [6] _require_prelocked(uuidkey::Base.PkgId)
│ @ Base ./loading.jl:1658
│ [7] run_extension_callbacks(extid::Base.ExtensionId)
│ @ Base ./loading.jl:1255
│ [8] run_extension_callbacks(pkgid::Base.PkgId)
│ @ Base ./loading.jl:1290
│ [9] run_package_callbacks(modkey::Base.PkgId)
│ @ Base ./loading.jl:1124
│ [10] _tryrequire_from_serialized(modkey::Base.PkgId, build_id::UInt128)
│ @ Base ./loading.jl:1357
│ [11] _tryrequire_from_serialized(pkg::Base.PkgId, path::String, ocachepath::String)
│ @ Base ./loading.jl:1435
│ [12] _require(pkg::Base.PkgId, env::String)
│ @ Base ./loading.jl:1816
│ [13] _require_prelocked(uuidkey::Base.PkgId, env::String)
│ @ Base ./loading.jl:1660
│ [14] macro expansion
│ @ ./loading.jl:1648 [inlined]
│ [15] macro expansion
│ @ ./lock.jl:267 [inlined]
│ [16] require(into::Module, mod::Symbol)
│ @ Base ./loading.jl:1611
│ [17] eval
│ @ ./boot.jl:370 [inlined]
│ [18] #58
│ @ ~/.julia/packages/Documenter/CJeWX/src/expander_pipeline.jl:754 [inlined]
│ [19] cd(f::Documenter.var"#58#60"{Module, Expr}, dir::String)
│ @ Base.Filesystem ./file.jl:112
│ [20] (::Documenter.var"#57#59"{Documenter.Page, Module, Expr})()
│ @ Documenter ~/.julia/packages/Documenter/CJeWX/src/expander_pipeline.jl:753
│ [21] (::IOCapture.var"#4#7"{DataType, Documenter.var"#57#59"{Documenter.Page, Module, Expr}, IOContext{Base.PipeEndpoint}, IOContext{Base.PipeEndpoint}, IOContext{Base.PipeEndpoint}, IOContext{Base.PipeEndpoint}})()
│ @ IOCapture ~/.julia/packages/IOCapture/Rzdxd/src/IOCapture.jl:161
│ [22] with_logstate(f::Function, logstate::Any)
│ @ Base.CoreLogging ./logging.jl:514
│ [23] with_logger
│ @ ./logging.jl:626 [inlined]
│ [24] capture(f::Documenter.var"#57#59"{Documenter.Page, Module, Expr}; rethrow::Type, color::Bool, passthrough::Bool, capture_buffer::IOBuffer)
│ @ IOCapture ~/.julia/packages/IOCapture/Rzdxd/src/IOCapture.jl:158
│ [25] runner(#unused#::Type{Documenter.Expanders.ExampleBlocks}, node::MarkdownAST.Node{Nothing}, page::Documenter.Page, doc::Documenter.Document)
│ @ Documenter ~/.julia/packages/Documenter/CJeWX/src/expander_pipeline.jl:752
│ [26] dispatch(::Type{Documenter.Expanders.ExpanderPipeline}, ::MarkdownAST.Node{Nothing}, ::Vararg{Any})
│ @ Documenter.Selectors ~/.julia/packages/Documenter/CJeWX/src/utilities/Selectors.jl:170
│ [27] expand(doc::Documenter.Document)
│ @ Documenter ~/.julia/packages/Documenter/CJeWX/src/expander_pipeline.jl:22
│ [28] runner(#unused#::Type{Documenter.Builder.ExpandTemplates}, doc::Documenter.Document)
│ @ Documenter ~/.julia/packages/Documenter/CJeWX/src/builder_pipeline.jl:222
│ [29] dispatch(#unused#::Type{Documenter.Builder.DocumentPipeline}, x::Documenter.Document)
│ @ Documenter.Selectors ~/.julia/packages/Documenter/CJeWX/src/utilities/Selectors.jl:170
│ [30] #86
│ @ ~/.julia/packages/Documenter/CJeWX/src/makedocs.jl:248 [inlined]
│ [31] withenv(::Documenter.var"#86#88"{Documenter.Document}, ::Pair{String, Nothing}, ::Vararg{Pair{String, Nothing}})
│ @ Base ./env.jl:197
│ [32] #85
│ @ ~/.julia/packages/Documenter/CJeWX/src/makedocs.jl:247 [inlined]
│ [33] cd(f::Documenter.var"#85#87"{Documenter.Document}, dir::String)
│ @ Base.Filesystem ./file.jl:112
│ [34] #makedocs#84
│ @ ~/.julia/packages/Documenter/CJeWX/src/makedocs.jl:247 [inlined]
│ [35] top-level scope
│ @ ~/work/Oceanostics.jl/Oceanostics.jl/docs/make.jl:44
│ [36] include(mod::Module, _path::String)
│ @ Base ./Base.jl:457
│ [37] exec_options(opts::Base.JLOptions)
│ @ Base ./client.jl:307
│ [38] _start()
│ @ Base ./client.jl:522
└ @ Base loading.jl:1261
Next we use Observable
s to lift the values and plot heatmaps and their colorbars
n = Observable(1)
Riₙ = @lift set(ds.Ri[Ti=$n, yC=Near(0)], :xC => X, :zF => 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, yC=Near(0)], :xC => X, :zC => 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, yC=Near(0)], :xC => X, :zC => Z)
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 = dims(ds, :Ti)
lines!(axb, Array(times), ds.∫χ, label = "∫χdV")
lines!(axb, Array(times), 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.