Oceananigans.jl

Appearance

Horizontal convection example ​

In "horizontal convection", a non-uniform buoyancy is imposed on top of an initially resting fluid.

This example demonstrates:

  • How to use computed Fields for output.

  • How to post-process saved output using FieldTimeSeries.

Install dependencies ​

First let's make sure we have all required packages installed.

julia
using Pkg
pkg"add Oceananigans, CairoMakie"

Horizontal convection ​

We consider two-dimensional horizontal convection of an incompressible flow   on the -plane (   and   ). The flow evolves under the effect of gravity. The only forcing on the fluid comes from a prescribed, non-uniform buoyancy at the top-surface of the domain.

We start by importing Oceananigans and Printf.

julia
using Oceananigans
using Printf

The grid ​

julia
H = 1            # vertical domain extent
Lx = 2H          # horizontal domain extent
Nx, Nz = 128, 64 # horizontal, vertical resolution

grid = RectilinearGrid(size = (Nx, Nz),
                       x = (-Lx/2, Lx/2),
                       z = (-H, 0),
                       topology = (Bounded, Flat, Bounded))
128×1×64 RectilinearGrid{Float64, Bounded, Flat, Bounded} on CPU with 3×0×3 halo
├── Bounded  x ∈ [-1.0, 1.0] regularly spaced with Δx=0.015625
├── Flat y                   
└── Bounded  z ∈ [-1.0, 0.0] regularly spaced with Δz=0.015625

Boundary conditions ​

At the surface, the imposed buoyancy is

while zero-flux boundary conditions are imposed on all other boundaries. We use free-slip boundary conditions on and everywhere.

julia
b★ = 1
@inline bˢ(x, t, p) = - p.b★ * cos(2π * x / p.Lx)

b_bcs = FieldBoundaryConditions(top = ValueBoundaryCondition(bˢ, parameters=(; b★, Lx)))
Oceananigans.FieldBoundaryConditions, with boundary conditions
├── west: DefaultBoundaryCondition (FluxBoundaryCondition: nothing)
├── east: DefaultBoundaryCondition (FluxBoundaryCondition: nothing)
├── south: DefaultBoundaryCondition (FluxBoundaryCondition: nothing)
├── north: DefaultBoundaryCondition (FluxBoundaryCondition: nothing)
├── bottom: DefaultBoundaryCondition (FluxBoundaryCondition: nothing)
├── top: ValueBoundaryCondition: ContinuousBoundaryFunction bˢ at (Nothing, Nothing, Nothing)
└── immersed: DefaultBoundaryCondition (FluxBoundaryCondition: nothing)

Non-dimensional control parameters and Turbulence closure ​

The problem is characterized by three non-dimensional parameters. The first is the domain's aspect ratio, and the other two are the Rayleigh () and Prandtl () numbers:

The Prandtl number expresses the ratio of momentum over heat diffusion; the Rayleigh number is a measure of the relative importance of gravity over viscosity in the momentum equation.

For a domain with a given extent, the nondimensional values of and uniquely determine the viscosity and diffusivity, i.e.,

We use isotropic viscosity and diffusivity, ν and κ, whose values are obtained from the prescribed and numbers. Here, we use   and  :

julia
Pr = 1
Ra = 1e8

ν = sqrt(Pr * b★ * Lx^3 / Ra)
κ = sqrt(b★ * Lx^3 / (Pr * Ra))

Model instantiation ​

We instantiate the model with the fifth-order WENO advection scheme, a 3rd order Runge-Kutta time-stepping scheme, and a BuoyancyTracer.

julia
model = NonhydrostaticModel(grid;
                            advection = WENO(),
                            timestepper = :RungeKutta3,
                            tracers = :b,
                            buoyancy = BuoyancyTracer(),
                            closure = ScalarDiffusivity(; ν, κ),
                            boundary_conditions = (; b=b_bcs))
NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 128×1×64 RectilinearGrid{Float64, Bounded, Flat, Bounded} on CPU with 3×0×3 halo
├── timestepper: RungeKutta3TimeStepper
├── advection scheme:
│   ├── momentum: WENO{3, Float64, Oceananigans.Utils.BackendOptimizedDivision}(order=5)
│   └── b: WENO{3, Float64, Oceananigans.Utils.BackendOptimizedDivision}(order=5)
├── tracers: b
├── closure: ScalarDiffusivity{ExplicitTimeDiscretization}(ν=0.000282843, κ=(b=0.000282843,))
├── buoyancy: BuoyancyTracer with ĝ = NegativeZDirection()
└── coriolis: Nothing

Simulation set-up ​

We set up a simulation that runs up to   with a JLD2Writer that saves the flow speed, , the buoyancy, , and the vorticity,  .

julia
simulation = Simulation(model, Δt=1e-2, stop_time=40)
Simulation of NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── Next time step: 10 ms
├── run_wall_time: 0 seconds
├── run_wall_time / iteration: NaN days
├── stop_time: 40 seconds
├── stop_iteration: Inf
├── wall_time_limit: Inf
├── minimum_relative_step: 0.0
├── callbacks: OrderedDict with 4 entries:
│   ├── stop_time_exceeded => Callback of stop_time_exceeded on IterationInterval(1)
│   ├── stop_iteration_exceeded => Callback of stop_iteration_exceeded on IterationInterval(1)
│   ├── wall_time_limit_exceeded => Callback of wall_time_limit_exceeded on IterationInterval(1)
│   └── nan_checker => Callback of NaNChecker for u on IterationInterval(100)
└── output_writers: OrderedDict with no entries

The TimeStepWizard ​

The TimeStepWizard manages the time-step adaptively, keeping the Courant-Friedrichs-Lewy (CFL) number close to 0.7.

julia
conjure_time_step_wizard!(simulation, IterationInterval(50), cfl=0.7, max_Δt=1e-1)

A progress messenger ​

We write a function that prints out a helpful progress message while the simulation runs.

julia
progress(sim) = @printf("Iter: % 6d, sim time: % 1.3f, wall time: % 10s, Δt: % 1.4f, advective CFL: %.2e, diffusive CFL: %.2e\n",
                        iteration(sim), time(sim), prettytime(sim.run_wall_time),
                        sim.Δt, AdvectiveCFL(sim.Δt)(sim.model), DiffusiveCFL(sim.Δt)(sim.model))

simulation.callbacks[:progress] = Callback(progress, IterationInterval(50))
Callback of progress on IterationInterval(50)

Output ​

We use computed Fields to diagnose and output the total flow speed, the vorticity, , and the buoyancy, . Note that computed Fields take "AbstractOperations" on Fields as input:

julia
u, v, w = model.velocities
b = model.tracers.b

s = @at (Center, Center, Center) sqrt(u^2 + w^2)
ζ = ∂z(u) - ∂x(w)

We add to the simulation a JLD2Writer that saves the speed, the buoyancy, and the vorticity. We want to post-process buoyancy to compute the buoyancy variance dissipation (which is proportional to ). By default JLD2Writer saves the fields' halos, so when we load the buoyancy back it comes with the proper boundary conditions for computing gradients.

julia
filename = "horizontal_convection.jld2"

simulation.output_writers[:fields] = JLD2Writer(model, (; s, b, ζ),
                                                schedule = TimeInterval(0.5),
                                                filename = filename,
                                                overwrite_files = true)

Ready to press the big red button:

julia
run!(simulation)
[ Info: Initializing simulation...
Iter:      0, sim time:  0.000, wall time:  0 seconds, Δt:  0.0110, advective CFL: 0.00e+00, diffusive CFL: 1.27e-02
[ Info:     ... simulation initialization complete (2.704 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (1.113 seconds).
Iter:     50, sim time:  0.544, wall time: 3.930 seconds, Δt:  0.0121, advective CFL: 4.38e-03, diffusive CFL: 1.40e-02
Iter:    100, sim time:  1.145, wall time: 4.046 seconds, Δt:  0.0133, advective CFL: 1.66e-02, diffusive CFL: 1.54e-02
Iter:    150, sim time:  1.806, wall time: 4.143 seconds, Δt:  0.0146, advective CFL: 3.74e-02, diffusive CFL: 1.70e-02
Iter:    200, sim time:  2.515, wall time: 4.245 seconds, Δt:  0.0161, advective CFL: 6.85e-02, diffusive CFL: 1.87e-02
Iter:    250, sim time:  3.306, wall time: 4.342 seconds, Δt:  0.0177, advective CFL: 1.16e-01, diffusive CFL: 2.05e-02
Iter:    300, sim time:  4.177, wall time: 4.445 seconds, Δt:  0.0195, advective CFL: 2.05e-01, diffusive CFL: 2.26e-02
Iter:    350, sim time:  5.136, wall time: 4.568 seconds, Δt:  0.0214, advective CFL: 5.86e-01, diffusive CFL: 2.48e-02
Iter:    400, sim time:  6.193, wall time: 4.676 seconds, Δt:  0.0165, advective CFL: 7.00e-01, diffusive CFL: 1.91e-02
Iter:    450, sim time:  7.000, wall time: 4.775 seconds, Δt:  0.0160, advective CFL: 7.00e-01, diffusive CFL: 1.85e-02
Iter:    500, sim time:  7.788, wall time: 4.881 seconds, Δt:  0.0157, advective CFL: 7.00e-01, diffusive CFL: 1.82e-02
Iter:    550, sim time:  8.563, wall time: 4.988 seconds, Δt:  0.0155, advective CFL: 7.00e-01, diffusive CFL: 1.80e-02
Iter:    600, sim time:  9.326, wall time: 5.103 seconds, Δt:  0.0166, advective CFL: 7.00e-01, diffusive CFL: 1.93e-02
Iter:    650, sim time:  10.133, wall time: 5.211 seconds, Δt:  0.0155, advective CFL: 7.00e-01, diffusive CFL: 1.80e-02
Iter:    700, sim time:  10.904, wall time: 5.317 seconds, Δt:  0.0159, advective CFL: 7.00e-01, diffusive CFL: 1.85e-02
Iter:    750, sim time:  11.675, wall time: 5.429 seconds, Δt:  0.0174, advective CFL: 7.00e-01, diffusive CFL: 2.01e-02
Iter:    800, sim time:  12.535, wall time: 5.534 seconds, Δt:  0.0179, advective CFL: 7.00e-01, diffusive CFL: 2.07e-02
Iter:    850, sim time:  13.430, wall time: 5.642 seconds, Δt:  0.0180, advective CFL: 7.00e-01, diffusive CFL: 2.09e-02
Iter:    900, sim time:  14.325, wall time: 5.756 seconds, Δt:  0.0185, advective CFL: 7.00e-01, diffusive CFL: 2.14e-02
Iter:    950, sim time:  15.222, wall time: 5.860 seconds, Δt:  0.0200, advective CFL: 7.00e-01, diffusive CFL: 2.31e-02
Iter:   1000, sim time:  16.200, wall time: 5.961 seconds, Δt:  0.0199, advective CFL: 7.00e-01, diffusive CFL: 2.31e-02
Iter:   1050, sim time:  17.159, wall time: 6.064 seconds, Δt:  0.0199, advective CFL: 7.00e-01, diffusive CFL: 2.31e-02
Iter:   1100, sim time:  18.119, wall time: 6.167 seconds, Δt:  0.0203, advective CFL: 7.00e-01, diffusive CFL: 2.35e-02
Iter:   1150, sim time:  19.122, wall time: 6.280 seconds, Δt:  0.0211, advective CFL: 7.00e-01, diffusive CFL: 2.45e-02
Iter:   1200, sim time:  20.169, wall time: 6.394 seconds, Δt:  0.0223, advective CFL: 7.00e-01, diffusive CFL: 2.59e-02
Iter:   1250, sim time:  21.268, wall time: 6.500 seconds, Δt:  0.0216, advective CFL: 7.00e-01, diffusive CFL: 2.50e-02
Iter:   1300, sim time:  22.324, wall time: 6.606 seconds, Δt:  0.0203, advective CFL: 7.00e-01, diffusive CFL: 2.35e-02
Iter:   1350, sim time:  23.325, wall time: 6.713 seconds, Δt:  0.0209, advective CFL: 7.00e-01, diffusive CFL: 2.42e-02
Iter:   1400, sim time:  24.355, wall time: 6.824 seconds, Δt:  0.0227, advective CFL: 7.00e-01, diffusive CFL: 2.63e-02
Iter:   1450, sim time:  25.453, wall time: 6.924 seconds, Δt:  0.0243, advective CFL: 7.00e-01, diffusive CFL: 2.81e-02
Iter:   1500, sim time:  26.646, wall time: 7.027 seconds, Δt:  0.0267, advective CFL: 6.88e-01, diffusive CFL: 3.09e-02
Iter:   1550, sim time:  27.954, wall time: 7.120 seconds, Δt:  0.0212, advective CFL: 7.00e-01, diffusive CFL: 2.46e-02
Iter:   1600, sim time:  28.988, wall time: 7.211 seconds, Δt:  0.0171, advective CFL: 7.00e-01, diffusive CFL: 1.98e-02
Iter:   1650, sim time:  29.824, wall time: 7.303 seconds, Δt:  0.0157, advective CFL: 7.00e-01, diffusive CFL: 1.82e-02
Iter:   1700, sim time:  30.594, wall time: 7.395 seconds, Δt:  0.0145, advective CFL: 7.00e-01, diffusive CFL: 1.68e-02
Iter:   1750, sim time:  31.319, wall time: 7.483 seconds, Δt:  0.0158, advective CFL: 7.00e-01, diffusive CFL: 1.83e-02
Iter:   1800, sim time:  32.095, wall time: 7.579 seconds, Δt:  0.0174, advective CFL: 6.48e-01, diffusive CFL: 2.02e-02
Iter:   1850, sim time:  32.953, wall time: 7.675 seconds, Δt:  0.0191, advective CFL: 7.00e-01, diffusive CFL: 2.21e-02
Iter:   1900, sim time:  33.882, wall time: 7.774 seconds, Δt:  0.0208, advective CFL: 7.00e-01, diffusive CFL: 2.41e-02
Iter:   1950, sim time:  34.895, wall time: 7.870 seconds, Δt:  0.0214, advective CFL: 7.00e-01, diffusive CFL: 2.48e-02
Iter:   2000, sim time:  35.950, wall time: 7.969 seconds, Δt:  0.0221, advective CFL: 7.00e-01, diffusive CFL: 2.56e-02
Iter:   2050, sim time:  37.022, wall time: 8.071 seconds, Δt:  0.0226, advective CFL: 7.00e-01, diffusive CFL: 2.62e-02
Iter:   2100, sim time:  38.113, wall time: 8.173 seconds, Δt:  0.0231, advective CFL: 7.00e-01, diffusive CFL: 2.68e-02
Iter:   2150, sim time:  39.254, wall time: 8.276 seconds, Δt:  0.0236, advective CFL: 7.00e-01, diffusive CFL: 2.73e-02
[ Info: Simulation is stopping after running for 8.342 seconds.
[ Info: Simulation time 40 seconds equals or exceeds stop time 40 seconds.

Load saved output, process, visualize ​

We animate the results by loading the saved output and plotting the saved fields. From the saved buoyancy field we compute the buoyancy dissipation,  , and plot that also.

To start we load the saved fields as FieldTimeSeries.

julia
using CairoMakie

s_timeseries = FieldTimeSeries(filename, "s")
b_timeseries = FieldTimeSeries(filename, "b")
ζ_timeseries = FieldTimeSeries(filename, "ζ")

times = b_timeseries.times
Nt = length(times)

χ_timeseries = deepcopy(b_timeseries)

for n in 1:Nt
    b = b_timeseries[n]
    χ_timeseries[n] .= @at (Center, Center, Center) κ * (∂x(b)^2 + ∂z(b)^2)
end
┌ Warning: Assignment to `b` in soft scope is ambiguous because a global variable by the same name exists: `b` will be treated as a new local. Disambiguate by using `local b` to suppress this warning or `global b` to assign to the existing global variable.
└ @ ~/Oceananigans.jl-36221/docs/src/literated/horizontal_convection.md:13

Now we're ready to animate using Makie.

julia
@info "Making an animation from saved data..."

n = Observable(1)

title = @lift @sprintf("t=%1.2f", times[$n])

sn = @lift s_timeseries[$n]
ζn = @lift ζ_timeseries[$n]
bn = @lift b_timeseries[$n]
χn = @lift χ_timeseries[$n]

slim = 0.6
blim = 0.6
ζlim = 9
χlim = 0.025

axis_kwargs = (xlabel = L"x / H",
               ylabel = L"z / H",
               limits = ((-Lx/2, Lx/2), (-H, 0)),
               aspect = Lx / H,
               titlesize = 20)

fig = Figure(size = (600, 1100))

ax_s = Axis(fig[2, 1]; title = L"speed, $(u^2+w^2)^{1/2} / (L_x b_*)^{1/2}$", axis_kwargs...)
ax_b = Axis(fig[3, 1]; title = L"buoyancy, $b / b_*$", axis_kwargs...)
ax_ζ = Axis(fig[4, 1]; axis_kwargs...,
            title = L"vorticity, $(∂u/∂z - ∂w/∂x) \, (L_x / b_*)^{1/2}$")

ax_χ = Axis(fig[5, 1]; axis_kwargs...,
            title = L"buoyancy dissipation, $κ |\mathbf{\nabla}b|^2 \, (L_x / {b_*}^5)^{1/2}$")

fig[1, :] = Label(fig, title, fontsize=24, tellwidth=false)

hm_s = heatmap!(ax_s, sn; colorrange=(0, slim), colormap=:speed)
Colorbar(fig[2, 2], hm_s)

hm_b = heatmap!(ax_b, bn; colorrange=(-blim, blim), colormap=:thermal)
Colorbar(fig[3, 2], hm_b)

hm_ζ = heatmap!(ax_ζ, ζn; colorrange=(-ζlim, ζlim), colormap=:balance)
Colorbar(fig[4, 2], hm_ζ)

hm_χ = heatmap!(ax_χ, χn; colorrange=(0, χlim), colormap=:dense)
Colorbar(fig[5, 2], hm_χ)
Colorbar()

And, finally, we record a movie.

julia
record(fig, "horizontal_convection.mp4", 1:Nt, framerate=8) do i
    n[] = i
end

At higher Rayleigh numbers the flow becomes much more vigorous. See, for example, an animation of the vorticity of the fluid at   on vimeo.

The Nusselt number ​

Often we are interested in how much the flow enhances mixing. This is quantified by the Nusselt number, which compares the mixing achieved by the flow with the mixing that diffusion alone would achieve. The Nusselt number is given by

where angle brackets above denote both a volume and time average and is the buoyancy dissipation that we get without any flow, i.e., the buoyancy dissipation associated with the buoyancy distribution that satisfies

with the same boundary conditions as our setup. In this case, we can readily find that

where is the surface boundary condition. The diffusive solution implies

We use the loaded FieldTimeSeries to compute the instantaneous Nusselt number and the volume-averaged kinetic energy of the fluid.

First we compute the diffusive buoyancy dissipation, (which is just a scalar):

julia
χᵈⁱᶠᶠ = κ * b★^2 * π * tanh(2π * H / Lx) / (Lx * H)

Then we compute the volume averages with Average, and plot.

julia
kinetic_energy = [Field(Average(s_timeseries[n]^2 / 2))[1, 1, 1] for n in 1:Nt]
Nu = [Field(Average(χ_timeseries[n]))[1, 1, 1] / χᵈⁱᶠᶠ for n in 1:Nt]

fig = Figure(size = (850, 450))

ax_KE = Axis(fig[1, 1], xlabel = L"t \, (b_* / L_x)^{1/2}", ylabel = L"KE $ / (L_x b_*)$")
lines!(ax_KE, times, kinetic_energy; linewidth = 3)

ax_Nu = Axis(fig[2, 1], xlabel = L"t \, (b_* / L_x)^{1/2}", ylabel = L"Nu")
lines!(ax_Nu, times, Nu; linewidth = 3)

fig

Julia version and environment information ​

This example was executed with the following version of Julia:

julia
using InteractiveUtils: versioninfo
versioninfo()
Julia Version 1.13.2
Commit 45b6d15e34b (2026-10-14 20:33 UTC)
Build Info:
  Official https://julialang.org release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 128 × AMD EPYC 9374F 32-Core Processor
  WORD_SIZE: 64
  LLVM: libLLVM-20.1.8 (ORCJIT, znver4)
  GC: Built with stock GC
Threads: 1 default, 1 interactive, 1 GC (on 128 virtual cores)
Environment:
  JULIA_VERSION_ENZYME = 1.13.1
  JULIA_DEPOT_PATH = /var/lib/buildkite-agent/.julia:/var/lib/buildkite-agent/.julia/juliaup/julia-1.13.2+0.x64.linux.gnu/local/share/julia:/var/lib/buildkite-agent/.julia/juliaup/julia-1.13.2+0.x64.linux.gnu/share/julia
  JULIA_VERSION = 1.13.2
  JULIA_PKG_SERVER_REGISTRY_PREFERENCE = eager
  LD_LIBRARY_PATH = 
  JULIA_LOAD_PATH = @:@v#.#:@stdlib
  JULIA_MAX_NUM_PRECOMPILE_FILES = 24
  JULIA_CUDA_USE_COMPAT = false
  JULIA_PROJECT = /var/lib/buildkite-agent/Oceananigans.jl-36221/docs/
  JULIA_DEBUG = Literate

These were the top-level packages installed in the environment:

julia
import Pkg
Pkg.status()
Status `~/Oceananigans.jl-36221/docs/Project.toml`
  [79e6a3ab] Adapt v4.8.0
⌃ [052768ef] CUDA v6.4.1
  [13f3f980] CairoMakie v0.15.17
  [8e50ac2c] ConservativeRegridding v0.2.21
  [e30172f5] Documenter v1.20.0
  [daee34ce] DocumenterCitations v1.5.0
  [4710194d] DocumenterVitepress v0.3.9
  [7da242da] Enzyme v0.13.219
  [033835bb] JLD2 v0.6.7
⌅ [63c18a36] KernelAbstractions v0.9.47
  [98b081ad] Literate v2.21.1
  [da04e1cc] MPI v0.20.29
  [85f8d34a] NCDatasets v0.14.15
  [9e8cae18] Oceananigans v0.118.1 `..`
  [429524aa] Optim v2.3.2
  [f27b6e38] Polynomials v4.1.6
  [3c362404] Reactant v0.2.308
  [6038ab10] Rotations v1.7.1
  [d496a93d] SeawaterPolynomials v0.3.11
  [09ab397b] StructArrays v0.7.4
  [bdfc003b] TimesDates v0.3.3
  [0a941bbe] Zarr v0.11.2
  [b77e0a4c] InteractiveUtils v1.11.0
  [37e2e46d] LinearAlgebra v1.13.0
  [44cfe95a] Pkg v1.13.0
Info Packages marked with ⌃ and ⌅ have new versions available. Those with ⌃ may be upgradable, but those with ⌅ are restricted by compatibility constraints from upgrading. To see why use `status --outdated`

This page was generated using Literate.jl.