Two dimensional turbulence example

In this example, we initialize a random velocity field and observe its turbulent decay in a two-dimensional domain. This example demonstrates:

• How to run a model with no tracers and no buoyancy model.

• How to use AbstractOperations to generate output.

Install dependencies

First let's make sure we have all required packages installed.

using Pkg
pkg"add Oceananigans, CairoMakie"

Model setup

We use a grid with 128² points, a fifth-order upwind-biased advection scheme, third-order Runge-Kutta time-stepping (the default), and a small isotropic viscosity. Note that we assign Flat to the z direction.

using Oceananigans
using Random

Random.seed!(404) # for reproducible results

grid = RectilinearGrid(size=(128, 128), extent=(2π, 2π), topology=(Periodic, Periodic, Flat))

model = NonhydrostaticModel(grid;
                            advection = UpwindBiased(order=5),
                            closure = ScalarDiffusivity(ν=1e-5))

NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 128×128×1 RectilinearGrid{Float32, Periodic, Periodic, Flat} on CPU with 3×3×0 halo
├── timestepper: RungeKutta3TimeStepper
├── advection scheme:
│   └── momentum: UpwindBiased(order=5)
├── tracers: ()
├── closure: ScalarDiffusivity{ExplicitTimeDiscretization}(ν=1.0e-5, κ=nothing)
├── buoyancy: Nothing
└── coriolis: Nothing

Random initial conditions

Our initial condition randomizes model.velocities.u and model.velocities.v. We ensure that both have zero mean for aesthetic reasons.

using Statistics

u, v, w = model.velocities

uᵢ = rand(size(u)...)
vᵢ = rand(size(v)...)

uᵢ .-= mean(uᵢ)
vᵢ .-= mean(vᵢ)

set!(model, u=uᵢ, v=vᵢ)

Setting up a simulation

We set up a simulation that stops at 50 time units, with an initial time-step of 0.2, and with adaptive time-stepping and progress printing.

simulation = Simulation(model, Δt=0.2, stop_time=50)

Simulation of NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── Next time step: 200.000 ms
├── run_wall_time: 0 seconds
├── run_wall_time / iteration: NaN days
├── stop_time: 50 seconds
├── stop_iteration: Inf
├── wall_time_limit: Inf
├── minimum_relative_step: 0.0
├── callbacks: OrderedDict with 4 entries:
│   ├── stop_time_exceeded => Callback of stop_time_exceeded on IterationInterval(1)
│   ├── stop_iteration_exceeded => Callback of stop_iteration_exceeded on IterationInterval(1)
│   ├── wall_time_limit_exceeded => Callback of wall_time_limit_exceeded on IterationInterval(1)
│   └── nan_checker => Callback of NaNChecker for u on IterationInterval(100)
└── output_writers: OrderedDict with no entries

A TimeStepWizard helps ensure stable time-stepping with a Courant-Friedrichs-Lewy (CFL) number of 0.7.

conjure_time_step_wizard!(simulation, IterationInterval(10), cfl=0.7, max_change=1.1, max_Δt=0.5)

Logging simulation progress

We set up a callback that logs the simulation iteration and time every 100 iterations.

using Printf

function progress_message(sim)
    u = sim.model.velocities.u
    @info @sprintf("Iteration: %04d, time: %1.3f, Δt: %.2e, max(|u|) = %.1e, wall time: %s",
                   iteration(sim), time(sim), sim.Δt, maximum(abs, u), prettytime(sim.run_wall_time))
    return nothing
end

add_callback!(simulation, progress_message, IterationInterval(100))

Output

We set up an output writer for the simulation that saves vorticity and speed every 0.6 time units.

Computing vorticity and speed

We create two abstract operations that calculate (i) vorticity that measures the rate at which the fluid rotates and is defined as

ω = ∂x(v) - ∂y(u)

BinaryOperation at (Face, Face, Center)
├── grid: 128×128×1 RectilinearGrid{Float32, Periodic, Periodic, Flat} on CPU with 3×3×0 halo
└── tree: 
    - at (Face, Face, Center)
    ├── ∂xᶠᶠᶜ at (Face, Face, Center) via identity
    │   └── 128×128×1 Field{Center, Face, Center} on RectilinearGrid on CPU
    └── ∂yᶠᶠᶜ at (Face, Face, Center) via identity
        └── 128×128×1 Field{Face, Center, Center} on RectilinearGrid on CPU

We also calculate (ii) the speed of the flow,

s = sqrt(u^2 + v^2)

UnaryOperation at (Face, Center, Center)
├── grid: 128×128×1 RectilinearGrid{Float32, Periodic, Periodic, Flat} on CPU with 3×3×0 halo
└── tree: 
    sqrt at (Face, Center, Center) via identity
    └── + at (Face, Center, Center)
        ├── ^ at (Face, Center, Center)
        │   ├── 128×128×1 Field{Face, Center, Center} on RectilinearGrid on CPU
        │   └── 2
        └── ^ at (Center, Face, Center)
            ├── 128×128×1 Field{Center, Face, Center} on RectilinearGrid on CPU
            └── 2

We pass these operations to an output writer, which computes and saves them during the simulation.

filename = "two_dimensional_turbulence.jld2"

simulation.output_writers[:fields] = JLD2Writer(model, (; ω, s); filename,
                                                schedule = TimeInterval(0.6),
                                                overwrite_files = true)

JLD2Writer scheduled on TimeInterval(600 ms):
├── filepath: two_dimensional_turbulence.jld2
├── 2 outputs: (ω, s)
├── array_type: Array{Float32}
├── including: [:coriolis, :buoyancy, :closure]
├── file_splitting: NoFileSplitting
└── file size: 0 bytes (file not yet created)

Running the simulation

Pretty much just

run!(simulation)

[ Info: Initializing simulation...
[ Info: Iteration: 0000, time: 0.000, Δt: 1.00e-01, max(|u|) = 7.3e-01, wall time: 0 seconds
[ Info:     ... simulation initialization complete (5.248 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (509.991 ms).
[ Info: Iteration: 0100, time: 6.900, Δt: 7.07e-02, max(|u|) = 2.9e-01, wall time: 5.941 seconds
[ Info: Iteration: 0200, time: 14.102, Δt: 7.48e-02, max(|u|) = 3.1e-01, wall time: 6.135 seconds
[ Info: Iteration: 0300, time: 21.338, Δt: 8.27e-02, max(|u|) = 2.9e-01, wall time: 6.331 seconds
[ Info: Iteration: 0400, time: 28.357, Δt: 7.90e-02, max(|u|) = 2.8e-01, wall time: 6.528 seconds
[ Info: Iteration: 0500, time: 36.085, Δt: 8.46e-02, max(|u|) = 2.7e-01, wall time: 6.740 seconds
[ Info: Iteration: 0600, time: 44.262, Δt: 7.89e-02, max(|u|) = 2.3e-01, wall time: 6.959 seconds
[ Info: Simulation is stopping after running for 7.127 seconds.
[ Info: Simulation time 50 seconds equals or exceeds stop time 50 seconds.

Visualizing the results

We load the output

ω_timeseries = FieldTimeSeries(filename, "ω")
s_timeseries = FieldTimeSeries(filename, "s")

times = ω_timeseries.times

and animate the vorticity and fluid speed.

using CairoMakie
set_theme!(Theme(fontsize = 20))

fig = Figure(size = (800, 500))

axis_kwargs = (xlabel = "x",
               ylabel = "y",
               limits = ((0, 2π), (0, 2π)),
               aspect = AxisAspect(1))

ax_ω = Axis(fig[2, 1]; title = "Vorticity", axis_kwargs...)
ax_s = Axis(fig[2, 2]; title = "Speed", axis_kwargs...)

We use Makie's Observable to animate the data. To dive into how Observables work we refer to Makie.jl's Documentation.

n = Observable(1)

Observable(1)

Now let's plot the vorticity and speed.

ω = @lift ω_timeseries[$n]
s = @lift s_timeseries[$n]

heatmap!(ax_ω, ω; colormap = :balance, colorrange = (-2, 2))
heatmap!(ax_s, s; colormap = :speed, colorrange = (0, 0.2))

title = @lift "t = $(round(times[$n], digits=2))"
Label(fig[1, 1:2], title, fontsize=24, tellwidth=false)

fig

Finally, we record a movie.

frames = 1:length(times)

@info "Making a neat animation of vorticity and speed..."

record(fig, "two_dimensional_turbulence.mp4", frames, framerate=24) do i
    n[] = i
end

[ Info: Making a neat animation of vorticity and speed...

---

Julia version and environment information

This example was executed with the following version of Julia:

using InteractiveUtils: versioninfo
versioninfo()

Julia Version 1.13.0
Commit d1c37793dd2 (2026-09-09 19:00 UTC)
Build Info:
  Official https://julialang.org release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 128 × AMD EPYC 9374F 32-Core Processor
  WORD_SIZE: 64
  LLVM: libLLVM-20.1.8 (ORCJIT, znver4)
  GC: Built with stock GC
Threads: 1 default, 1 interactive, 1 GC (on 128 virtual cores)
Environment:
  JULIA_LOAD_PATH = @:@v#.#:@stdlib
  JULIA_DEPOT_PATH = /var/lib/buildkite-agent/.julia:/var/lib/buildkite-agent/.julia/juliaup/julia-1.13.0+0.x64.linux.gnu/local/share/julia:/var/lib/buildkite-agent/.julia/juliaup/julia-1.13.0+0.x64.linux.gnu/share/julia
  JULIA_VERSION_ENZYME = 1.13.0
  JULIA_PKG_SERVER_REGISTRY_PREFERENCE = eager
  LD_LIBRARY_PATH = 
  JULIA_MAX_NUM_PRECOMPILE_FILES = 24
  JULIA_VERSION = 1.13.0
  JULIA_CUDA_USE_COMPAT = false
  JULIA_PROJECT = /var/lib/buildkite-agent/Oceananigans.jl-34931/docs/
  JULIA_DEBUG = Literate

These were the top-level packages installed in the environment:

import Pkg
Pkg.status()

Status `~/Oceananigans.jl-34931/docs/Project.toml`
  [79e6a3ab] Adapt v4.7.4
⌃ [052768ef] CUDA v6.2.2
  [13f3f980] CairoMakie v0.15.16
  [8e50ac2c] ConservativeRegridding v0.2.16
⌅ [e30172f5] Documenter v1.17.0
  [daee34ce] DocumenterCitations v1.5.0
  [4710194d] DocumenterVitepress v0.3.8
  [7da242da] Enzyme v0.13.216
  [033835bb] JLD2 v0.6.7
⌅ [63c18a36] KernelAbstractions v0.9.45
  [98b081ad] Literate v2.21.0
  [da04e1cc] MPI v0.20.29
  [85f8d34a] NCDatasets v0.14.15
  [9e8cae18] Oceananigans v0.114.0 `..`
  [429524aa] Optim v2.3.2
  [f27b6e38] Polynomials v4.1.3
⌅ [3c362404] Reactant v0.2.292
  [6038ab10] Rotations v1.7.1
  [d496a93d] SeawaterPolynomials v0.3.11
  [09ab397b] StructArrays v0.7.4
  [bdfc003b] TimesDates v0.3.3
  [0a941bbe] Zarr v0.11.2
  [b77e0a4c] InteractiveUtils v1.11.0
  [37e2e46d] LinearAlgebra v1.13.0
  [44cfe95a] Pkg v1.13.0
Info Packages marked with ⌃ and ⌅ have new versions available. Those with ⌃ may be upgradable, but those with ⌅ are restricted by compatibility constraints from upgrading. To see why use `status --outdated`

---

This page was generated using Literate.jl.