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 Pkgpkg"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_k$ and $\chi$ alone.

using Oceananigansgrid = 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:
│   ├── momentum: Centered(order=4)
│   └── c: 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, Statisticsu, v, w = model.velocitiesc = model.tracers.cRandom.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 uamp_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.0409281653898634
 -0.12951936644761403
 -0.5538074230529054
 -0.8782796313850649
 -0.6095931282713188
  0.3215956600054153
  1.6761977141306177
 -0.1706308193969656
  0.8717488531798974
  0.8020525821907217
 -0.9533722295902304
 -0.7776191950259158
  1.2655708980723255
 -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=-2.77819e-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 Δtsimulation = 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 Oceanosticsprogress = 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 eₖ — and the dissipation rates εₖ and χ, which we will use to close KE and tracer variance budgets.

using Oceananigans.AbstractOperations: @atspeed     = @at (Center, Center, Center) √(u^2 + v^2)vorticity = ∂x(v) - ∂y(u)eₖ        = 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·qᶜⱼ

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 e_k\, \mathrm{d}V = -\int \varepsilon_k\, \mathrm{d}V,\qquad \frac{d}{dt} \int c^2\, \mathrm{d}V = -\int \chi\, \mathrm{d}V.\]

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_k$ from the current model state and difference $\int e_k\, \mathrm{d}V$ across a time step 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.

∫eₖ = Integral(eₖ)∫c² = Integral(c^2)∫εₖ = Integral(εₖ)∫χ  = Integral(χ)∂ₜ∫eₖ = TimeDerivative(∫eₖ)∂ₜ∫c² = TimeDerivative(∫c²)

TimeDerivative of 1×1×1 Field{Nothing, Nothing, Nothing} reduced over dims = (1, 2, 3) on RectilinearGrid on CPU

We use two NetCDF writers. A visualization writer outputs the 2D snapshot fields and a budget writer only the (cheap) integrated scalars, both on TimeInterval(0.6). Separating the two avoids writing the heavy 2D fields twice per output time.

Each tendency (∂ₜ∫eₖ and ∂ₜ∫c²) is its own entry in the budget writer's outputs, which is how the writer recognizes a TimeDerivative and registers the callback that updates it on the iteration before each output as well as at the output; a TimeDerivative inside a larger operation gets no such callback. The difference therefore spans one model time step. It is a backward difference centered at t - Δt/2, while the source terms are evaluated at t, so the budget closes to first order in the model time step.

using NCDatasetsfilename = "two_dimensional_turbulence"simulation.output_writers[:nc] = NetCDFWriter(model, (; speed, vorticity, eₖ, c),                                              filename = joinpath(@__DIR__, filename),                                              schedule = TimeInterval(0.6),                                              overwrite_files = true)simulation.output_writers[:budget] = NetCDFWriter(model, (; ∂ₜ∫eₖ, ∂ₜ∫c², ∫εₖ, ∫χ),                                                  filename = joinpath(@__DIR__, filename * "_budget"),                                                  schedule = TimeInterval(0.6),                                                  overwrite_files = true)

NetCDFWriter scheduled on TimeInterval(600 ms):
├── filepath: two_dimensional_turbulence_budget.nc
├── dimensions: time(0), x_faa(256), x_caa(256), y_afa(256), y_aca(256)
├── 4 outputs: (∂ₜ∫eₖ, ∂ₜ∫c², ∫εₖ, ∫χ)
├── array_type: Array{Float32}
├── file_splitting: NoFileSplitting
└── file size: (file not yet created)

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 = 41.448 seconds,  advective CFL = 0.25,  diffusive CFL = 0.0051
[ Info:     ... simulation initialization complete (24.053 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (870.328 ms).
[ Info: [002.32%] time = 1.860 seconds,  Δt = 12.110 ms,  walltime = 51.470 seconds,  advective CFL = 0.8,  diffusive CFL = 0.02
[ Info: [005.02%] time = 4.013 seconds,  Δt = 10.075 ms,  walltime = 55.719 seconds,  advective CFL = 0.8,  diffusive CFL = 0.017
[ Info: [007.93%] time = 6.343 seconds,  Δt = 12.454 ms,  walltime = 1.007 minutes,  advective CFL = 0.8,  diffusive CFL = 0.021
[ Info: [011.18%] time = 8.943 seconds,  Δt = 14.071 ms,  walltime = 1.083 minutes,  advective CFL = 0.8,  diffusive CFL = 0.023
[ Info: [014.96%] time = 11.971 seconds,  Δt = 15.045 ms,  walltime = 1.168 minutes,  advective CFL = 0.8,  diffusive CFL = 0.025
[ Info: [018.52%] time = 14.814 seconds,  Δt = 14.067 ms,  walltime = 1.254 minutes,  advective CFL = 0.8,  diffusive CFL = 0.023
[ Info: [022.21%] time = 17.764 seconds,  Δt = 15.116 ms,  walltime = 1.337 minutes,  advective CFL = 0.8,  diffusive CFL = 0.025
[ Info: [025.99%] time = 20.793 seconds,  Δt = 14.938 ms,  walltime = 1.422 minutes,  advective CFL = 0.8,  diffusive CFL = 0.025
[ Info: [029.52%] time = 23.617 seconds,  Δt = 13.491 ms,  walltime = 1.508 minutes,  advective CFL = 0.8,  diffusive CFL = 0.022
[ Info: [032.86%] time = 26.285 seconds,  Δt = 13.473 ms,  walltime = 1.585 minutes,  advective CFL = 0.8,  diffusive CFL = 0.022
[ Info: [036.13%] time = 28.906 seconds,  Δt = 13.213 ms,  walltime = 1.669 minutes,  advective CFL = 0.8,  diffusive CFL = 0.022
[ Info: [039.53%] time = 31.626 seconds,  Δt = 14.253 ms,  walltime = 1.747 minutes,  advective CFL = 0.8,  diffusive CFL = 0.024
[ Info: [043.10%] time = 34.478 seconds,  Δt = 14.727 ms,  walltime = 1.832 minutes,  advective CFL = 0.8,  diffusive CFL = 0.024
[ Info: [046.85%] time = 37.478 seconds,  Δt = 15.403 ms,  walltime = 1.917 minutes,  advective CFL = 0.8,  diffusive CFL = 0.026
[ Info: [050.58%] time = 40.466 seconds,  Δt = 14.727 ms,  walltime = 2.001 minutes,  advective CFL = 0.8,  diffusive CFL = 0.024
[ Info: [054.19%] time = 43.349 seconds,  Δt = 14.895 ms,  walltime = 2.090 minutes,  advective CFL = 0.8,  diffusive CFL = 0.025
[ Info: [057.92%] time = 46.336 seconds,  Δt = 15.077 ms,  walltime = 2.173 minutes,  advective CFL = 0.8,  diffusive CFL = 0.025
[ Info: [061.66%] time = 49.330 seconds,  Δt = 14.308 ms,  walltime = 2.257 minutes,  advective CFL = 0.8,  diffusive CFL = 0.024
[ Info: [065.04%] time = 52.031 seconds,  Δt = 13.446 ms,  walltime = 2.335 minutes,  advective CFL = 0.8,  diffusive CFL = 0.022
[ Info: [068.28%] time = 54.625 seconds,  Δt = 12.375 ms,  walltime = 2.419 minutes,  advective CFL = 0.8,  diffusive CFL = 0.021
[ Info: [071.18%] time = 56.947 seconds,  Δt = 11.736 ms,  walltime = 2.489 minutes,  advective CFL = 0.8,  diffusive CFL = 0.019
[ Info: [074.25%] time = 59.398 seconds,  Δt = 13.889 ms,  walltime = 2.569 minutes,  advective CFL = 0.8,  diffusive CFL = 0.023
[ Info: [078.16%] time = 1.042 minutes,  Δt = 16.072 ms,  walltime = 2.659 minutes,  advective CFL = 0.8,  diffusive CFL = 0.027
[ Info: [082.02%] time = 1.094 minutes,  Δt = 15.622 ms,  walltime = 2.741 minutes,  advective CFL = 0.8,  diffusive CFL = 0.026
[ Info: [086.05%] time = 1.147 minutes,  Δt = 17.709 ms,  walltime = 2.830 minutes,  advective CFL = 0.8,  diffusive CFL = 0.029
[ Info: [090.47%] time = 1.206 minutes,  Δt = 17.994 ms,  walltime = 2.925 minutes,  advective CFL = 0.8,  diffusive CFL = 0.03
[ Info: [094.50%] time = 1.260 minutes,  Δt = 15.147 ms,  walltime = 3.011 minutes,  advective CFL = 0.8,  diffusive CFL = 0.025
[ Info: [098.31%] time = 1.311 minutes,  Δt = 16.116 ms,  walltime = 3.099 minutes,  advective CFL = 0.8,  diffusive CFL = 0.027
[ Info: Simulation is stopping after running for 2.769 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].filepathspeed_t       = FieldTimeSeries(snap_filepath, "speed")vorticity_t   = FieldTimeSeries(snap_filepath, "vorticity")eₖ_t          = FieldTimeSeries(snap_filepath, "eₖ")c_t           = FieldTimeSeries(snap_filepath, "c")ds = NCDataset(snap_filepath)times = ds["time"][:]close(ds)

closed Dataset

Read the budget scalars from the :budget writer. Every record carries both tendencies and both source terms at the same time. The first record has no earlier state to difference against and is written as zero, so the budget starts from the second.

bud_filepath = simulation.output_writers[:budget].filepathds_bud = NCDataset(bud_filepath)nb     = 2:length(ds_bud["time"])t_bud    = ds_bud["time"][nb]deₖdt    = ds_bud["∂ₜ∫eₖ"][nb]dc²dt    = ds_bud["∂ₜ∫c²"][nb]εₖ_bud   = ds_bud["∫εₖ"][nb]χ_bud    = ds_bud["∫χ"][nb]close(ds_bud)

closed Dataset

Budget residuals in sum-to-zero form: the negative tendency plus the source term. Plotting every curve with these signs makes them add up to the residual, which stays near zero.

eₖ_resid = @. -deₖdt - εₖ_budc²_resid = @. -dc²dt - χ_bud

Plotting

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

using CairoMakieset_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_eₖ    = 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]eₖₙ    = @lift eₖ_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_eₖ = heatmap!(ax_eₖ, eₖₙ, colormap = :plasma, colorrange=(0, 0.5))Colorbar(fig[3, 3], hm_eₖ; 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: the negative tendency -d(∫eₖ)/dt and -∫εₖ dV, which sum to the residual.

budget_kwargs = (height = 180, width = 1080)ax_eₖbud = Axis(fig[4, 1:4]; title = "Volume-integrated kinetic energy budget", budget_kwargs...)lines!(ax_eₖbud, t_bud, -deₖdt,   label = "-d(∫eₖ)/dt")lines!(ax_eₖbud, t_bud, -εₖ_bud,  label = "-∫εₖ dV")lines!(ax_eₖbud, t_bud, eₖ_resid, label = "residual", color = :black, linestyle = :dash)axislegend(ax_eₖbud; labelsize = 10, position = :rb)

Makie.Legend()

Volume-integrated c² budget: the negative tendency -d(∫c²)/dt and -∫χ dV, which sum to the residual.

ax_c²bud = Axis(fig[5, 1:4]; title = "Volume-integrated tracer variance budget", xlabel = "Time", budget_kwargs...)lines!(ax_c²bud, t_bud, -dc²dt,  label = "-d(∫c²)/dt")lines!(ax_c²bud, t_bud, -χ_bud,   label = "-∫χ dV")lines!(ax_c²bud, t_bud, 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_eₖbud, 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[] = iend

[ Info: Animating...

The two bottom panels show the volume-integrated KE and $c^2$ budgets. Each plots the negative tendency $-d/dt$ of the integrated quantity alongside $-\int \varepsilon_k\, \mathrm{d}V$ (respectively $-\int \chi\, \mathrm{d}V$), the only source term that survives volume-integration for a periodic incompressible flow with a centered advection scheme. With the tendency negated, the two curves add up to the residual, which stays near zero.

This page was generated using Literate.jl.