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, Rasters"
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 = UpwindBiasedFifthOrder(),
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: Upwind Biased reconstruction 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: ZeroFlux
└── 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.69819, min=-0.645101, mean=5.52862e-16
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 = KineticEnergy(model)
ε = KineticEnergyDissipationRate(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", "(u=128×128×1 Field{Face, Center, Center} on RectilinearGrid on CPU, v=128×128×1 Field{Center, Face, Center} on RectilinearGrid on CPU, w=128×128×1 Field{Center, Center, Face} on RectilinearGrid on CPU)", "(closure=ScalarDiffusivity{ExplicitTimeDiscretization}(ν=1.0e-5), clock=Clock(time=0 seconds, iteration=0, last_Δt=Inf days), buoyancy=Nothing)")
And we can define their volume-integrals
∫KE = Integral(KE)
∫ε = Integral(ε)
sum! over dims (1, 2, 3) of BinaryOperation at (Center, Center, Center)
└── 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 KineticEnergyStressTerm
, 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(KineticEnergyStressTerm(model))
sum! over dims (1, 2, 3) of BinaryOperation at (Center, Center, Center)
└── 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, ∫ε, ∫εᴰ)
filename = "two_dimensional_turbulence"
simulation.output_writers[:nc] = NetCDFOutputWriter(model, output_fields,
filename = joinpath(@__DIR__, filename),
schedule = TimeInterval(0.6),
overwrite_existing = true)
NetCDFOutputWriter scheduled on TimeInterval(600 ms):
├── filepath: /home/runner/work/Oceanostics.jl/Oceanostics.jl/docs/build/generated/two_dimensional_turbulence.nc
├── dimensions: zC(1), zF(1), xC(128), yF(128), xF(128), yC(128), time(0)
├── 5 outputs: (ε, ∫εᴰ, ∫ε, KE, ∫KE)
└── array type: Array{Float64}
├── file_splitting: NoFileSplitting
└── file size: 19.4 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 = 18.514 seconds, advective CFL = 2.3, diffusive CFL = 0.00042
[ Info: ... simulation initialization complete (11.669 seconds)
[ Info: Executing initial time step...
[ Info: ... initial time step complete (3.974 seconds).
[ Info: [014.40%] time = 7.200 seconds, Δt = 83.829 ms, walltime = 33.790 seconds, advective CFL = 0.8, diffusive CFL = 0.00035
[ Info: [030.92%] time = 15.461 seconds, Δt = 93.847 ms, walltime = 36.388 seconds, advective CFL = 0.8, diffusive CFL = 0.00039
[ Info: [048.00%] time = 24.000 seconds, Δt = 103.104 ms, walltime = 38.982 seconds, advective CFL = 0.8, diffusive CFL = 0.00043
[ Info: [064.80%] time = 32.400 seconds, Δt = 82.237 ms, walltime = 41.573 seconds, advective CFL = 0.8, diffusive CFL = 0.00034
[ Info: [082.65%] time = 41.323 seconds, Δt = 103.834 ms, walltime = 44.306 seconds, advective CFL = 0.8, diffusive CFL = 0.00043
[ Info: [098.40%] time = 49.200 seconds, Δt = 79.356 ms, walltime = 46.784 seconds, advective CFL = 0.8, diffusive CFL = 0.00033
[ Info: Simulation is stopping after running for 31.302 seconds.
[ Info: Simulation time 50 seconds equals or exceeds stop time 50 seconds.
Now we'll read the results using Rasters.jl, which works somewhat similarly to Python's Xarray and can speed-up the workflow
using Rasters
ds = RasterStack(simulation.output_writers[:nc].filepath)
RasterStack with dimensions:
Dim{:xC} Sampled{Float64} Float64[0.02454369260617025, 0.07363107781851076, …, 6.209554229361075, 6.258641614573415] ForwardOrdered Regular Points,
Dim{:yC} Sampled{Float64} Float64[0.02454369260617025, 0.07363107781851076, …, 6.209554229361075, 6.258641614573415] ForwardOrdered Regular Points,
Dim{:zC} Sampled{Float64} Float64[1.0] ForwardOrdered Regular Points,
Ti Sampled{Float64} Float64[0.0, 0.6, …, 49.20000000000001, 49.80000000000001] ForwardOrdered Regular Points
and 5 layers:
:ε Float64 dims: Dim{:xC}, Dim{:yC}, Dim{:zC}, Ti (128×128×1×84)
:∫εᴰ Float64 dims: Ti (84)
:∫ε Float64 dims: Ti (84)
:KE Float64 dims: Dim{:xC}, Dim{:yC}, Dim{:zC}, Ti (128×128×1×84)
:∫KE Float64 dims: Ti (84)
with metadata Metadata{Rasters.NCDsource} of Dict{String, Any} with 6 entries:
"interval" => 0.6
"Oceananigans" => "This file was generated using "
"Julia" => "This file was generated using Julia Version 1.9.4\…
"output time interval" => "Output was saved every 600 ms."
"date" => "This file was generated on 2024-05-09T17:29:48.863…
"schedule" => "TimeInterval"
In order to plot results, we use Makie.jl, for which Rasters.jl already has some recipes
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...)
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 set(ds.KE[zC=1, Ti=$n], :xC => X, :yC => Y);
εₙ = @lift set(ds.ε[zC=1, Ti=$n], :xC => X, :yC => Y);
Note that, in Rasters, the time
coordinate gets shortened to Ti
.
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 = dims(ds, :Ti)
lines!(ax3, Array(times), ds.∫KE)
ax4 = Axis(fig[4, 2]; axis_kwargs...)
lines!(ax4, Array(times), ds.∫ε, label="∫εdV")
lines!(ax4, Array(times), ds.∫εᴰ, 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)
MakieCore.Plot{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
![Example block output](f561f8b1.png)
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 KineticEnergyStressTerm
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.