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.000 s, iteration = 0 
├── Next time step (Float64): 100.000 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.000 s, wall time limit: Inf
├── Stop time: Inf day, 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, (ω = ω,), iteration_interval = 20,
                     prefix = "2d_turbulence_vorticity", force = true)

JLD2OutputWriter{Int64,Nothing,NamedTuple{(:ω,),Tuple{Oceananigans.Fields.ComputedField{Face,Face,Cell,Oceananigans.Fields.FieldStatus{Float64},OffsetArrays.OffsetArray{Float64,3,Array{Float64,3}},RegularCartesianGrid{Float64,Periodic,Periodic,Bounded,OffsetArrays.OffsetArray{Float64,1,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}},Oceananigans.AbstractOperations.BinaryOperation{Face,Face,Cell,typeof(-),Oceananigans.AbstractOperations.Derivative{Face,Face,Cell,typeof(∂xᶠᵃᵃ),OffsetArrays.OffsetArray{Float64,3,Array{Float64,3}},typeof(identity),RegularCartesianGrid{Float64,Periodic,Periodic,Bounded,OffsetArrays.OffsetArray{Float64,1,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}}},Oceananigans.AbstractOperations.Derivative{Face,Face,Cell,typeof(∂yᵃᶠᵃ),OffsetArrays.OffsetArray{Float64,3,Array{Float64,3}},typeof(identity),RegularCartesianGrid{Float64,Periodic,Periodic,Bounded,OffsetArrays.OffsetArray{Float64,1,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}}},typeof(identity),typeof(identity),typeof(identity),RegularCartesianGrid{Float64,Periodic,Periodic,Bounded,OffsetArrays.OffsetArray{Float64,1,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}}}}}},FieldSlicer{Colon,Colon,Colon,Bool},UnionAll,typeof(Oceananigans.OutputWriters.noinit),Array{Symbol,1},Dict{Symbol,Any}}("./2d_turbulence_vorticity.jld2", (ω = Oceananigans.Fields.ComputedField{Face,Face,Cell,Oceananigans.Fields.FieldStatus{Float64},OffsetArrays.OffsetArray{Float64,3,Array{Float64,3}},RegularCartesianGrid{Float64,Periodic,Periodic,Bounded,OffsetArrays.OffsetArray{Float64,1,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}},Oceananigans.AbstractOperations.BinaryOperation{Face,Face,Cell,typeof(-),Oceananigans.AbstractOperations.Derivative{Face,Face,Cell,typeof(∂xᶠᵃᵃ),OffsetArrays.OffsetArray{Float64,3,Array{Float64,3}},typeof(identity),RegularCartesianGrid{Float64,Periodic,Periodic,Bounded,OffsetArrays.OffsetArray{Float64,1,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}}},Oceananigans.AbstractOperations.Derivative{Face,Face,Cell,typeof(∂yᵃᶠᵃ),OffsetArrays.OffsetArray{Float64,3,Array{Float64,3}},typeof(identity),RegularCartesianGrid{Float64,Periodic,Periodic,Bounded,OffsetArrays.OffsetArray{Float64,1,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}}},typeof(identity),typeof(identity),typeof(identity),RegularCartesianGrid{Float64,Periodic,Periodic,Bounded,OffsetArrays.OffsetArray{Float64,1,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}}}}([0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0… ; ]

[0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0… ; ]

[0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0… ; ]

[0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0… ; ]

[0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0… ; ], RegularCartesianGrid{Float64, Periodic, Periodic, Bounded}
                   domain: x ∈ [2.4577456270327724e-17, 6.283185307179588], y ∈ [2.4577456270327724e-17, 6.283185307179588], z ∈ [-6.283185307179586, 0.0]
                 topology: (Periodic, Periodic, Bounded)
  resolution (Nx, Ny, Nz): (128, 128, 1)
   halo size (Hx, Hy, Hz): (2, 2, 2)
grid spacing (Δx, Δy, Δz): (0.049087385212340524, 0.049087385212340524, 6.283185307179586), BinaryOperation at (Face, Face, Cell)
├── grid: RegularCartesianGrid{Float64,Periodic,Periodic,Bounded,OffsetArrays.OffsetArray{Float64,1,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}}
│   ├── size: (128, 128, 1)
│   └── domain: x ∈ [2.4577456270327724e-17, 6.283185307179588], y ∈ [2.4577456270327724e-17, 6.283185307179588], z ∈ [-6.283185307179586, 0.0]
└── tree: 
    - at (Face, Face, Cell) via identity
    ├── ∂xᶠᵃᵃ at (Face, Face, Cell) via identity
    │   └── OffsetArray{Float64, 3, Array{Float64,3}}
    └── ∂yᵃᶠᵃ at (Face, Face, Cell) via identity
        └── OffsetArray{Float64, 3, Array{Float64,3}}, Oceananigans.Fields.FieldStatus{Float64}(0.0)),), 20, nothing, FieldSlicer{Colon,Colon,Colon,Bool}(Colon(), Colon(), Colon(), false), Array{Float32,N} where N, Oceananigans.OutputWriters.noinit, [:grid, :coriolis, :buoyancy, :closure], 0.0, 1, Inf, true, false, Dict{Symbol,Any}())

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.