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 (8.102 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (2.491 seconds).
Iter:     50, sim time:  0.544, wall time: 10.729 seconds, Δt:  0.0121, advective CFL: 4.38e-03, diffusive CFL: 1.40e-02
Iter:    100, sim time:  1.145, wall time: 10.945 seconds, Δt:  0.0133, advective CFL: 1.66e-02, diffusive CFL: 1.54e-02
Iter:    150, sim time:  1.806, wall time: 11.100 seconds, Δt:  0.0146, advective CFL: 3.74e-02, diffusive CFL: 1.70e-02
Iter:    200, sim time:  2.515, wall time: 11.258 seconds, Δt:  0.0161, advective CFL: 6.85e-02, diffusive CFL: 1.87e-02
Iter:    250, sim time:  3.306, wall time: 11.403 seconds, Δt:  0.0177, advective CFL: 1.16e-01, diffusive CFL: 2.05e-02
Iter:    300, sim time:  4.177, wall time: 11.600 seconds, Δt:  0.0195, advective CFL: 2.05e-01, diffusive CFL: 2.26e-02
Iter:    350, sim time:  5.136, wall time: 11.803 seconds, Δt:  0.0214, advective CFL: 5.84e-01, diffusive CFL: 2.48e-02
Iter:    400, sim time:  6.193, wall time: 12.043 seconds, Δt:  0.0165, advective CFL: 7.00e-01, diffusive CFL: 1.91e-02
Iter:    450, sim time:  7.000, wall time: 12.208 seconds, Δt:  0.0160, advective CFL: 7.00e-01, diffusive CFL: 1.85e-02
Iter:    500, sim time:  7.788, wall time: 12.378 seconds, Δt:  0.0157, advective CFL: 7.00e-01, diffusive CFL: 1.82e-02
Iter:    550, sim time:  8.563, wall time: 12.579 seconds, Δt:  0.0155, advective CFL: 7.00e-01, diffusive CFL: 1.80e-02
Iter:    600, sim time:  9.326, wall time: 12.777 seconds, Δt:  0.0167, advective CFL: 7.00e-01, diffusive CFL: 1.93e-02
Iter:    650, sim time:  10.150, wall time: 12.992 seconds, Δt:  0.0155, advective CFL: 7.00e-01, diffusive CFL: 1.80e-02
Iter:    700, sim time:  10.920, wall time: 13.190 seconds, Δt:  0.0159, advective CFL: 7.00e-01, diffusive CFL: 1.85e-02
Iter:    750, sim time:  11.691, wall time: 13.351 seconds, Δt:  0.0174, advective CFL: 7.00e-01, diffusive CFL: 2.02e-02
Iter:    800, sim time:  12.552, wall time: 13.583 seconds, Δt:  0.0179, advective CFL: 7.00e-01, diffusive CFL: 2.07e-02
Iter:    850, sim time:  13.429, wall time: 13.765 seconds, Δt:  0.0181, advective CFL: 7.00e-01, diffusive CFL: 2.09e-02
Iter:    900, sim time:  14.325, wall time: 13.928 seconds, Δt:  0.0185, advective CFL: 7.00e-01, diffusive CFL: 2.15e-02
Iter:    950, sim time:  15.241, wall time: 14.089 seconds, Δt:  0.0200, advective CFL: 7.00e-01, diffusive CFL: 2.31e-02
Iter:   1000, sim time:  16.220, wall time: 14.236 seconds, Δt:  0.0199, advective CFL: 7.00e-01, diffusive CFL: 2.31e-02
Iter:   1050, sim time:  17.179, wall time: 14.365 seconds, Δt:  0.0199, advective CFL: 7.00e-01, diffusive CFL: 2.31e-02
Iter:   1100, sim time:  18.139, wall time: 14.488 seconds, Δt:  0.0203, advective CFL: 7.00e-01, diffusive CFL: 2.35e-02
Iter:   1150, sim time:  19.142, wall time: 14.702 seconds, Δt:  0.0211, advective CFL: 7.00e-01, diffusive CFL: 2.45e-02
Iter:   1200, sim time:  20.190, wall time: 14.849 seconds, Δt:  0.0223, advective CFL: 7.00e-01, diffusive CFL: 2.59e-02
Iter:   1250, sim time:  21.290, wall time: 15.003 seconds, Δt:  0.0215, advective CFL: 7.00e-01, diffusive CFL: 2.49e-02
Iter:   1300, sim time:  22.344, wall time: 15.157 seconds, Δt:  0.0202, advective CFL: 7.00e-01, diffusive CFL: 2.34e-02
Iter:   1350, sim time:  23.344, wall time: 15.362 seconds, Δt:  0.0207, advective CFL: 7.00e-01, diffusive CFL: 2.40e-02
Iter:   1400, sim time:  24.352, wall time: 15.590 seconds, Δt:  0.0226, advective CFL: 7.00e-01, diffusive CFL: 2.62e-02
Iter:   1450, sim time:  25.453, wall time: 15.738 seconds, Δt:  0.0241, advective CFL: 7.00e-01, diffusive CFL: 2.79e-02
Iter:   1500, sim time:  26.645, wall time: 15.926 seconds, Δt:  0.0265, advective CFL: 6.86e-01, diffusive CFL: 3.07e-02
Iter:   1550, sim time:  27.951, wall time: 16.082 seconds, Δt:  0.0213, advective CFL: 7.00e-01, diffusive CFL: 2.47e-02
Iter:   1600, sim time:  28.989, wall time: 16.300 seconds, Δt:  0.0171, advective CFL: 7.00e-01, diffusive CFL: 1.98e-02
Iter:   1650, sim time:  29.824, wall time: 16.512 seconds, Δt:  0.0157, advective CFL: 7.00e-01, diffusive CFL: 1.82e-02
Iter:   1700, sim time:  30.594, wall time: 16.732 seconds, Δt:  0.0144, advective CFL: 7.00e-01, diffusive CFL: 1.67e-02
Iter:   1750, sim time:  31.303, wall time: 16.869 seconds, Δt:  0.0157, advective CFL: 7.00e-01, diffusive CFL: 1.82e-02
Iter:   1800, sim time:  32.078, wall time: 17.089 seconds, Δt:  0.0173, advective CFL: 6.43e-01, diffusive CFL: 2.00e-02
Iter:   1850, sim time:  32.932, wall time: 17.253 seconds, Δt:  0.0190, advective CFL: 6.97e-01, diffusive CFL: 2.20e-02
Iter:   1900, sim time:  33.861, wall time: 17.443 seconds, Δt:  0.0207, advective CFL: 7.00e-01, diffusive CFL: 2.40e-02
Iter:   1950, sim time:  34.873, wall time: 17.610 seconds, Δt:  0.0214, advective CFL: 7.00e-01, diffusive CFL: 2.48e-02
Iter:   2000, sim time:  35.928, wall time: 17.822 seconds, Δt:  0.0221, advective CFL: 7.00e-01, diffusive CFL: 2.56e-02
Iter:   2050, sim time:  37.000, wall time: 17.978 seconds, Δt:  0.0226, advective CFL: 7.00e-01, diffusive CFL: 2.61e-02
Iter:   2100, sim time:  38.090, wall time: 18.177 seconds, Δt:  0.0230, advective CFL: 7.00e-01, diffusive CFL: 2.66e-02
Iter:   2150, sim time:  39.230, wall time: 18.319 seconds, Δt:  0.0235, advective CFL: 7.00e-01, diffusive CFL: 2.72e-02
[ Info: Simulation is stopping after running for 18.436 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-34931/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.0
Commit d1c37793dd2 (2026-09-09 19:00 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_LOAD_PATH = @:@v#.#:@stdlib
  JULIA_DEPOT_PATH = /var/lib/buildkite-agent/.julia:/var/lib/buildkite-agent/.julia/juliaup/julia-1.13.0+0.x64.linux.gnu/local/share/julia:/var/lib/buildkite-agent/.julia/juliaup/julia-1.13.0+0.x64.linux.gnu/share/julia
  JULIA_VERSION_ENZYME = 1.13.0
  JULIA_PKG_SERVER_REGISTRY_PREFERENCE = eager
  LD_LIBRARY_PATH = 
  JULIA_MAX_NUM_PRECOMPILE_FILES = 24
  JULIA_VERSION = 1.13.0
  JULIA_CUDA_USE_COMPAT = false
  JULIA_PROJECT = /var/lib/buildkite-agent/Oceananigans.jl-34931/docs/
  JULIA_DEBUG = Literate

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

julia
import Pkg
Pkg.status()
Status `~/Oceananigans.jl-34931/docs/Project.toml`
  [79e6a3ab] Adapt v4.7.4
⌃ [052768ef] CUDA v6.2.2
  [13f3f980] CairoMakie v0.15.16
  [8e50ac2c] ConservativeRegridding v0.2.16
⌅ [e30172f5] Documenter v1.17.0
  [daee34ce] DocumenterCitations v1.5.0
  [4710194d] DocumenterVitepress v0.3.8
  [7da242da] Enzyme v0.13.216
  [033835bb] JLD2 v0.6.7
⌅ [63c18a36] KernelAbstractions v0.9.45
  [98b081ad] Literate v2.21.0
  [da04e1cc] MPI v0.20.29
  [85f8d34a] NCDatasets v0.14.15
  [9e8cae18] Oceananigans v0.114.0 `..`
  [429524aa] Optim v2.3.2
  [f27b6e38] Polynomials v4.1.3
⌅ [3c362404] Reactant v0.2.292
  [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.