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$ 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: Centered(order=4)
├── tracers: c
├── closure: ScalarDiffusivity{ExplicitTimeDiscretization}(ν=0.0001, κ=(c=0.001,))
├── buoyancy: Nothing
└── coriolis: NothingGrid-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 v32-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.431663181261639Sum 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.44446e-17We 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 KE — 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)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 haloWe 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 NCDatasetsfilename = "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: 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: (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 = 33.967 seconds, advective CFL = 0.25, diffusive CFL = 0.0051
[ Info: ... simulation initialization complete (39.891 seconds)
[ Info: Executing initial time step...
[ Info: ... initial time step complete (2.546 seconds).
[ Info: [002.32%] time = 1.860 seconds, Δt = 12.110 ms, walltime = 1.348 minutes, advective CFL = 0.8, diffusive CFL = 0.02
[ Info: [005.02%] time = 4.013 seconds, Δt = 10.075 ms, walltime = 1.460 minutes, advective CFL = 0.8, diffusive CFL = 0.017
[ Info: [007.93%] time = 6.343 seconds, Δt = 12.454 ms, walltime = 1.580 minutes, advective CFL = 0.8, diffusive CFL = 0.021
[ Info: [011.18%] time = 8.943 seconds, Δt = 14.071 ms, walltime = 1.698 minutes, advective CFL = 0.8, diffusive CFL = 0.023
[ Info: [014.96%] time = 11.971 seconds, Δt = 15.045 ms, walltime = 1.826 minutes, advective CFL = 0.8, diffusive CFL = 0.025
[ Info: [018.52%] time = 14.814 seconds, Δt = 14.067 ms, walltime = 1.959 minutes, advective CFL = 0.8, diffusive CFL = 0.023
[ Info: [022.21%] time = 17.764 seconds, Δt = 15.116 ms, walltime = 2.087 minutes, advective CFL = 0.8, diffusive CFL = 0.025
[ Info: [025.99%] time = 20.793 seconds, Δt = 14.938 ms, walltime = 2.215 minutes, advective CFL = 0.8, diffusive CFL = 0.025
[ Info: [029.52%] time = 23.617 seconds, Δt = 13.491 ms, walltime = 2.344 minutes, advective CFL = 0.8, diffusive CFL = 0.022
[ Info: [032.86%] time = 26.285 seconds, Δt = 13.473 ms, walltime = 2.466 minutes, advective CFL = 0.8, diffusive CFL = 0.022
[ Info: [036.13%] time = 28.906 seconds, Δt = 13.213 ms, walltime = 2.594 minutes, advective CFL = 0.8, diffusive CFL = 0.022
[ Info: [039.53%] time = 31.626 seconds, Δt = 14.253 ms, walltime = 2.716 minutes, advective CFL = 0.8, diffusive CFL = 0.024
[ Info: [043.10%] time = 34.478 seconds, Δt = 14.727 ms, walltime = 2.843 minutes, advective CFL = 0.8, diffusive CFL = 0.024
[ Info: [046.85%] time = 37.478 seconds, Δt = 15.403 ms, walltime = 2.973 minutes, advective CFL = 0.8, diffusive CFL = 0.026
[ Info: [050.58%] time = 40.466 seconds, Δt = 14.727 ms, walltime = 3.102 minutes, advective CFL = 0.8, diffusive CFL = 0.024
[ Info: [054.19%] time = 43.349 seconds, Δt = 14.895 ms, walltime = 3.232 minutes, advective CFL = 0.8, diffusive CFL = 0.025
[ Info: [057.92%] time = 46.336 seconds, Δt = 15.077 ms, walltime = 3.363 minutes, advective CFL = 0.8, diffusive CFL = 0.025
[ Info: [061.66%] time = 49.330 seconds, Δt = 14.308 ms, walltime = 3.491 minutes, advective CFL = 0.8, diffusive CFL = 0.024
[ Info: [065.04%] time = 52.031 seconds, Δt = 13.446 ms, walltime = 3.612 minutes, advective CFL = 0.8, diffusive CFL = 0.022
[ Info: [068.28%] time = 54.625 seconds, Δt = 12.375 ms, walltime = 3.741 minutes, advective CFL = 0.8, diffusive CFL = 0.021
[ Info: [071.18%] time = 56.947 seconds, Δt = 11.736 ms, walltime = 3.853 minutes, advective CFL = 0.8, diffusive CFL = 0.019
[ Info: [074.25%] time = 59.398 seconds, Δt = 13.889 ms, walltime = 3.975 minutes, advective CFL = 0.8, diffusive CFL = 0.023
[ Info: [078.16%] time = 1.042 minutes, Δt = 16.072 ms, walltime = 4.111 minutes, advective CFL = 0.8, diffusive CFL = 0.027
[ Info: [082.02%] time = 1.094 minutes, Δt = 15.622 ms, walltime = 4.245 minutes, advective CFL = 0.8, diffusive CFL = 0.026
[ Info: [086.05%] time = 1.147 minutes, Δt = 17.709 ms, walltime = 4.377 minutes, advective CFL = 0.8, diffusive CFL = 0.029
[ Info: [090.47%] time = 1.206 minutes, Δt = 17.994 ms, walltime = 4.517 minutes, advective CFL = 0.8, diffusive CFL = 0.03
[ Info: [094.50%] time = 1.260 minutes, Δt = 15.147 ms, walltime = 4.645 minutes, advective CFL = 0.8, diffusive CFL = 0.025
[ Info: [098.31%] time = 1.311 minutes, Δt = 16.116 ms, walltime = 4.780 minutes, advective CFL = 0.8, diffusive CFL = 0.027
[ Info: Simulation is stopping after running for 4.308 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")KE_t = FieldTimeSeries(snap_filepath, "KE")c_t = FieldTimeSeries(snap_filepath, "c")ds = NCDataset(snap_filepath)times = ds["time"][:]close(ds)closed DatasetRead 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].filepathds_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 snapshotsidx2 = 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;Source terms at the pair midpoint
ε_pair = @. 0.5 * (∫ε_t[idx1] + ∫ε_t[idx2]);χ_pair = @. 0.5 * (∫χ_t[idx1] + ∫χ_t[idx2]);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.
KE_resid = @. -dKEdt - ε_pairc²_resid = @. -dc²dt - χ_pairPlotting
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_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: the negative tendency -d(∫KE)/dt and -∫ε dV, which sum to 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: the negative tendency -d(∫c²)/dt and -∫χ dV, which sum to 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[] = 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\, dV$ (respectively $-\int \chi\, dV$), 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.