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 computed Fields to generate output.

Install dependencies

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

using Pkg
pkg"add Oceananigans, 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 = RectilinearGrid(size=(128, 128), extent=(2π, 2π),
                              topology=(Periodic, Periodic, Flat))

model = NonhydrostaticModel(timestepper = :RungeKutta3,
                              advection = UpwindBiasedFifthOrder(),
                                   grid = grid,
                               buoyancy = nothing,
                                tracers = nothing,
                                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, κ=NamedTuple())
├── 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 "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.

simulation.output_writers[:fields] = JLD2OutputWriter(model, (; ω, s),
                                                      schedule = TimeInterval(2),
                                                      prefix = "two_dimensional_turbulence",
                                                      force = true)

JLD2OutputWriter scheduled on TimeInterval(2 seconds):
├── 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 (10.724 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (33.751 seconds).
[ Info: Iteration: 100, time: 19.199999999999996
[ Info: Iteration: 200, time: 37.800000000000026
[ 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.

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

128×128×1×26 FieldTimeSeries{InMemory} located at (Face, Center, Center) on CPU
├── grid: 128×128×1 RectilinearGrid{Float64, Periodic, Periodic, Flat} on CPU with 3×3×0 halo
├── indices: (1:128, 1:128, 1:1)
└── data: 128×128×1×26 OffsetArray(::Array{Float64, 4}, 1:128, 1:128, 1:1, 1:26) with eltype Float64 with indices 1:128×1:128×1:1×1:26
    └── max=0.854176, min=0.00177511, mean=0.110725

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

xω, yω, zω = nodes(ω_timeseries)
xs, ys, zs = nodes(s_timeseries)

and animate the vorticity and fluid speed.

using Plots

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

anim = @animate for (i, t) in enumerate(ω_timeseries.times)

    @info "Plotting frame $i of $(length(ω_timeseries.times))..."

    ωi = interior(ω_timeseries[i], :, :, 1)
    si = interior(s_timeseries[i], :, :, 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.(ωi', -ω_lim, ω_lim);
                       color = :balance,
                      levels = ω_levels,
                       clims = (-ω_lim, ω_lim),
                      kwargs...)

    s_plot = contourf(xs, ys, clamp.(si', 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.