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
ComputedFields to generate output.
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.
using Oceananigans, Oceananigans.Advection
grid = RegularCartesianGrid(size=(128, 128, 1), extent=(2π, 2π, 2π))
model = IncompressibleModel(timestepper = :RungeKutta3,
advection = UpwindBiasedFifthOrder(),
grid = grid,
buoyancy = nothing,
tracers = nothing,
closure = IsotropicDiffusivity(ν=1e-5)
)IncompressibleModel{CPU, Float64}(time = 0 seconds, iteration = 0)
├── grid: RegularCartesianGrid{Float64, Periodic, Periodic, Bounded}(Nx=128, Ny=128, Nz=1)
├── tracers: ()
├── closure: IsotropicDiffusivity{Float64,NamedTuple{(),Tuple{}}}
├── buoyancy: Nothing
└── coriolis: NothingRandom 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₀ = rand(size(model.grid)...)
u₀ .-= mean(u₀)
set!(model, u=u₀, v=u₀)Computing vorticity and speed
using Oceananigans.Fields, Oceananigans.AbstractOperations
# To make our equations prettier, we unpack `u`, `v`, and `w` from
# the `NamedTuple` model.velocities:
u, v, w = model.velocities(u = Field located at (Face, Cell, Cell)
├── data: OffsetArrays.OffsetArray{Float64,3,Array{Float64,3}}, size: (134, 134, 7)
├── grid: RegularCartesianGrid{Float64, Periodic, Periodic, Bounded}(Nx=128, Ny=128, Nz=1)
└── boundary conditions: x=(west=Periodic, east=Periodic), y=(south=Periodic, north=Periodic), z=(bottom=ZeroFlux, top=ZeroFlux), v = Field located at (Cell, Face, Cell)
├── data: OffsetArrays.OffsetArray{Float64,3,Array{Float64,3}}, size: (134, 134, 7)
├── grid: RegularCartesianGrid{Float64, Periodic, Periodic, Bounded}(Nx=128, Ny=128, Nz=1)
└── boundary conditions: x=(west=Periodic, east=Periodic), y=(south=Periodic, north=Periodic), z=(bottom=ZeroFlux, top=ZeroFlux), w = Field located at (Cell, Cell, Face)
├── data: OffsetArrays.OffsetArray{Float64,3,Array{Float64,3}}, size: (134, 134, 8)
├── grid: RegularCartesianGrid{Float64, Periodic, Periodic, Bounded}(Nx=128, Ny=128, Nz=1)
└── boundary conditions: x=(west=Periodic, east=Periodic), y=(south=Periodic, north=Periodic), z=(bottom=NormalFlow, top=NormalFlow))Next we create two objects called ComputedFields that calculate (i) vorticity, defined as
\[ω = ∂_x v - ∂_y u \, ,\]
ω = ∂x(v) - ∂y(u)
ω_field = ComputedField(ω)ComputedField located at (Face, Face, Cell) of BinaryOperation at (Face, Face, Cell)
├── data: OffsetArrays.OffsetArray{Float64,3,Array{Float64,3}}, size: (134, 134, 7)
├── grid: RegularCartesianGrid{Float64, Periodic, Periodic, Bounded}(Nx=128, Ny=128, Nz=1)
├── operand: BinaryOperation at (Face, Face, Cell)
└── status: time=0.0
"Vorticity" measures the rate at which the fluid rotates. We also calculate (ii) the speed of the flow,
\[s = \sqrt{u^2 + v^2} \, .\]
s = sqrt(u^2 + v^2)
s_field = ComputedField(s)ComputedField located at (Face, Cell, Cell) of UnaryOperation at (Face, Cell, Cell)
├── data: OffsetArrays.OffsetArray{Float64,3,Array{Float64,3}}, size: (134, 134, 7)
├── grid: RegularCartesianGrid{Float64, Periodic, Periodic, Bounded}(Nx=128, Ny=128, Nz=1)
├── operand: UnaryOperation at (Face, Cell, Cell)
└── status: time=0.0
We'll pass these ComputedFields to an output writer below to calculate them during the simulation. Now we construct a simulation that prints out the iteration and model time as it runs.
progress(sim) = @info "Iteration: $(sim.model.clock.iteration), time: $(round(Int, sim.model.clock.time))"
simulation = Simulation(model, Δt=0.2, stop_time=50, iteration_interval=100, progress=progress)Simulation{IncompressibleModel{CPU, Float64}}
├── Model clock: time = 0 seconds, iteration = 0
├── Next time step (Float64): 200 ms
├── Iteration interval: 100
├── 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: 50 seconds, stop iteration: Inf
├── Diagnostics: OrderedCollections.OrderedDict with no entries
└── Output writers: OrderedCollections.OrderedDict with no entriesOutput
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, (ω=ω_field, s=s_field),
schedule = TimeInterval(2),
prefix = "two_dimensional_turbulence",
force = true)JLD2OutputWriter scheduled on TimeInterval(2 seconds):
├── filepath: ./two_dimensional_turbulence.jld2
├── 2 outputs: (:ω, :s)
├── field slicer: FieldSlicer(:, :, :, with_halos=false)
├── array type: Array{Float32}
├── including: [:grid, :coriolis, :buoyancy, :closure]
└── max filesize: Inf YiBRunning the simulation
Pretty much just
run!(simulation)[ Info: Iteration: 100, time: 20 [ Info: Iteration: 200, time: 40 [ Info: Iteration: 300, time: 60 [ Info: Simulation is stopping. Model time 1.000 minutes has hit or exceeded simulation stop time 50 seconds.
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"]))31-element Array{Int64,1}:
0
11
20
31
41
51
60
70
80
90
⋮
220
230
240
250
260
270
280
290
300Construct the $x, y$ grid for plotting purposes,
using Oceananigans.Grids
xω, yω, zω = nodes(ω_field)
xs, ys, zs = nodes(s_field)([-7.512787384681048e-18, 0.04908738521234051, 0.09817477042468102, 0.14726215563702152, 0.19634954084936204, 0.24543692606170256, 0.2945243112740431, 0.3436116964863836, 0.39269908169872414, 0.4417864669110646 … 5.792311455056181, 5.841398840268521, 5.890486225480862, 5.939573610693202, 5.988660995905543, 6.037748381117884, 6.086835766330224, 6.135923151542564, 6.185010536754905, 6.234097921967245], [0.02454369260617027, 0.07363107781851079, 0.12271846303085131, 0.17180584824319184, 0.22089323345553236, 0.2699806186678729, 0.3190680038802134, 0.3681553890925539, 0.4172427743048944, 0.46633015951723494 … 5.816855147662352, 5.865942532874692, 5.915029918087033, 5.964117303299373, 6.013204688511713, 6.062292073724055, 6.111379458936395, 6.160466844148735, 6.209554229361076, 6.258641614573416], [-3.141592653589793])
and animate the vorticity.
using Plots
@info "Making a neat movie of vorticity and speed..."
anim = @animate for (i, iteration) in enumerate(iterations)
@info "Plotting frame $i from iteration $iteration..."
t = file["timeseries/t/$iteration"]
ω_snapshot = file["timeseries/ω/$iteration"][:, :, 1]
s_snapshot = file["timeseries/s/$iteration"][:, :, 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.(ω_snapshot, -ω_lim, ω_lim)';
color = :balance,
levels = ω_levels,
clims = (-ω_lim, ω_lim),
kwargs...)
s_plot = contourf(xs, ys, clamp.(s_snapshot', 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