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,
                            advection = UpwindBiasedFifthOrder(),
                            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
├── 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ᵢ)

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 years
├── 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

Logging simulation progress

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

progress(sim) = @info string("Iteration: ", iteration(sim), ", time: ", time(sim))
simulation.callbacks[:progress] = Callback(progress, IterationInterval(100))
Callback of progress on 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{Float32}
├── including: [:grid, :coriolis, :buoyancy, :closure]
└── max filesize: Inf YiB

Running the simulation

Pretty much just

run!(simulation)
[ Info: Initializing simulation...
[ Info: Iteration: 0, time: 0.0
[ Info:     ... simulation initialization complete (11.222 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (30.689 seconds).
[ Info: Iteration: 100, time: 18.0
[ Info: Iteration: 200, time: 33.000000000000036
[ Info: Iteration: 300, time: 48.00000000000007
[ Info: Simulation is stopping. Model time 50 seconds has hit or exceeded simulation stop time 50 seconds.

Visualizing the results

We load the output.

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

times = ω_timeseries.times
84-element Vector{Float64}:
  0.0
  0.6
  1.2
  1.7999999999999998
  2.4
  3.0
  3.6
  4.2
  4.8
  5.3999999999999995
  ⋮
 45.000000000000064
 45.600000000000065
 46.20000000000007
 46.80000000000007
 47.40000000000007
 48.00000000000007
 48.60000000000007
 49.200000000000074
 49.800000000000075

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))

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

fig = Figure(resolution = (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...)
[ Info: Making a neat movie of vorticity and speed...

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{Int64} with 0 listeners. Value:
1

Now let's plot the vorticity and speed.

ω = @lift interior(ω_timeseries[$n], :, :, 1)
s = @lift interior(s_timeseries[$n], :, :, 1)

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

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

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
    msg = string("Plotting frame ", i, " of ", frames[end])
    print(msg * " \r")
    n[] = i
end
[ Info: Making a neat animation of vorticity and speed...
Plotting frame 1 of 84 
Plotting frame 2 of 84 
Plotting frame 3 of 84 
Plotting frame 4 of 84 
Plotting frame 5 of 84 
Plotting frame 6 of 84 
Plotting frame 7 of 84 
Plotting frame 8 of 84 
Plotting frame 9 of 84 
Plotting frame 10 of 84 
Plotting frame 11 of 84 
Plotting frame 12 of 84 
Plotting frame 13 of 84 
Plotting frame 14 of 84 
Plotting frame 15 of 84 
Plotting frame 16 of 84 
Plotting frame 17 of 84 
Plotting frame 18 of 84 
Plotting frame 19 of 84 
Plotting frame 20 of 84 
Plotting frame 21 of 84 
Plotting frame 22 of 84 
Plotting frame 23 of 84 
Plotting frame 24 of 84 
Plotting frame 25 of 84 
Plotting frame 26 of 84 
Plotting frame 27 of 84 
Plotting frame 28 of 84 
Plotting frame 29 of 84 
Plotting frame 30 of 84 
Plotting frame 31 of 84 
Plotting frame 32 of 84 
Plotting frame 33 of 84 
Plotting frame 34 of 84 
Plotting frame 35 of 84 
Plotting frame 36 of 84 
Plotting frame 37 of 84 
Plotting frame 38 of 84 
Plotting frame 39 of 84 
Plotting frame 40 of 84 
Plotting frame 41 of 84 
Plotting frame 42 of 84 
Plotting frame 43 of 84 
Plotting frame 44 of 84 
Plotting frame 45 of 84 
Plotting frame 46 of 84 
Plotting frame 47 of 84 
Plotting frame 48 of 84 
Plotting frame 49 of 84 
Plotting frame 50 of 84 
Plotting frame 51 of 84 
Plotting frame 52 of 84 
Plotting frame 53 of 84 
Plotting frame 54 of 84 
Plotting frame 55 of 84 
Plotting frame 56 of 84 
Plotting frame 57 of 84 
Plotting frame 58 of 84 
Plotting frame 59 of 84 
Plotting frame 60 of 84 
Plotting frame 61 of 84 
Plotting frame 62 of 84 
Plotting frame 63 of 84 
Plotting frame 64 of 84 
Plotting frame 65 of 84 
Plotting frame 66 of 84 
Plotting frame 67 of 84 
Plotting frame 68 of 84 
Plotting frame 69 of 84 
Plotting frame 70 of 84 
Plotting frame 71 of 84 
Plotting frame 72 of 84 
Plotting frame 73 of 84 
Plotting frame 74 of 84 
Plotting frame 75 of 84 
Plotting frame 76 of 84 
Plotting frame 77 of 84 
Plotting frame 78 of 84 
Plotting frame 79 of 84 
Plotting frame 80 of 84 
Plotting frame 81 of 84 
Plotting frame 82 of 84 
Plotting frame 83 of 84 
Plotting frame 84 of 84


This page was generated using Literate.jl.