# 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
`ComputedField`

s 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,
advection = UpwindBiasedFifthOrder(),
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 `ComputedField`

s 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 `ComputedField`

s 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
```