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 buoyancy equation or tracers;
• How to create user-defined fields
• How to use differentiation functions
• How to save a computed field

Model setup

We instantiate the model with a simple isotropic diffusivity. We also use a 4-th order advection scheme and Runge-Kutta 3rd order time-stepping scheme.

using Oceananigans, Oceananigans.Advection

model = IncompressibleModel(
        grid = RegularCartesianGrid(size=(128, 128, 1), halo=(2, 2, 2), extent=(2π, 2π, 2π)),
 timestepper = :RungeKutta3,
   advection = CenteredFourthOrder(),
    buoyancy = nothing,
     tracers = nothing,
     closure = IsotropicDiffusivity(ν=1e-4)
)

Setting initial conditions

Our initial condition randomizes u and v. We also ensure that both have zero mean for purely aesthetic reasons.

using Statistics

u₀ = rand(size(model.grid)...)
u₀ .-= mean(u₀)

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

Calculating vorticity

using Oceananigans.Fields, Oceananigans.AbstractOperations

Next we create an object called an ComputedField that calculates vorticity. We'll use this object to calculate vorticity on-line and output it as the simulation progresses.

u, v, w = model.velocities

ω = ComputedField(∂x(v) - ∂y(u))

Now we construct a simulation.

simulation = Simulation(model, Δt=0.1, stop_iteration=2000)

Simulation{IncompressibleModel{CPU, Float64}}
├── Model clock: time = 0 seconds, iteration = 0 
├── Next time step (Float64): 100 ms 
├── Iteration interval: 1
├── 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: Inf years, stop iteration: 2000
├── Diagnostics: OrderedCollections.OrderedDict with no entries
└── 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, (ω = ω,), schedule=IterationInterval(20),
                     prefix = "2d_turbulence_vorticity", force = true)

JLD2OutputWriter scheduled on IterationInterval(20):
├── filepath: ./2d_turbulence_vorticity.jld2
├── 1 outputs: (:ω,)
├── field slicer: FieldSlicer(:, :, :, with_halos=false)
├── array type: Array{Float32}
├── including: [:grid, :coriolis, :buoyancy, :closure]
└── max filesize: Inf YiB

Running the simulation

Finally, we run the simulation.

run!(simulation)

[ Info: Simulation is stopping. Model iteration 2000 has hit or exceeded simulation stop iteration 2000.

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

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

using Oceananigans.Grids

x, y = xnodes(ω)[:], ynodes(ω)[:]

and animate the vorticity.

using Plots

anim = @animate for iteration in iterations

    ω_snapshot = file["timeseries/ω/$iteration"][:, :, 1]

    heatmap(x, y, ω_snapshot',
            xlabel="x", ylabel="y",
            aspectratio=1,
            color=:balance,
            clims=(-0.2, 0.2),
            xlims=(0, model.grid.Lx),
            ylims=(0, model.grid.Ly))
end

This page was generated using Literate.jl.