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.