Two-dimensional turbulence example

In this example we simulate a 2D flow initialized with random-noise velocities and a passive tracer $c$ with a smooth sine/cosine initial condition. We then use Oceanostics to close the volume-integrated kinetic energy and tracer variance ($c^2$) budgets.

Before starting, make sure you have the required packages installed for this example, which can be done with

using Pkg
pkg"add Oceananigans, Oceanostics, CairoMakie"

Model and simulation setup

We begin by creating a model with an isotropic diffusivity and a fourth-order centered advection scheme on a 256² grid, with one passive tracer c. Using a centered scheme avoids numerical dissipation, so the volume-integrated KE and $c^2$ budgets reduce to purely dissipative balances and we can close them against $\varepsilon$ and $\chi$ alone.

using Oceananigans

grid = RectilinearGrid(size=(256, 256), extent=(2π, 2π), topology=(Periodic, Periodic, Flat))

model = NonhydrostaticModel(grid; timestepper = :RungeKutta3,
                            advection = Centered(order=4),
                            tracers = :c,
                            closure = ScalarDiffusivity(ν=1e-4, κ=1e-3))

NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 256×256×1 RectilinearGrid{Float64, Periodic, Periodic, Flat} on CPU with 3×3×0 halo
├── timestepper: RungeKutta3TimeStepper
├── advection scheme: Centered(order=4)
├── tracers: c
├── closure: ScalarDiffusivity{ExplicitTimeDiscretization}(ν=0.0001, κ=(c=0.001,))
├── buoyancy: Nothing
└── coriolis: Nothing

Grid-scale white noise is not really resolved by the grid, so instead we build a randomized but well-resolved velocity initial condition as a sum of N_blobs Gaussian bumps with random centers and random amplitudes. Each bump is $\sigma_b \approx 10\Delta x$ wide and the periodic copies of each center are summed in so the resulting field is smooth across the periodic boundary. The tracer keeps a smooth sine/cosine pattern.

using Random, Statistics

u, v, w = model.velocities
c = model.tracers.c

Random.seed!(772)
N_blobs = 32
σ_blob  = 10 * minimum_xspacing(grid)
xc      = grid.Lx * rand(N_blobs)
yc      = grid.Ly * rand(N_blobs)
amp_u   = randn(N_blobs) # random Gaussian amplitudes for u
amp_v   = randn(N_blobs) # ... and for v

32-element Vector{Float64}:
 -0.8353820358062429
  1.7747312451395405
 -1.9487671556166866
  0.16949972192372834
 -0.18212397359825486
 -0.3390149109706559
  0.5024285562453012
 -0.7023462100880767
 -0.6899087629358344
 -1.7794241370329904
  ⋮
 -1.5689986344020759
  0.48606157251275944
 -0.043133230163655625
  1.0790526075733167
  1.04042304875431
 -0.26015686321021925
 -0.2621007538287886
 -1.419160151826568
  0.431663181261639

Sum of blobs and their periodic images at (dx, dy) ∈ {-Lx, 0, Lx} × {-Ly, 0, Ly}

blob_sum(x, y, amp) = sum(amp[k] * exp(-((x - xc[k] - dx)^2 + (y - yc[k] - dy)^2) / σ_blob^2)
                          for k  in 1:N_blobs,
                              dx in (-grid.Lx, 0, grid.Lx),
                              dy in (-grid.Ly, 0, grid.Ly))

uᵢ(x, y) = blob_sum(x, y, amp_u)
vᵢ(x, y) = blob_sum(x, y, amp_v)
cᵢ(x, y) = sin(2x) * cos(3y) + cos(x) * sin(2y)

set!(model, u=uᵢ, v=vᵢ, c=cᵢ)

u .-= mean(u)
v .-= mean(v)

256×256×1 Field{Center, Face, Center} on RectilinearGrid on CPU
├── grid: 256×256×1 RectilinearGrid{Float64, Periodic, Periodic, Flat} on CPU with 3×3×0 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: Nothing, top: Nothing, immersed: Nothing
└── data: 262×262×1 OffsetArray(::Array{Float64, 3}, -2:259, -2:259, 1:1) with eltype Float64 with indices -2:259×-2:259×1:1
    └── max=0.865115, min=-1.148, mean=-3.7222e-17

We use this model to create a simulation with a TimeStepWizard to maximize the Δt

u_max = max(maximum(abs, u), maximum(abs, v)) # peak speed magnitude (not signed max)
Δt = 0.2 * minimum_xspacing(grid) / u_max      # Start with a conservative Δt
simulation = Simulation(model; Δt, stop_time=80)

wizard = TimeStepWizard(cfl=0.8, diffusive_cfl=0.8)
simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(5))

Callback of TimeStepWizard(cfl=0.8, max_Δt=Inf, min_Δt=0.0) on IterationInterval(5)

Model diagnostics

Up until now we have only used Oceananigans, but we can make use of Oceanostics for the first diagnostic we'll set-up: a progress messenger. Here we use a BasicMessenger, which, as the name suggests, displays only basic information about the simulation

using Oceanostics

progress = ProgressMessengers.BasicMessenger()
simulation.callbacks[:progress] = Callback(progress, IterationInterval(200))

Callback of Oceanostics.ProgressMessengers.BasicMessenger{Oceanostics.ProgressMessengers.AbstractProgressMessenger} on IterationInterval(200)

We define the visualization fields — speed, vorticity, kinetic energy KE — and the dissipation rates ε and χ, which we will use to close KE and tracer variance budgets.

using Oceananigans.AbstractOperations: @at

speed     = @at (Center, Center, Center) √(u^2 + v^2)
vorticity = ∂x(v) - ∂y(u)
KE        = KineticEnergyEquation.KineticEnergy(model)
ε         = KineticEnergyEquation.DissipationRate(model)
χ         = TracerVarianceEquation.TracerVarianceDissipationRate(model, :c)

TracerVarianceDissipationRate KernelFunctionOperation at (Center, Center, Center)
├── grid: 256×256×1 RectilinearGrid{Float64, Periodic, Periodic, Flat} on CPU with 3×3×0 halo
├── kernel_function: tracer_variance_dissipation_rate_ccc (generic function with 1 method)
└── arguments: ("ScalarDiffusivity", "Nothing", "Val", "Field", "Clock", "NamedTuple", "Nothing")
└── computes: tracer variance dissipation rate  χ = 2 ∂ⱼc·Fⱼ

Note that KineticEnergyEquation.DissipationRate (ε) and TracerVarianceEquation.TracerVarianceDissipationRate (χ) –- which can also be called as KineticEnergyDissipationRate and TracerVarianceDissipationRate –- are implemented using the same kernels as Oceananigans (and therefore use the same interpolations and discretizations).

To close the budgets we also define the relevant volume integrals as scalar outputs. For a 2D periodic domain with no forcing or buoyancy, advection and pressure-redistribution terms volume-integrate to zero due to incompressibility, so the volume-integrated KE and $c^2$ evolution equations reduce to

\[\frac{d}{dt} \int \mathrm{KE}\, dV = -\int \varepsilon\, dV,\qquad \frac{d}{dt} \int c^2\, dV = -\int \chi\, dV.\]

A caveat: a discretized version of the continuum KE equation (such as the one above) is not guaranteed to exactly conserve energy at the discrete level. To get strict discrete conservation of energy one would have to derive a discrete KE equation directly from the discrete momentum equations — using both the current and previous time-step velocities. We are not doing that here: we compute $\varepsilon$ from the current model state and finite-difference snapshots of $\int \mathrm{KE}\, dV$ independently. The two relations are consistent in the continuum limit but only approximately at the discrete level for a well-resolved flow, so we expect the KE budget to close only approximately.

∫KE = Integral(KE)
∫c² = Integral(c^2)
∫ε  = Integral(ε)
∫χ  = Integral(χ)

Integral of BinaryOperation at (Center, Center, Center) over dims (1, 2, 3)
└── operand: BinaryOperation at (Center, Center, Center)
    └── grid: 256×256×1 RectilinearGrid{Float64, Periodic, Periodic, Flat} on CPU with 3×3×0 halo

We use two NetCDF writers. A visualization writer outputs the 2D snapshot fields on a plain TimeInterval(0.6). A budget writer outputs only the (cheap) integrated scalars on ConsecutiveIterations(TimeInterval(0.6)) — i.e. a second sample one model step after each scheduled time — which lets us finite-difference the integrated quantities across that single step to estimate $d/dt$ without time-integration accumulators. Separating the two avoids writing the heavy 2D fields twice per output time.

using NCDatasets
filename = "two_dimensional_turbulence"

simulation.output_writers[:nc] = NetCDFWriter(model, (; speed, vorticity, KE, c),
                                              filename = joinpath(@__DIR__, filename),
                                              schedule = TimeInterval(0.6),
                                              overwrite_existing = true)

simulation.output_writers[:budget] = NetCDFWriter(model, (; ∫KE, ∫c², ∫ε, ∫χ),
                                                  filename = joinpath(@__DIR__, filename * "_budget"),
                                                  schedule = ConsecutiveIterations(TimeInterval(0.6)),
                                                  overwrite_existing = true)

NetCDFWriter scheduled on ConsecutiveIterations(TimeInterval(600 ms), 1):
├── filepath: build/generated/two_dimensional_turbulence_budget.nc
├── dimensions: time(0), y_afa(256), x_faa(256), x_caa(256), y_aca(256)
├── 4 outputs: (∫ε, ∫c², ∫KE, ∫χ)
├── array_type: Array{Float32}
├── file_splitting: NoFileSplitting
└── file size: 34.4 KiB

Run the simulation and process results

To run the simulation:

run!(simulation)

[ Info: Initializing simulation...
[ Info: [000.00%] time = 0 seconds,  Δt = 3.069 ms,  walltime = 22.052 seconds,  advective CFL = 0.25,  diffusive CFL = 0.0051
[ Info:     ... simulation initialization complete (3.609 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (283.796 ms).
[ Info: [002.32%] time = 1.860 seconds,  Δt = 12.110 ms,  walltime = 31.459 seconds,  advective CFL = 0.8,  diffusive CFL = 0.02
[ Info: [005.02%] time = 4.013 seconds,  Δt = 10.075 ms,  walltime = 37.265 seconds,  advective CFL = 0.8,  diffusive CFL = 0.017
[ Info: [007.93%] time = 6.343 seconds,  Δt = 12.454 ms,  walltime = 43.826 seconds,  advective CFL = 0.8,  diffusive CFL = 0.021
[ Info: [011.18%] time = 8.943 seconds,  Δt = 14.071 ms,  walltime = 50.669 seconds,  advective CFL = 0.8,  diffusive CFL = 0.023
[ Info: [014.96%] time = 11.971 seconds,  Δt = 15.045 ms,  walltime = 57.971 seconds,  advective CFL = 0.8,  diffusive CFL = 0.025
[ Info: [018.52%] time = 14.814 seconds,  Δt = 14.067 ms,  walltime = 1.087 minutes,  advective CFL = 0.8,  diffusive CFL = 0.023
[ Info: [022.21%] time = 17.764 seconds,  Δt = 15.116 ms,  walltime = 1.209 minutes,  advective CFL = 0.8,  diffusive CFL = 0.025
[ Info: [025.99%] time = 20.793 seconds,  Δt = 14.938 ms,  walltime = 1.332 minutes,  advective CFL = 0.8,  diffusive CFL = 0.025
[ Info: [029.52%] time = 23.617 seconds,  Δt = 13.491 ms,  walltime = 1.455 minutes,  advective CFL = 0.8,  diffusive CFL = 0.022
[ Info: [032.86%] time = 26.285 seconds,  Δt = 13.473 ms,  walltime = 1.569 minutes,  advective CFL = 0.8,  diffusive CFL = 0.022
[ Info: [036.13%] time = 28.906 seconds,  Δt = 13.213 ms,  walltime = 1.691 minutes,  advective CFL = 0.8,  diffusive CFL = 0.022
[ Info: [039.53%] time = 31.626 seconds,  Δt = 14.253 ms,  walltime = 1.802 minutes,  advective CFL = 0.8,  diffusive CFL = 0.024
[ Info: [043.10%] time = 34.478 seconds,  Δt = 14.727 ms,  walltime = 1.922 minutes,  advective CFL = 0.8,  diffusive CFL = 0.024
[ Info: [046.85%] time = 37.478 seconds,  Δt = 15.403 ms,  walltime = 2.044 minutes,  advective CFL = 0.8,  diffusive CFL = 0.026
[ Info: [050.58%] time = 40.466 seconds,  Δt = 14.727 ms,  walltime = 2.166 minutes,  advective CFL = 0.8,  diffusive CFL = 0.024
[ Info: [054.19%] time = 43.349 seconds,  Δt = 14.895 ms,  walltime = 2.287 minutes,  advective CFL = 0.8,  diffusive CFL = 0.025
[ Info: [057.92%] time = 46.336 seconds,  Δt = 15.077 ms,  walltime = 2.408 minutes,  advective CFL = 0.8,  diffusive CFL = 0.025
[ Info: [061.66%] time = 49.330 seconds,  Δt = 14.308 ms,  walltime = 2.525 minutes,  advective CFL = 0.8,  diffusive CFL = 0.024
[ Info: [065.04%] time = 52.031 seconds,  Δt = 13.446 ms,  walltime = 2.637 minutes,  advective CFL = 0.8,  diffusive CFL = 0.022
[ Info: [068.28%] time = 54.625 seconds,  Δt = 12.375 ms,  walltime = 2.756 minutes,  advective CFL = 0.8,  diffusive CFL = 0.021
[ Info: [071.18%] time = 56.947 seconds,  Δt = 11.736 ms,  walltime = 2.859 minutes,  advective CFL = 0.8,  diffusive CFL = 0.019
[ Info: [074.25%] time = 59.398 seconds,  Δt = 13.889 ms,  walltime = 2.970 minutes,  advective CFL = 0.8,  diffusive CFL = 0.023
[ Info: [078.16%] time = 1.042 minutes,  Δt = 16.072 ms,  walltime = 3.100 minutes,  advective CFL = 0.8,  diffusive CFL = 0.027
[ Info: [082.02%] time = 1.094 minutes,  Δt = 15.622 ms,  walltime = 3.214 minutes,  advective CFL = 0.8,  diffusive CFL = 0.026
[ Info: [086.05%] time = 1.147 minutes,  Δt = 17.709 ms,  walltime = 3.327 minutes,  advective CFL = 0.8,  diffusive CFL = 0.029
[ Info: [090.47%] time = 1.206 minutes,  Δt = 17.994 ms,  walltime = 3.448 minutes,  advective CFL = 0.8,  diffusive CFL = 0.03
[ Info: [094.50%] time = 1.260 minutes,  Δt = 15.147 ms,  walltime = 3.561 minutes,  advective CFL = 0.8,  diffusive CFL = 0.025
[ Info: [098.31%] time = 1.311 minutes,  Δt = 16.116 ms,  walltime = 3.683 minutes,  advective CFL = 0.8,  diffusive CFL = 0.027
[ Info: Simulation is stopping after running for 3.366 minutes.
[ Info: Simulation time 1.333 minutes equals or exceeds stop time 1.333 minutes.

Read visualization snapshots from the :nc writer.

snap_filepath = simulation.output_writers[:nc].filepath
speed_t       = FieldTimeSeries(snap_filepath, "speed")
vorticity_t   = FieldTimeSeries(snap_filepath, "vorticity")
KE_t          = FieldTimeSeries(snap_filepath, "KE")
c_t           = FieldTimeSeries(snap_filepath, "c")

ds = NCDataset(snap_filepath)
times = ds["time"][:]
close(ds)

closed Dataset

Read the integrated-quantity scalars from the :budget writer. These come in consecutive- iteration pairs: (t₀, t₀+Δt_model, t₀+0.6, t₀+0.6+Δt_model, …). Pair k has indices (2k-1, 2k); we obtain $d/dt$ from a one-step finite difference inside each pair.

bud_filepath = simulation.output_writers[:budget].filepath
ds_bud = NCDataset(bud_filepath)
times_bud = ds_bud["time"][:]
∫KE_t = ds_bud["∫KE"][:]
∫c²_t = ds_bud["∫c²"][:]
∫ε_t  = ds_bud["∫ε"][:]
∫χ_t  = ds_bud["∫χ"][:]
close(ds_bud)

idx1     = 1:2:length(times_bud) - 1   # primary snapshots
idx2     = 2:2:length(times_bud)       # consecutive-iteration snapshots
Δt_pair  = times_bud[idx2] .- times_bud[idx1]
t_pair   = @. 0.5 * (times_bud[idx1] + times_bud[idx2])

dKEdt    = (∫KE_t[idx2] .- ∫KE_t[idx1]) ./ Δt_pair
dc²dt    = (∫c²_t[idx2] .- ∫c²_t[idx1]) ./ Δt_pair

134-element Vector{Float64}:
 -0.3548989047406253
 -0.4617221626657696
 -0.7438569497464611
 -1.0578007294339566
 -1.2609977148921634
 -1.3585533075002856
 -1.4201250020239262
 -1.4699955056855198
 -1.4516182558510475
 -1.3581433652232064
  ⋮
 -0.002569251363929729
 -0.0025276220686924426
 -0.002476650844018774
 -0.002432534401397486
 -0.002394865684380938
 -0.0023563589462259906
 -0.0023165350455735477
 -0.0022798143275940346
 -0.0022476122575162028

Source terms at the pair midpoint

ε_pair   = @. 0.5 * (∫ε_t[idx1] + ∫ε_t[idx2])
χ_pair   = @. 0.5 * (∫χ_t[idx1] + ∫χ_t[idx2])

KE_resid = @. dKEdt - (-ε_pair)
c²_resid = @. dc²dt - (-χ_pair)

Plotting

We use Makie.jl, which has recipes for Oceananigans Fields.

using CairoMakie

set_theme!(Theme(fontsize = 20))
fig = Figure()

axis_kwargs = (aspect = DataAspect(),
               height = 250, width = 250,
               xticksvisible = false, yticksvisible = false,
               xticklabelsvisible = false, yticklabelsvisible = false)

ax_speed = Axis(fig[2, 1]; title = "Speed",          axis_kwargs...)
ax_ω     = Axis(fig[2, 2]; title = "Vorticity",      axis_kwargs...)
ax_KE    = Axis(fig[2, 3]; title = "Kinetic energy", axis_kwargs...)
ax_c     = Axis(fig[2, 4]; title = "Tracer c",       axis_kwargs...)

Makie.Axis with 0 plots:

Each frame is one visualization snapshot.

n = Observable(1)

speedₙ = @lift speed_t[$n]
ωₙ     = @lift vorticity_t[$n]
KEₙ    = @lift KE_t[$n]
cₙ     = @lift c_t[$n]

hm_speed = heatmap!(ax_speed, speedₙ, colormap = :magma, colorrange=(0, 1.5))
Colorbar(fig[3, 1], hm_speed; vertical=false, height=8, ticklabelsize=12)

hm_ω = heatmap!(ax_ω, ωₙ, colormap = :balance, colorrange=(-10, 10))
Colorbar(fig[3, 2], hm_ω; vertical=false, height=8, ticklabelsize=12)

hm_KE = heatmap!(ax_KE, KEₙ, colormap = :plasma, colorrange=(0, 0.5))
Colorbar(fig[3, 3], hm_KE; vertical=false, height=8, ticklabelsize=12)

hm_c = heatmap!(ax_c, cₙ, colormap = :balance, colorrange=(-1.5, 1.5))
Colorbar(fig[3, 4], hm_c; vertical=false, height=8, ticklabelsize=12)

Makie.Colorbar()

Volume-integrated KE budget — d(∫KE)/dt against -∫ε dV, with the residual.

budget_kwargs = (height = 180, width = 1080)

ax_KEbud = Axis(fig[4, 1:4]; title = "Volume-integrated KE budget", budget_kwargs...)
lines!(ax_KEbud, t_pair, dKEdt,   label = "d(∫KE)/dt")
lines!(ax_KEbud, t_pair, -ε_pair, label = "-∫ε dV")
lines!(ax_KEbud, t_pair, KE_resid, label = "residual", color = :black, linestyle = :dash)
axislegend(ax_KEbud; labelsize = 10, position = :rb)

Makie.Legend()

Volume-integrated c² budget — d(∫c²)/dt against -∫χ dV, with the residual.

ax_c²bud = Axis(fig[5, 1:4]; title = "Volume-integrated ∫c² budget", xlabel = "Time", budget_kwargs...)
lines!(ax_c²bud, t_pair, dc²dt,   label = "d(∫c²)/dt")
lines!(ax_c²bud, t_pair, -χ_pair, label = "-∫χ dV")
lines!(ax_c²bud, t_pair, c²_resid, label = "residual", color = :black, linestyle = :dash)
axislegend(ax_c²bud; labelsize = 10, position = :rb)

Makie.Legend()

Time marker on both budget panels (using the snapshot time shown in the heatmaps)

tₙ = @lift times[$n]
vlines!(ax_KEbud, tₙ, color = :black, linestyle = :dot)
vlines!(ax_c²bud, tₙ, color = :black, linestyle = :dot)

title = @lift "Time = " * string(round(times[$n], digits=2))
Label(fig[1, 1:4], title, fontsize=24, tellwidth=false);

Adjust the total figure size based on our panels and record a movie.

resize_to_layout!(fig)
@info "Animating..."
record(fig, filename * ".mp4", 1:length(times), framerate=24) do i
    n[] = i
end

[ Info: Animating...

The two bottom panels show the volume-integrated KE and $c^2$ budgets: $d/dt$ of the integrated quantity is compared against $-\int \varepsilon\, dV$ (respectively $-\int \chi\, dV$), the only term that survives volume-integration for a periodic incompressible flow with a centered advection scheme. The residual shows the gap between them.

This page was generated using Literate.jl.