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"

  Resolving package versions...
No Changes to `/storage7/buildkite-agent/builds/tartarus-mit-edu-6/clima/oceananigans/docs/Project.toml`
No Changes to `/storage7/buildkite-agent/builds/tartarus-mit-edu-6/clima/oceananigans/docs/Manifest.toml`

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.

using Oceananigans, Oceananigans.Advection

grid = RegularCartesianGrid(size=(128, 128, 1), extent=(2π, 2π, 2π))

model = IncompressibleModel(timestepper = :RungeKutta3,
                              advection = UpwindBiasedFifthOrder(),
                                   grid = grid,
                               buoyancy = nothing,
                                tracers = nothing,
                                closure = IsotropicDiffusivity(ν=1e-5)
                           )

IncompressibleModel{CPU, Float64}(time = 0 seconds, iteration = 0) 
├── grid: RegularCartesianGrid{Float64, Periodic, Periodic, Bounded}(Nx=128, Ny=128, Nz=1)
├── tracers: ()
├── closure: IsotropicDiffusivity{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₀ = rand(size(model.grid)...)
u₀ .-= mean(u₀)

set!(model, u=u₀, v=u₀)

Computing vorticity and speed

using Oceananigans.Fields, Oceananigans.AbstractOperations

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, Cell, Cell)
├── data: OffsetArrays.OffsetArray{Float64,3,Array{Float64,3}}, size: (134, 134, 7)
├── grid: RegularCartesianGrid{Float64, Periodic, Periodic, Bounded}(Nx=128, Ny=128, Nz=1)
└── boundary conditions: x=(west=Periodic, east=Periodic), y=(south=Periodic, north=Periodic), z=(bottom=ZeroFlux, top=ZeroFlux), v = Field located at (Cell, Face, Cell)
├── data: OffsetArrays.OffsetArray{Float64,3,Array{Float64,3}}, size: (134, 134, 7)
├── grid: RegularCartesianGrid{Float64, Periodic, Periodic, Bounded}(Nx=128, Ny=128, Nz=1)
└── boundary conditions: x=(west=Periodic, east=Periodic), y=(south=Periodic, north=Periodic), z=(bottom=ZeroFlux, top=ZeroFlux), w = Field located at (Cell, Cell, Face)
├── data: OffsetArrays.OffsetArray{Float64,3,Array{Float64,3}}, size: (134, 134, 8)
├── grid: RegularCartesianGrid{Float64, Periodic, Periodic, Bounded}(Nx=128, Ny=128, Nz=1)
└── boundary conditions: x=(west=Periodic, east=Periodic), y=(south=Periodic, north=Periodic), z=(bottom=NormalFlow, top=NormalFlow))

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, Cell) of BinaryOperation at (Face, Face, Cell)
├── data: OffsetArrays.OffsetArray{Float64,3,Array{Float64,3}}, size: (134, 134, 7)
├── grid: RegularCartesianGrid{Float64, Periodic, Periodic, Bounded}(Nx=128, Ny=128, Nz=1)
├── operand: BinaryOperation at (Face, Face, Cell)
└── 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, Cell, Cell) of UnaryOperation at (Face, Cell, Cell)
├── data: OffsetArrays.OffsetArray{Float64,3,Array{Float64,3}}, size: (134, 134, 7)
├── grid: RegularCartesianGrid{Float64, Periodic, Periodic, Bounded}(Nx=128, Ny=128, Nz=1)
├── operand: UnaryOperation at (Face, Cell, Cell)
└── 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{IncompressibleModel{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: OrderedCollections.OrderedDict with 1 entry:
│   └── nan_checker => NaNChecker
└── Output writers: OrderedCollections.OrderedDict with no entries

Output

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

using Oceananigans.OutputWriters

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: 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 Array{Int64,1}:
   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,

using Oceananigans.Grids

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

This page was generated using Literate.jl.