Edit on GitHub

Two dimensional turbulence example

In this example, we initialize a random velocity field and observe its viscous, 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

For this example, we need PyPlot for plotting and Statistics for setting up a random initial condition with zero mean velocity.

using Oceananigans, Oceananigans.AbstractOperations, PyPlot, Statistics

In addition to importing plotting and statistics packages, we import some types from Oceananigans that will aid in the calculation and visualization of voriticty.

using Oceananigans: Face, Cell

Face and Cell represent "locations" on the staggered grid. We instantiate the model with a simple isotropic diffusivity.

model = Model(
        grid = RegularCartesianGrid(size=(128, 128, 1), length=(2π, 2π, 2π)),
    buoyancy = nothing,
     tracers = nothing,
     closure = ConstantIsotropicDiffusivity(ν=1e-3, κ=1e-3)
)
Oceananigans.Model on a CPU architecture (time = 0.000 ns, iteration = 0) 
├── grid: RegularCartesianGrid{Float64,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}
├── tracers: ()
├── closure: ConstantIsotropicDiffusivity{Float64,NamedTuple{(),Tuple{}}}
├── buoyancy: Nothing
├── coriolis: Nothing
├── output writers: OrderedCollections.OrderedDict{Symbol,Oceananigans.AbstractOutputWriter} with no entries
└── diagnostics: OrderedCollections.OrderedDict{Symbol,Oceananigans.AbstractDiagnostic} with no entries

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

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

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

Next we create an object called an Operation that represents a vorticity calculation. We'll use this object to calculate vorticity on-line as the simulation progresses.

u, v, w = model.velocities
(u = Field at (Oceananigans.Face, Oceananigans.Cell, Oceananigans.Cell)
├── data: OffsetArrays.OffsetArray{Float64,3,Array{Float64,3}}
└── grid: RegularCartesianGrid{Float64,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}
    ├── size: (128, 128, 1)
    └── domain: x ∈ [0.0, 6.283185307179586], y ∈ [0.0, 6.283185307179586], z ∈ [-6.283185307179586, 0.0]
, v = Field at (Oceananigans.Cell, Oceananigans.Face, Oceananigans.Cell)
├── data: OffsetArrays.OffsetArray{Float64,3,Array{Float64,3}}
└── grid: RegularCartesianGrid{Float64,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}
    ├── size: (128, 128, 1)
    └── domain: x ∈ [0.0, 6.283185307179586], y ∈ [0.0, 6.283185307179586], z ∈ [-6.283185307179586, 0.0]
, w = Field at (Oceananigans.Cell, Oceananigans.Cell, Oceananigans.Face)
├── data: OffsetArrays.OffsetArray{Float64,3,Array{Float64,3}}
└── grid: RegularCartesianGrid{Float64,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}
    ├── size: (128, 128, 1)
    └── domain: x ∈ [0.0, 6.283185307179586], y ∈ [0.0, 6.283185307179586], z ∈ [-6.283185307179586, 0.0]
)

Create an object that represents the 'operation' required to compute vorticity.

vorticity_operation = ∂x(v) - ∂y(u)
BinaryOperation at (Oceananigans.Face, Oceananigans.Face, Oceananigans.Cell)
├── grid: RegularCartesianGrid{Float64,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}
│   ├── size: (128, 128, 1)
│   └── domain: x ∈ [0.0, 6.283185307179586], y ∈ [0.0, 6.283185307179586], z ∈ [-6.283185307179586, 0.0]
└── tree: 

- at (Oceananigans.Face, Oceananigans.Face, Oceananigans.Cell) via identity
├── ∂x_faa at (Oceananigans.Face, Oceananigans.Face, Oceananigans.Cell) via identity
│   └── OffsetArrays.OffsetArray{Float64,3,Array{Float64,3}}
└── ∂y_afa at (Oceananigans.Face, Oceananigans.Face, Oceananigans.Cell) via identity
    └── OffsetArrays.OffsetArray{Float64,3,Array{Float64,3}}

The instance vorticity_operation is a binary subtraction between two derivative operations acting on OffsetArrays (the underyling representation of u, and v). In order to use vorticity_operation we create a field ω to store the result of the operation, and a Computation object for coordinate the computation of vorticity and storage in ω:

ω = Field(Face, Face, Cell, model.architecture, model.grid)

vorticity_computation = Computation(vorticity_operation, ω)
Oceananigans.AbstractOperations.Computation{Field{Oceananigans.Face,Oceananigans.Face,Oceananigans.Cell,OffsetArrays.OffsetArray{Float64,3,Array{Float64,3}},RegularCartesianGrid{Float64,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}},UnionAll,Oceananigans.AbstractOperations.BinaryOperation{Oceananigans.Face,Oceananigans.Face,Oceananigans.Cell,typeof(-),Oceananigans.AbstractOperations.Derivative{Oceananigans.Face,Oceananigans.Face,Oceananigans.Cell,typeof(Oceananigans.TurbulenceClosures.∂x_faa),OffsetArrays.OffsetArray{Float64,3,Array{Float64,3}},typeof(identity),RegularCartesianGrid{Float64,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}},Oceananigans.AbstractOperations.Derivative{Oceananigans.Face,Oceananigans.Face,Oceananigans.Cell,typeof(Oceananigans.TurbulenceClosures.∂y_afa),OffsetArrays.OffsetArray{Float64,3,Array{Float64,3}},typeof(identity),RegularCartesianGrid{Float64,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}},typeof(identity),typeof(identity),typeof(identity),RegularCartesianGrid{Float64,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}},RegularCartesianGrid{Float64,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}}(BinaryOperation at (Oceananigans.Face, Oceananigans.Face, Oceananigans.Cell)
├── grid: RegularCartesianGrid{Float64,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}
│   ├── size: (128, 128, 1)
│   └── domain: x ∈ [0.0, 6.283185307179586], y ∈ [0.0, 6.283185307179586], z ∈ [-6.283185307179586, 0.0]
└── tree: 

- at (Oceananigans.Face, Oceananigans.Face, Oceananigans.Cell) via identity
├── ∂x_faa at (Oceananigans.Face, Oceananigans.Face, Oceananigans.Cell) via identity
│   └── OffsetArrays.OffsetArray{Float64,3,Array{Float64,3}}
└── ∂y_afa at (Oceananigans.Face, Oceananigans.Face, Oceananigans.Cell) via identity
    └── OffsetArrays.OffsetArray{Float64,3,Array{Float64,3}}, Field at (Oceananigans.Face, Oceananigans.Face, Oceananigans.Cell)
├── data: OffsetArrays.OffsetArray{Float64,3,Array{Float64,3}}
└── grid: RegularCartesianGrid{Float64,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}
    ├── size: (128, 128, 1)
    └── domain: x ∈ [0.0, 6.283185307179586], y ∈ [0.0, 6.283185307179586], z ∈ [-6.283185307179586, 0.0]
, RegularCartesianGrid{Float64}
domain: x ∈ [0.0, 6.283185307179586], y ∈ [0.0, 6.283185307179586], z ∈ [0.0, -6.283185307179586]
  resolution (Nx, Ny, Nz) = (128, 128, 1)
   halo size (Hx, Hy, Hz) = (1, 1, 1)
grid spacing (Δx, Δy, Δz) = (0.04908738521234052, 0.04908738521234052, 6.283185307179586), Array)

We ask for computation of vorticity by writing

compute!(vorticity_computation),

as shown below.

Finally, we run the model.

fig, ax = subplots()

for i = 1:10
    time_step!(model, Nt=100, Δt=1e-1)

    compute!(vorticity_computation)

    cla()
    imshow(interior(ω)[:, :, 1])
    ax.axis("off")
    pause(0.1)
end

We plot out the final vorticity field.

gcf()

This page was generated using Literate.jl.