# 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 computed Fields 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 instantiate the model with an isotropic diffusivity. We use a grid with 128² points, a fifth-order advection scheme, third-order Runge-Kutta time-stepping, and a small isotropic viscosity. Note that we assign Flat to the z direction.

using Oceananigans

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

model = NonhydrostaticModel(; grid,
timestepper = :RungeKutta3,
closure = ScalarDiffusivity(ν=1e-5))
NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 128×128×1 RectilinearGrid{Float64, Periodic, Periodic, Flat} on CPU with 3×3×0 halo
├── timestepper: RungeKutta3TimeStepper
├── advection scheme: Upwind Biased reconstruction order 5
├── tracers: ()
├── closure: ScalarDiffusivity{ExplicitTimeDiscretization}(ν=1.0e-5)
├── 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.1, 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 ms
├── Elapsed wall time: 0 seconds
├── Wall time per iteration: NaN days
├── Stop time: 50 seconds
├── Stop iteration : Inf
├── Wall time limit: Inf
├── 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
└── Diagnostics: OrderedDict with no entries

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

wizard = TimeStepWizard(cfl=0.7, max_change=1.1, max_Δt=0.5)
simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(10))
Callback of TimeStepWizard(cfl=0.7, max_Δt=0.5, min_Δt=0.0) on IterationInterval(10)

## Logging simulation progress

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

using Printf

function progress_message(sim)
max_abs_u = maximum(abs, sim.model.velocities.u)
walltime = prettytime(sim.run_wall_time)

return @info @sprintf("Iteration: %04d, time: %1.3f, Δt: %.2e, max(|u|) = %.1e, wall time: %s\n",
iteration(sim), time(sim), sim.Δt, max_abs_u, walltime)
end

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

## Output

We set up an output writer for the simulation that saves vorticity and speed every 20 iterations.

### Computing vorticity and speed

To make our equations prettier, we unpack u, v, and w from the NamedTuple model.velocities:

u, v, w = model.velocities
NamedTuple with 3 Fields on 128×128×1 RectilinearGrid{Float64, Periodic, Periodic, Flat} on CPU with 3×3×0 halo:
├── u: 128×128×1 Field{Face, Center, Center} on RectilinearGrid on CPU
├── v: 128×128×1 Field{Center, Face, Center} on RectilinearGrid on CPU
└── w: 128×128×1 Field{Center, Center, Face} on RectilinearGrid on CPU

Next we create two Fields that calculate (i) vorticity that measures the rate at which the fluid rotates and is defined as

$$$ω = ∂_x v - ∂_y u \, ,$$$
ω = ∂x(v) - ∂y(u)
BinaryOperation at (Face, Face, Center)
├── grid: 128×128×1 RectilinearGrid{Float64, 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} \, .$$$
s = sqrt(u^2 + v^2)
UnaryOperation at (Face, Center, Center)
├── grid: 128×128×1 RectilinearGrid{Float64, 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 below to calculate and output them during the simulation.

filename = "two_dimensional_turbulence"

simulation.output_writers[:fields] = JLD2OutputWriter(model, (; ω, s),
schedule = TimeInterval(0.6),
filename = filename * ".jld2",
overwrite_existing = true)
JLD2OutputWriter scheduled on TimeInterval(600 ms):
├── filepath: ./two_dimensional_turbulence.jld2
├── 2 outputs: (ω, s)
├── array type: Array{Float64}
├── including: [:grid, :coriolis, :buoyancy, :closure]
├── file_splitting: NoFileSplitting
└── file size: 27.8 KiB

## Running the simulation

Pretty much just

run!(simulation)
[ Info: Initializing simulation...
[ Info: Iteration: 0000, time: 0.000, Δt: 1.00e-01, max(|u|) = 6.8e-01, wall time: 0 seconds
[ Info:     ... simulation initialization complete (7.844 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (4.352 seconds).
[ Info: Iteration: 0100, time: 6.413, Δt: 6.76e-02, max(|u|) = 3.0e-01, wall time: 13.653 seconds
[ Info: Iteration: 0200, time: 13.758, Δt: 7.51e-02, max(|u|) = 3.1e-01, wall time: 14.895 seconds
[ Info: Iteration: 0300, time: 20.625, Δt: 6.98e-02, max(|u|) = 2.9e-01, wall time: 16.388 seconds
[ Info: Iteration: 0400, time: 27.789, Δt: 9.50e-02, max(|u|) = 2.5e-01, wall time: 17.862 seconds
[ Info: Iteration: 0500, time: 36.253, Δt: 8.09e-02, max(|u|) = 2.6e-01, wall time: 19.177 seconds
[ Info: Iteration: 0600, time: 44.400, Δt: 8.21e-02, max(|u|) = 2.5e-01, wall time: 20.364 seconds
[ Info: Simulation is stopping after running for 21.447 seconds.
[ Info: Simulation time 50.000 seconds equals or exceeds stop time 50 seconds.


## Visualizing the results

ω_timeseries = FieldTimeSeries(filename * ".jld2", "ω")
s_timeseries = FieldTimeSeries(filename * ".jld2", "s")

times = ω_timeseries.times
0.0:0.6:49.8

Construct the $x, y, z$ grid for plotting purposes,

xω, yω, zω = nodes(ω_timeseries)
xs, ys, zs = nodes(s_timeseries)

and animate the vorticity and fluid speed.

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

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 = " * string(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, filename * ".mp4", frames, framerate=24) do i
n[] = i
end
[ Info: Making a neat animation of vorticity and speed...