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: NothingLet'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.673004, min=-0.702298, mean=-4.79488e-17We 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 haloWe 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 haloThe idea in calculating this term is that, in integrated form, all transport contributions in it should equal zero and ∫εᴰ should equal ∫ε.
To illustrate the effect of spatial low-pass filtering, we also compute two related quantities: the GaussianFilter applied to KE, and the KineticEnergy computed from the GaussianFiltered velocities. In general these two quantities differ, since filtering and the (nonlinear) KE operator do not commute.
σ = π/8
KE_filt = GaussianFilter(KE, dims=(1, 2), σ=σ)
u_filt = GaussianFilter(u, dims=(1, 2), σ=σ)
v_filt = GaussianFilter(v, dims=(1, 2), σ=σ)
KE_of_filt = KineticEnergyEquation.KineticEnergy(model, u_filt, v_filt, w)KernelFunctionOperation at (Center, Center, Center)
├── grid: 128×128×1 RectilinearGrid{Float64, Periodic, Periodic, Flat} on CPU with 3×3×0 halo
├── kernel_function: kinetic_energy_ccc (generic function with 1 method)
└── arguments: ("KernelFunctionOperation", "KernelFunctionOperation", "Field")We output the previous quantities to a NetCDF file
output_fields = (; KE, ε, KE_filt, KE_of_filt, ∫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)
├── 7 outputs: (KE_filt, ε, ∫εᴰ, ∫ε, KE, KE_of_filt, ∫KE)
├── array_type: Array{Float32}
├── file_splitting: NoFileSplitting
└── file size: 32.0 KiBRun the simulation and process results
To run the simulation:
run!(simulation)[ Info: Initializing simulation...
[ Info: [000.00%] time = 0 seconds, Δt = 100 ms, walltime = 22.336 seconds, advective CFL = 2.3, diffusive CFL = 0.00042
[ Info: ... simulation initialization complete (4.756 seconds)
[ Info: Executing initial time step...
[ Info: ... initial time step complete (1.903 seconds).
[ Info: [014.95%] time = 7.477 seconds, Δt = 91.939 ms, walltime = 30.825 seconds, advective CFL = 0.8, diffusive CFL = 0.00038
[ Info: [032.11%] time = 16.055 seconds, Δt = 89.107 ms, walltime = 33.987 seconds, advective CFL = 0.8, diffusive CFL = 0.00037
[ Info: [049.20%] time = 24.600 seconds, Δt = 91.011 ms, walltime = 37.146 seconds, advective CFL = 0.8, diffusive CFL = 0.00038
[ Info: [066.53%] time = 33.267 seconds, Δt = 88.187 ms, walltime = 40.505 seconds, advective CFL = 0.8, diffusive CFL = 0.00037
[ Info: [083.68%] time = 41.842 seconds, Δt = 89.026 ms, walltime = 43.670 seconds, advective CFL = 0.8, diffusive CFL = 0.00037
[ Info: Simulation is stopping after running for 25.583 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, "ε")
KE_filt_t = FieldTimeSeries(filepath, "KE_filt")
KE_of_filt_t = FieldTimeSeries(filepath, "KE_of_filt")128×128×1×84 FieldTimeSeries{InMemory} located at (Center, Center, Center) of KE_of_filt 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: KE_of_filt
└── 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.0100097, min=2.06041e-9, mean=0.000723009Volume-integrated quantities are scalar time series, so we read them directly with NCDatasets:
ds = NCDataset(filepath)
∫KE = ds["∫KE"][:]
∫ε = ds["∫ε"][:]
∫εᴰ = ds["∫εᴰ"][:]
close(ds)closed DatasetIn 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...)
ax5 = Axis(fig[4, 1]; title = "GaussianFilter(KE)", axis_kwargs...)
ax6 = Axis(fig[4, 2]; title = "KE(GaussianFiltered velocities)", 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]
KE_filtₙ = @lift KE_filt_t[$n]
KE_of_filtₙ = @lift KE_of_filt_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.00052236, min=1.57507e-8, mean=7.70095e-5)
Now we plot the heatmaps, each with its own colorbar below. The filtered KE panels use the same colormap as KE but a tighter colorrange to emphasize their differences.
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)
hm_KE_filt = heatmap!(ax5, KE_filtₙ, colormap = :plasma, colorrange=(0, 2e-2))
Colorbar(fig[5, 1], hm_KE_filt; vertical=false, height=8, ticklabelsize=12)
hm_KE_of_filt = heatmap!(ax6, KE_of_filtₙ, colormap = :plasma, colorrange=(0, 2e-2))
Colorbar(fig[5, 2], hm_KE_of_filt; 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[6, 1]; axis_kwargs...)
times = KE_t.times
lines!(ax3, times, ∫KE)
ax4 = Axis(fig[6, 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 and we 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...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.