Two-dimensional turbulence example

In this example (based on the homonymous Oceananigans one) we simulate a 2D flow initialized with random noise and observe the flow evolve.

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 fifth-order advection on a 128² grid.

using Oceananigans

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

model = NonhydrostaticModel(grid; timestepper = :RungeKutta3,
                            advection = UpwindBiased(order=5),
                            closure = ScalarDiffusivity(ν=1e-5))

NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 128×128×1 RectilinearGrid{Float64, Periodic, Periodic, Flat} on CPU with 3×3×0 halo
├── timestepper: RungeKutta3TimeStepper
├── advection scheme: UpwindBiased(order=5)
├── tracers: ()
├── closure: ScalarDiffusivity{ExplicitTimeDiscretization}(ν=1.0e-5)
├── buoyancy: Nothing
└── coriolis: Nothing

Let's give the model zero-mean grid-scale white noise as the initial condition

using Statistics

u, v, w = model.velocities

noise(x, y) = rand()
set!(model, u=noise, v=noise)

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

128×128×1 Field{Center, Face, Center} on RectilinearGrid on CPU
├── grid: 128×128×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: 134×134×1 OffsetArray(::Array{Float64, 3}, -2:131, -2:131, 1:1) with eltype Float64 with indices -2:131×-2:131×1:1
    └── max=0.641862, min=-0.794923, mean=4.0705e-17

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

simulation = Simulation(model, Δt=0.2, stop_time=50)

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 basic information about the simulation

using Oceanostics

progress = ProgressMessengers.BasicMessenger()

simulation.callbacks[:progress] = Callback(progress, IterationInterval(100))

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

Using Oceanostics we can easily calculate two important diagnostics, the kinetic energy KE and its dissipation rate ε

KE = KineticEnergyEquation.KineticEnergy(model)
ε = KineticEnergyEquation.DissipationRate(model)

KernelFunctionOperation at (Center, Center, Center)
├── grid: 128×128×1 RectilinearGrid{Float64, Periodic, Periodic, Flat} on CPU with 3×3×0 halo
├── kernel_function: viscous_dissipation_rate_ccc (generic function with 1 method)
└── arguments: ("Nothing", "NamedTuple", "NamedTuple")

And we can define their volume-integrals

∫KE = Integral(KE)
∫ε = Integral(ε)

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

We also create another integrated quantity that appears in the TKE evolution equation: the KineticEnergyStress, which in our case is

\[\varepsilon^D = u_i \partial_j \tau_{ij}\]

where $\tau_{ij}$ is the diffusive flux of $i$ momentum in the $j$-th direction.

∫εᴰ = Integral(KineticEnergyEquation.Stress(model))

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

The idea in calculating this term is that, in integrated form, all transport contributions in it should equal zero and ∫εᴰ should equal ∫ε.

We output the previous quantities to a NetCDF file

output_fields = (; KE, ε, ∫KE, ∫ε, ∫εᴰ)

using NCDatasets
filename = "two_dimensional_turbulence"
simulation.output_writers[:nc] = NetCDFWriter(model, output_fields,
                                              filename = joinpath(@__DIR__, filename),
                                              schedule = TimeInterval(0.6),
                                              overwrite_existing = true)

NetCDFWriter scheduled on TimeInterval(600 ms):
├── filepath: build/generated/two_dimensional_turbulence.nc
├── dimensions: time(0), y_afa(128), x_faa(128), x_caa(128), y_aca(128)
├── 5 outputs: (ε, ∫εᴰ, ∫ε, KE, ∫KE)
├── array_type: Array{Float32}
├── file_splitting: NoFileSplitting
└── file size: 30.6 KiB

Run the simulation and process results

To run the simulation:

run!(simulation)

[ Info: Initializing simulation...
[ Info: [000.00%] time = 0 seconds,  Δt = 100 ms,  walltime = 20.415 seconds,  advective CFL = 2.4,  diffusive CFL = 0.00042
[ Info:     ... simulation initialization complete (5.512 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (2.154 seconds).
[ Info: [014.73%] time = 7.365 seconds,  Δt = 90.641 ms,  walltime = 26.887 seconds,  advective CFL = 0.79,  diffusive CFL = 0.00038
[ Info: [030.70%] time = 15.348 seconds,  Δt = 92.870 ms,  walltime = 28.311 seconds,  advective CFL = 0.8,  diffusive CFL = 0.00039
[ Info: [048.00%] time = 24 seconds,  Δt = 100.636 ms,  walltime = 29.782 seconds,  advective CFL = 0.8,  diffusive CFL = 0.00042
[ Info: [066.00%] time = 33 seconds,  Δt = 98.999 ms,  walltime = 31.294 seconds,  advective CFL = 0.8,  diffusive CFL = 0.00041
[ Info: [083.94%] time = 41.970 seconds,  Δt = 95.806 ms,  walltime = 32.805 seconds,  advective CFL = 0.8,  diffusive CFL = 0.0004
[ Info: Simulation is stopping after running for 16.603 seconds.
[ Info: Simulation time 50 seconds equals or exceeds stop time 50 seconds.

Now we'll read the results using FieldTimeSeries

filepath = simulation.output_writers[:nc].filepath
KE_t = FieldTimeSeries(filepath, "KE")
ε_t  = FieldTimeSeries(filepath, "ε")

128×128×1×84 FieldTimeSeries{InMemory} located at (Center, Center, Center) of ε at /home/runner/work/Oceanostics.jl/Oceanostics.jl/docs/build/generated/two_dimensional_turbulence.nc
├── grid: 128×128×1 RectilinearGrid{Float64, Periodic, Periodic, Flat} on CPU with 3×3×0 halo
├── indices: (1:128, 1:128, 1:1)
├── time_indexing: Linear()
├── backend: InMemory()
├── path: /home/runner/work/Oceanostics.jl/Oceanostics.jl/docs/build/generated/two_dimensional_turbulence.nc
├── name: ε
└── data: 128×128×1×84 OffsetArray(::Array{Float64, 4}, 1:128, 1:128, 1:1, 1:84) with eltype Float64 with indices 1:128×1:128×1:1×1:84
    └── max=0.0128719, min=5.79802e-10, mean=3.52058e-5

Volume-integrated quantities are scalar time series, so we read them directly with NCDatasets:

ds = NCDataset(filepath)
∫KE = ds["∫KE"][:]
∫ε  = ds["∫ε"][:]
∫εᴰ = ds["∫εᴰ"][:]
close(ds)

closed Dataset

In order to plot results, we use Makie.jl, which has recipes for Oceananigans Fields.

using CairoMakie

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

axis_kwargs = (xlabel = "x", ylabel = "y",
               aspect = DataAspect(),
               height = 300, width = 300)

ax1 = Axis(fig[2, 1]; title = "Kinetic energy", axis_kwargs...)
ax2 = Axis(fig[2, 2]; title = "Kinetic energy dissip rate", axis_kwargs...)

Makie.Axis with 0 plots:

Now we plot the snapshots and set the title

n = Observable(1)

Observable(1)

n above is a Makie.Observable, which allows us to animate things easily. Creating observable KE and ε can be done simply with

KEₙ = @lift KE_t[$n]
εₙ  = @lift ε_t[$n]

Observable(128×128×1 Field{Center, Center, Center} on RectilinearGrid on CPU
├── grid: 128×128×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
├── indices: (1:128, 1:128, 1:1)
└── data: 128×128×1 OffsetArray(view(::Array{Float64, 4}, :, :, :, 1), 1:128, 1:128, 1:1) with eltype Float64 with indices 1:128×1:128×1:1
    └── max=0.0128719, min=1.32868e-5, mean=0.00137467)

Now we plot the heatmaps, each with its own colorbar below

hm_KE = heatmap!(ax1, KEₙ, colormap = :plasma, colorrange=(0, 5e-2))
Colorbar(fig[3, 1], hm_KE; vertical=false, height=8, ticklabelsize=12)

hm_ε = heatmap!(ax2, εₙ, colormap = :inferno, colorrange=(0, 5e-5))
Colorbar(fig[3, 2], hm_ε; vertical=false, height=8, ticklabelsize=12)

Makie.Colorbar()

We now plot the time evolution of our integrated quantities

axis_kwargs = (xlabel = "Time",
               height=150, width=300)

ax3 = Axis(fig[4, 1]; axis_kwargs...)
times = KE_t.times
lines!(ax3, times, ∫KE)

ax4 = Axis(fig[4, 2]; axis_kwargs...)
lines!(ax4, times, ∫ε,  label="∫εdV")
lines!(ax4, times, ∫εᴰ, label="∫εᴰdV", linestyle=:dash)
axislegend(ax4, labelsize=14)

Makie.Legend()

Now we mark the time by placing a vertical line in the bottom plots:

tₙ = @lift times[$n]
vlines!(ax3, tₙ, color=:black, linestyle=:dash)
vlines!(ax4, tₙ, color=:black, linestyle=:dash)

Makie.VLines{Tuple{Float64}}

and by creating a useful title

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

Next we adjust the total figure size based on our panels, which makes it look like this

resize_to_layout!(fig)
fig

Finally, we record a movie.

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

[ Info: Animating...

Although simple, this example and the animation above already illustrate a couple of interesting things. First, the KE dissipation rate ε is distributed at much smaller scales than the KE, which is expected due to the second-order derivatives present in ε.

Second, again as expected, the volume-integrated KE dissipation rate is the same as the volume-integrated KE diffusion term (since all the non-dissipation parts of the term volume-integrate to zero). In fact, both KineticEnergyDissipationRate and KineticEnergyStress in Oceanostics are implemented in an energy-conserving form (i.e., they use the exact same discretization scheme and interpolations as used in Oceananigans), so they agree to machine-precision, and are great for closing budgets.

This page was generated using Literate.jl.