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"

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
├── 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.1 KiB

Running the simulation

Pretty much just

run!(simulation)

[ Info: Initializing simulation...
[ Info: Iteration: 0000, time: 0.000, Δt: 1.00e-01, max(|u|) = 7.2e-01, wall time: 0 seconds
[ Info:     ... simulation initialization complete (12.373 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (9.186 seconds).
[ Info: Iteration: 0100, time: 6.669, Δt: 7.29e-02, max(|u|) = 3.3e-01, wall time: 24.852 seconds
[ Info: Iteration: 0200, time: 13.615, Δt: 8.29e-02, max(|u|) = 3.0e-01, wall time: 29.956 seconds
[ Info: Iteration: 0300, time: 21.174, Δt: 8.47e-02, max(|u|) = 3.0e-01, wall time: 32.585 seconds
[ Info: Iteration: 0400, time: 28.800, Δt: 8.22e-02, max(|u|) = 3.2e-01, wall time: 36.181 seconds
[ Info: Iteration: 0500, time: 35.854, Δt: 7.37e-02, max(|u|) = 2.7e-01, wall time: 38.831 seconds
[ Info: Iteration: 0600, time: 43.150, Δt: 7.91e-02, max(|u|) = 2.6e-01, wall time: 40.681 seconds
[ Info: Simulation is stopping after running for 42.881 seconds.
[ Info: Simulation time 50 seconds equals or exceeds stop time 50 seconds.

Visualizing the results

We load the output.

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

times = ω_timeseries.times

85-element Vector{Float64}:
  0.0
  0.6
  1.2
  1.7999999999999998
  2.4
  3.0000000000000004
  3.6
  4.200000000000001
  4.8
  5.3999999999999995
  5.999999999999999
  6.599999999999999
  7.199999999999998
  7.799999999999998
  8.399999999999999
  8.999999999999998
  9.599999999999998
 10.199999999999998
 10.799999999999997
 11.399999999999997
 11.999999999999996
 12.599999999999996
 13.199999999999996
 13.799999999999995
 14.399999999999995
 14.999999999999995
 15.599999999999996
 15.6
 16.2
 16.8
 17.400000000000002
 18.000000000000004
 18.60000000000001
 19.2
 19.8
 20.400000000000002
 21.000000000000004
 21.600000000000005
 22.200000000000006
 22.800000000000008
 23.40000000000001
 24.00000000000001
 24.600000000000012
 25.200000000000014
 25.800000000000015
 26.400000000000016
 27.000000000000018
 27.60000000000002
 28.20000000000002
 28.800000000000022
 29.400000000000023
 30.000000000000025
 30.600000000000026
 31.200000000000028
 31.80000000000003
 32.40000000000003
 33.00000000000003
 33.60000000000003
 34.20000000000003
 34.80000000000003
 35.400000000000034
 36.000000000000036
 36.60000000000004
 37.20000000000004
 37.80000000000005
 38.4
 39.0
 39.6
 40.2
 40.800000000000004
 41.400000000000006
 42.00000000000001
 42.60000000000001
 43.20000000000001
 43.80000000000001
 44.40000000000001
 45.000000000000014
 45.600000000000016
 46.20000000000002
 46.80000000000002
 47.40000000000002
 48.00000000000003
 48.6
 49.2
 49.80000000000001

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 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, 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...

This page was generated using Literate.jl.