# 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.
• How to use ComputedFields to generate output.

## Install dependencies

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

using Pkg
pkg"add Oceananigans, JLD2, Plots"

## 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 = RegularRectilinearGrid(size=(128, 128), extent=(2π, 2π),
topology=(Periodic, Periodic, Flat))

model = NonhydrostaticModel(timestepper = :RungeKutta3,
grid = grid,
buoyancy = nothing,
tracers = nothing,
closure = IsotropicDiffusivity(ν=1e-5)
)
NonhydrostaticModel{CPU, Float64}(time = 0 seconds, iteration = 0)
├── grid: RegularRectilinearGrid{Float64, Periodic, Periodic, Flat}(Nx=128, Ny=128, Nz=1)
├── tracers: ()
├── closure: IsotropicDiffusivity{Oceananigans.TurbulenceClosures.ExplicitTimeDiscretization, Float64, NamedTuple{(), Tuple{}}}
├── 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ᵢ)

## 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
(u = Field located at (Face, Center, Center)
├── data: OffsetArrays.OffsetArray{Float64, 3, Array{Float64, 3}}, size: (128, 128, 1)
├── grid: RegularRectilinearGrid{Float64, Periodic, Periodic, Flat}(Nx=128, Ny=128, Nz=1)
└── boundary conditions: west=Periodic, east=Periodic, south=Periodic, north=Periodic, bottom=Nothing, top=Nothing, immersed=ZeroFlux, v = Field located at (Center, Face, Center)
├── data: OffsetArrays.OffsetArray{Float64, 3, Array{Float64, 3}}, size: (128, 128, 1)
├── grid: RegularRectilinearGrid{Float64, Periodic, Periodic, Flat}(Nx=128, Ny=128, Nz=1)
└── boundary conditions: west=Periodic, east=Periodic, south=Periodic, north=Periodic, bottom=Nothing, top=Nothing, immersed=ZeroFlux, w = Field located at (Center, Center, Face)
├── data: OffsetArrays.OffsetArray{Float64, 3, Array{Float64, 3}}, size: (128, 128, 1)
├── grid: RegularRectilinearGrid{Float64, Periodic, Periodic, Flat}(Nx=128, Ny=128, Nz=1)
└── boundary conditions: west=Periodic, east=Periodic, south=Periodic, north=Periodic, bottom=Nothing, top=Nothing, immersed=ZeroFlux)

Next we create two objects called ComputedFields 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)

ω_field = ComputedField(ω)
ComputedField located at (Face, Face, Center) of BinaryOperation at (Face, Face, Center)
├── data: OffsetArrays.OffsetArray{Float64, 3, Array{Float64, 3}}, size: (128, 128, 1)
├── grid: RegularRectilinearGrid{Float64, Periodic, Periodic, Flat}(Nx=128, Ny=128, Nz=1)
├── operand: BinaryOperation at (Face, Face, Center)
└── status: time=0.0

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

$$$s = \sqrt{u^2 + v^2} \, .$$$
s = sqrt(u^2 + v^2)

s_field = ComputedField(s)
ComputedField located at (Face, Center, Center) of UnaryOperation at (Face, Center, Center)
├── data: OffsetArrays.OffsetArray{Float64, 3, Array{Float64, 3}}, size: (128, 128, 1)
├── grid: RegularRectilinearGrid{Float64, Periodic, Periodic, Flat}(Nx=128, Ny=128, Nz=1)
├── operand: UnaryOperation at (Face, Center, Center)
└── status: time=0.0

We'll pass these ComputedFields to an output writer below to calculate them during the simulation. Now we construct a simulation that prints out the iteration and model time as it runs.

progress(sim) = @info "Iteration: $(sim.model.clock.iteration), time:$(round(Int, sim.model.clock.time))"

simulation = Simulation(model, Δt=0.2, stop_time=50, iteration_interval=100, progress=progress)
Simulation{typename(NonhydrostaticModel){typename(CPU), Float64}}
├── Model clock: time = 0 seconds, iteration = 0
├── Next time step (Float64): 200 ms
├── Iteration interval: 100
├── Stop criteria: Any[Oceananigans.Simulations.iteration_limit_exceeded, Oceananigans.Simulations.stop_time_exceeded, Oceananigans.Simulations.wall_time_limit_exceeded]
├── Run time: 0 seconds, wall time limit: Inf
├── Stop time: 50 seconds, stop iteration: Inf
├── Diagnostics: typename(OrderedCollections.OrderedDict) with 1 entry:
│   └── nan_checker => typename(NaNChecker)
└── Output writers: typename(OrderedCollections.OrderedDict) with no entries

## Output

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

simulation.output_writers[:fields] = JLD2OutputWriter(model, (ω=ω_field, s=s_field),
schedule = TimeInterval(2),
prefix = "two_dimensional_turbulence",
force = true)
JLD2OutputWriter scheduled on TimeInterval(2 seconds):
├── filepath: ./two_dimensional_turbulence.jld2
├── 2 outputs: (:ω, :s)
├── field slicer: FieldSlicer(:, :, :, with_halos=false)
├── array type: Array{Float32}
├── including: [:grid, :coriolis, :buoyancy, :closure]
└── max filesize: Inf YiB

## Running the simulation

Pretty much just

run!(simulation)
[ Info: Updating model auxiliary state before the first time step...
[ Info:     ... updated in 387.674 μs.
[ Info: Executing first time step...
[ Info: Iteration: 100, time: 19
[ Info: Iteration: 200, time: 38
[ Info: Iteration: 261, time: 50
[ Info: Simulation is stopping. Model time 50 seconds has hit or exceeded simulation stop time 50 seconds.

## Visualizing the results

We load the output and make a movie.

using JLD2

file = jldopen(simulation.output_writers[:fields].filepath)

iterations = parse.(Int, keys(file["timeseries/t"]))
26-element Vector{Int64}:
0
11
21
32
43
53
63
73
83
94
⋮
181
191
201
211
221
231
241
251
261

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

xω, yω, zω = nodes(ω_field)
xs, ys, zs = nodes(s_field)

and animate the vorticity and fluid speed.

using Plots

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

anim = @animate for (i, iteration) in enumerate(iterations)

@info "Plotting frame $i from iteration$iteration..."

t = file["timeseries/t/$iteration"] ω_snapshot = file["timeseries/ω/$iteration"][:, :, 1]
s_snapshot = file["timeseries/s/\$iteration"][:, :, 1]

ω_lim = 2.0
ω_levels = range(-ω_lim, stop=ω_lim, length=20)

s_lim = 0.2
s_levels = range(0, stop=s_lim, length=20)

kwargs = (xlabel="x", ylabel="y", aspectratio=1, linewidth=0, colorbar=true,
xlims=(0, model.grid.Lx), ylims=(0, model.grid.Ly))

ω_plot = contourf(xω, yω, clamp.(ω_snapshot', -ω_lim, ω_lim);
color = :balance,
levels = ω_levels,
clims = (-ω_lim, ω_lim),
kwargs...)

s_plot = contourf(xs, ys, clamp.(s_snapshot', 0, s_lim);
color = :thermal,
levels = s_levels,
clims = (0., s_lim),
kwargs...)

plot(ω_plot, s_plot, title=["Vorticity" "Speed"], layout=(1, 2), size=(1200, 500))
end