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: 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: 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.783905, min=-0.693481, mean=-5.62091e-16We 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 haloWe 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 haloThe 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 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 = 18.722 seconds, advective CFL = 2.4, diffusive CFL = 0.00042
[ Info: ... simulation initialization complete (11.725 seconds)
[ Info: Executing initial time step...
[ Info: ... initial time step complete (4.010 seconds).
[ Info: [014.56%] time = 7.280 seconds, Δt = 76.300 ms, walltime = 34.388 seconds, advective CFL = 0.8, diffusive CFL = 0.00032
[ Info: [030.38%] time = 15.188 seconds, Δt = 94.843 ms, walltime = 37.102 seconds, advective CFL = 0.8, diffusive CFL = 0.00039
[ Info: [048.41%] time = 24.207 seconds, Δt = 104.168 ms, walltime = 40.074 seconds, advective CFL = 0.8, diffusive CFL = 0.00043
[ Info: [067.07%] time = 33.537 seconds, Δt = 108.163 ms, walltime = 43.028 seconds, advective CFL = 0.8, diffusive CFL = 0.00045
[ Info: [085.79%] time = 42.895 seconds, Δt = 94.857 ms, walltime = 46.121 seconds, advective CFL = 0.8, diffusive CFL = 0.00039
[ Info: Simulation is stopping after running for 32.395 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-23T11:28:34.342…
"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
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.