# 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 ComputedFields for output.
• How to post-process saved output using FieldTimeSeries.

## Install dependencies

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

using Pkg
pkg"add Oceananigans, JLD2, Plots"

## Horizontal convection

We consider here two-dimensional horizontal convection of an incompressible flow $\boldsymbol{u} = (u, w)$ on the $(x, z)$-plane ($-L_x/2 \le x \le L_x/2$ and $-H \le z \le 0$). 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.

using Oceananigans
using Printf

### The grid

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

grid = RegularRectilinearGrid(size = (Nx, Nz),
x = (-Lx/2, Lx/2),
z = (-H, 0),
halo = (3, 3),
topology = (Bounded, Flat, Bounded))
RegularRectilinearGrid{Float64, Bounded, Flat, Bounded}
domain: x ∈ [-1.0, 1.0], y ∈ [0.0, 0.0], z ∈ [-1.0, 0.0]
topology: (Bounded, Flat, Bounded)
resolution (Nx, Ny, Nz): (128, 1, 64)
halo size (Hx, Hy, Hz): (3, 0, 3)
grid spacing (Δx, Δy, Δz): (0.015625, 0.0, 0.015625)

### Boundary conditions

At the surface, the imposed buoyancy is

$$$b(x, z = 0, t) = - b_* \cos (2 \pi x / L_x) \, ,$$$

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

b★ = 1.0

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

b_bcs = FieldBoundaryConditions(top = ValueBoundaryCondition(bˢ, parameters=(; b★, Lx)))
Oceananigans.FieldBoundaryConditions, with boundary conditions
├── west: Oceananigans.BoundaryConditions.DefaultPrognosticFieldBoundaryCondition
├── east: Oceananigans.BoundaryConditions.DefaultPrognosticFieldBoundaryCondition
├── south: Oceananigans.BoundaryConditions.DefaultPrognosticFieldBoundaryCondition
├── north: Oceananigans.BoundaryConditions.DefaultPrognosticFieldBoundaryCondition
├── bottom: Oceananigans.BoundaryConditions.DefaultPrognosticFieldBoundaryCondition
├── top: BoundaryCondition{Value, Oceananigans.BoundaryConditions.ContinuousBoundaryFunction{Nothing, Nothing, Nothing, Nothing, typeof(Main.bˢ), NamedTuple{(:b★, :Lx), Tuple{Float64, Float64}}, Tuple{}, Nothing, Nothing}}
└── immersed: BoundaryCondition{Flux, 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, $L_x / H$ and the other two are the Rayleigh ($Ra$) and Prandtl ($Pr$) numbers:

$$$Ra = \frac{b_* L_x^3}{\nu \kappa} \, , \quad \text{and}\, \quad Pr = \frac{\nu}{\kappa} \, .$$$

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 $Ra$ and $Pr$ uniquely determine the viscosity and diffusivity, i.e.,

$$$\nu = \sqrt{\frac{Pr b_* L_x^3}{Ra}} \quad \text{and} \quad \kappa = \sqrt{\frac{b_* L_x^3}{Pr Ra}} \, .$$$

We use isotropic viscosity and diffusivities, ν and κ whose values are obtain from the prescribed $Ra$ and $Pr$ numbers. Here, we use $Pr = 1$ and $Ra = 10^8$:

Pr = 1.0    # Prandtl number
Ra = 1e8    # Rayleigh number

ν = sqrt(Pr * b★ * Lx^3 / Ra)  # Laplacian viscosity
κ = ν * Pr                     # Laplacian diffusivity

## Model instantiation

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

model = NonhydrostaticModel(
architecture = CPU(),
grid = grid,
timestepper = :RungeKutta3,
tracers = :b,
buoyancy = BuoyancyTracer(),
closure = IsotropicDiffusivity(ν=ν, κ=κ),
boundary_conditions = (b=b_bcs,)
)
NonhydrostaticModel{CPU, Float64}(time = 0 seconds, iteration = 0)
├── grid: RegularRectilinearGrid{Float64, Bounded, Flat, Bounded}(Nx=128, Ny=1, Nz=64)
├── tracers: (:b,)
├── closure: IsotropicDiffusivity{Oceananigans.TurbulenceClosures.ExplicitTimeDiscretization, Float64, NamedTuple{(:b,), Tuple{Float64}}}
├── buoyancy: BuoyancyTracer
└── coriolis: Nothing

## Simulation set-up

We set up a simulation that runs up to $t = 40$ with a JLD2OutputWriter that saves the flow speed, $\sqrt{u^2 + w^2}$, the buoyancy, $b$, andthe vorticity, $\partial_z u - \partial_x w$.

### The TimeStepWizard

The TimeStepWizard manages the time-step adaptively, keeping the Courant-Freidrichs-Lewy (CFL) number close to 0.75 while ensuring the time-step does not increase beyond the maximum allowable value for numerical stability.

max_Δt = 1e-1
wizard = TimeStepWizard(cfl=0.75, Δt=1e-2, max_change=1.2, max_Δt=max_Δt)
TimeStepWizard{Float64, typeof(Oceananigans.Utils.cell_advection_timescale), typeof(Oceananigans.Simulations.infinite_diffusion_timescale)}(0.75, Inf, 1.2, 0.5, 0.1, 0.0, 0.01, Oceananigans.Utils.cell_advection_timescale, Oceananigans.Simulations.infinite_diffusion_timescale)

### A progress messenger

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

CFL = AdvectiveCFL(wizard)

start_time = time_ns()

progress(sim) = @printf("i: % 6d, sim time: % 1.3f, wall time: % 10s, Δt: % 1.4f, CFL: %.2e\n",
sim.model.clock.iteration,
sim.model.clock.time,
prettytime(1e-9 * (time_ns() - start_time)),
sim.Δt.Δt,
CFL(sim.model))

### Build the simulation

We're ready to build and run the simulation. We ask for a progress message and time-step update every 50 iterations,

simulation = Simulation(model, Δt = wizard, iteration_interval = 50,
stop_time = 40.0,
progress = progress)
Simulation{typename(NonhydrostaticModel){typename(CPU), Float64}}
├── Model clock: time = 0 seconds, iteration = 0
├── Next time step (TimeStepWizard{Float64, typeof(Oceananigans.Utils.cell_advection_timescale), typeof(Oceananigans.Simulations.infinite_diffusion_timescale)}): 10 ms
├── Iteration interval: 50
├── Stop criteria: Any[Oceananigans.Simulations.iteration_limit_exceeded, Oceananigans.Simulations.stop_time_exceeded, Oceananigans.Simulations.wall_time_limit_exceeded]
├── Run time: 0 seconds, wall time limit: Inf
├── Stop time: 40 seconds, stop iteration: Inf
├── Diagnostics: typename(OrderedCollections.OrderedDict) with 1 entry:
│   └── nan_checker => typename(NaNChecker)
└── Output writers: typename(OrderedCollections.OrderedDict) with no entries

### Output

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

u, v, w = model.velocities # unpack velocity Fields
b = model.tracers.b        # unpack buoyancy Field

# total flow speed
s = ComputedField(sqrt(u^2 + w^2))

# y-component of vorticity
ζ = ComputedField(∂z(u) - ∂x(w))

outputs = (s = s, b = b, ζ = ζ)

We create a JLD2OutputWriter that saves the speed, and the vorticity. We then add the JLD2OutputWriter to the simulation.

saved_output_prefix = "horizontal_convection"
saved_output_filename = saved_output_prefix * ".jld2"

simulation.output_writers[:fields] = JLD2OutputWriter(model, outputs,
field_slicer = nothing,
schedule = TimeInterval(0.5),
prefix = saved_output_prefix,
force = true)
┌ Warning: Cannot serialize timeseries/b/serialized/boundary_conditions as it contains functions. Will replace with missing. Function boundary conditions must be restored manually.
└ @ Oceananigans.OutputWriters ~/builds/tartarus-12/clima/oceananigans/src/OutputWriters/output_writer_utils.jl:65

Ready to press the big red button:

run!(simulation)
[ Info: Updating model auxiliary state before the first time step...
[ Info:     ... updated in 1.324 ms.
[ Info: Executing first time step...
i:     50, sim time:  0.500, wall time: 2.154 minutes, Δt:  0.0100, CFL: 3.12e-03
i:    100, sim time:  1.096, wall time: 2.204 minutes, Δt:  0.0120, CFL: 1.39e-02
i:    150, sim time:  1.802, wall time: 2.247 minutes, Δt:  0.0144, CFL: 3.65e-02
i:    200, sim time:  2.656, wall time: 2.291 minutes, Δt:  0.0173, CFL: 7.98e-02
i:    250, sim time:  3.666, wall time: 2.335 minutes, Δt:  0.0207, CFL: 1.63e-01
i:    300, sim time:  4.873, wall time: 2.379 minutes, Δt:  0.0249, CFL: 4.66e-01
i:    350, sim time:  6.328, wall time: 2.423 minutes, Δt:  0.0299, CFL: 1.27e+00
i:    400, sim time:  7.193, wall time: 2.466 minutes, Δt:  0.0176, CFL: 7.62e-01
i:    450, sim time:  8.052, wall time: 2.513 minutes, Δt:  0.0173, CFL: 7.74e-01
i:    500, sim time:  8.885, wall time: 2.556 minutes, Δt:  0.0168, CFL: 7.36e-01
i:    550, sim time:  9.722, wall time: 2.605 minutes, Δt:  0.0171, CFL: 7.58e-01
i:    600, sim time:  10.551, wall time: 2.653 minutes, Δt:  0.0169, CFL: 7.57e-01
i:    650, sim time:  11.385, wall time: 2.698 minutes, Δt:  0.0167, CFL: 7.07e-01
i:    700, sim time:  12.249, wall time: 2.742 minutes, Δt:  0.0178, CFL: 6.66e-01
i:    750, sim time:  13.220, wall time: 2.786 minutes, Δt:  0.0200, CFL: 7.79e-01
i:    800, sim time:  14.173, wall time: 2.836 minutes, Δt:  0.0193, CFL: 7.22e-01
i:    850, sim time:  15.160, wall time: 2.880 minutes, Δt:  0.0200, CFL: 7.00e-01
i:    900, sim time:  16.193, wall time: 2.927 minutes, Δt:  0.0214, CFL: 7.55e-01
i:    950, sim time:  17.234, wall time: 2.974 minutes, Δt:  0.0213, CFL: 7.51e-01
i:   1000, sim time:  18.277, wall time: 3.022 minutes, Δt:  0.0213, CFL: 7.30e-01
i:   1050, sim time:  19.350, wall time: 3.067 minutes, Δt:  0.0219, CFL: 7.19e-01
i:   1100, sim time:  20.479, wall time: 3.115 minutes, Δt:  0.0228, CFL: 7.15e-01
i:   1150, sim time:  21.668, wall time: 3.161 minutes, Δt:  0.0239, CFL: 7.66e-01
i:   1200, sim time:  22.805, wall time: 3.205 minutes, Δt:  0.0234, CFL: 8.39e-01
i:   1250, sim time:  23.835, wall time: 3.250 minutes, Δt:  0.0209, CFL: 6.79e-01
i:   1300, sim time:  24.962, wall time: 3.299 minutes, Δt:  0.0231, CFL: 7.10e-01
i:   1350, sim time:  26.146, wall time: 3.345 minutes, Δt:  0.0244, CFL: 6.72e-01
i:   1400, sim time:  27.491, wall time: 3.392 minutes, Δt:  0.0273, CFL: 8.04e-01
i:   1450, sim time:  28.729, wall time: 3.437 minutes, Δt:  0.0254, CFL: 9.03e-01
i:   1500, sim time:  29.775, wall time: 3.482 minutes, Δt:  0.0211, CFL: 8.76e-01
i:   1550, sim time:  30.663, wall time: 3.526 minutes, Δt:  0.0181, CFL: 6.07e-01
i:   1600, sim time:  31.717, wall time: 3.572 minutes, Δt:  0.0217, CFL: 8.21e-01
i:   1650, sim time:  32.678, wall time: 3.618 minutes, Δt:  0.0198, CFL: 7.27e-01
i:   1700, sim time:  33.664, wall time: 3.662 minutes, Δt:  0.0205, CFL: 7.03e-01
i:   1750, sim time:  34.718, wall time: 3.711 minutes, Δt:  0.0218, CFL: 7.12e-01
i:   1800, sim time:  35.845, wall time: 3.754 minutes, Δt:  0.0230, CFL: 7.08e-01
i:   1850, sim time:  37.024, wall time: 3.800 minutes, Δt:  0.0244, CFL: 7.29e-01
i:   1900, sim time:  38.275, wall time: 3.846 minutes, Δt:  0.0250, CFL: 6.86e-01
i:   1950, sim time:  39.527, wall time: 3.895 minutes, Δt:  0.0274, CFL: 8.38e-01
i:   1970, sim time:  40.000, wall time: 3.913 minutes, Δt:  0.0245, CFL: 7.56e-01
[ Info: Simulation is stopping. Model time 40 seconds has hit or exceeded simulation stop time 40 seconds.

## Load saved output, process, visualize

We animate the results by opening the JLD2 file, extracting data for the iterations we ended up saving at, and ploting the saved fields. From the saved buoyancy field we compute the buoyancy dissipation, $\chi = \kappa |\boldsymbol{\nabla} b|^2$, and plot that also.

To start we load the saved fields are FieldTimeSeries and prepare for animating the flow by creating coordinate arrays that each field lives on.

using JLD2, Plots
using Oceananigans
using Oceananigans.Fields
using Oceananigans.AbstractOperations: volume

saved_output_prefix = "horizontal_convection"
saved_output_filename = saved_output_prefix * ".jld2"

# Open the file with our data
file = jldopen(saved_output_filename)

κ = file["closure/κ/b"]

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

times = b_timeseries.times

# Coordinate arrays
xs, ys, zs = nodes(s_timeseries)
xb, yb, zb = nodes(b_timeseries)
xζ, yζ, zζ = nodes(ζ_timeseries)
WARNING: using Plots.grid in module Main conflicts with an existing identifier.

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

anim = @animate for i in 1:length(times)

# Load fields from FieldTimeSeries and compute χ
t = times[i]

s_snapshot = interior(s_timeseries[i])[:, 1, :]
ζ_snapshot = interior(ζ_timeseries[i])[:, 1, :]

b = b_timeseries[i]
χ = ComputedField(κ * (∂x(b)^2 + ∂z(b)^2))
compute!(χ)

b_snapshot = interior(b)[:, 1, :]
χ_snapshot = interior(χ)[:, 1, :]

# determine colorbar limits and contour levels
slim = 0.6
slevels = vcat(range(0, stop=slim, length=31), [slim])

blim = 0.6
blevels = vcat([-1], range(-blim, stop=blim, length=51), )

ζlim = 9
ζlevels = range(-ζlim, stop=ζlim, length=41)

χlim = 0.025
χlevels = range(0, stop=χlim, length=41)

@info @sprintf("Drawing frame %d:", i)

kwargs = (      xlims = (-Lx/2, Lx/2),
ylims = (-H, 0),
xlabel = "x / H",
ylabel = "z / H",
aspectratio = 1,
linewidth = 0)

s_plot = contourf(xs, zs, s_snapshot';
clims = (0, slim), levels = slevels,
color = :speed, kwargs...)
s_title = @sprintf("speed √[(u²+w²)/(b⋆H)] @ t=%1.2f", t)

b_plot = contourf(xb, zb, b_snapshot';
clims = (-blim, blim), levels = blevels,
color = :thermal, kwargs...)
b_title = @sprintf("buoyancy, b/b⋆ @ t=%1.2f", t)

ζ_plot = contourf(xζ, zζ, clamp.(ζ_snapshot', -ζlim, ζlim);
clims=(-ζlim, ζlim), levels = ζlevels,
color = :balance, kwargs...)
ζ_title = @sprintf("vorticity, (∂u/∂z - ∂w/∂x) √(H/b⋆) @ t=%1.2f", t)

χ_plot = contourf(xs, zs, χ_snapshot';
clims = (0, χlim), levels = χlevels,
color = :dense, kwargs...)
χ_title = @sprintf("buoyancy dissipation, κ|∇b|² √(H/b⋆⁵) @ t=%1.2f", t)

plot(s_plot, b_plot, ζ_plot, χ_plot,
size = (700, 1200),
layout = Plots.grid(4, 1),
title = [s_title b_title ζ_title χ_title])
end

At higher Rayleigh numbers the flow becomes much more vigorous. See, for example, an animation of the voricity of the fluid at $Ra = 10^{12}$ on vimeo.

### The Nusselt number

Often we are interested on how much the flow enhances mixing. This is quantified by the Nusselt number, which measures how much the flow enhances mixing compared if only diffusion was in operation. The Nusselt number is given by

$$$Nu = \frac{\langle \chi \rangle}{\langle \chi_{\rm diff} \rangle} \, ,$$$

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

$$$\kappa \nabla^2 b_{\rm diff} = 0 \, ,$$$

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

$$$b_{\rm diff}(x, z) = b_s(x) \frac{\cosh \left [2 \pi (H + z) / L_x \right ]}{\cosh(2 \pi H / L_x)} \, ,$$$

which implies $\langle \chi_{\rm diff} \rangle = \kappa b_*^2 \pi \tanh(2 \pi Η /Lx)$.

We use the loaded FieldTimeSeries to compute the Nusselt number from buoyancy and the volume average kinetic energy of the fluid.

First we compute the diffusive buoyancy dissipation, $\chi_{\rm diff}$ (which is just a scalar):

H  = file["grid/Lz"]
Lx = file["grid/Lx"]
κ  = file["closure/κ/b"]

χ_diff = κ * b★^2 * π * tanh(2π * H / Lx)

We then create two ReducedFields to store the volume integrals of the kinetic energy density and the buoyancy dissipation. We need the grid to do so; the grid can be recoverd from any FieldTimeSeries, e.g.,

grid = b_timeseries.grid

∫ⱽ_s² = ReducedField(Nothing, Nothing, Nothing, CPU(), grid, dims=(1, 2, 3))
∫ⱽ_mod²_∇b = ReducedField(Nothing, Nothing, Nothing, CPU(), grid, dims=(1, 2, 3))
Oceananigans.Fields.ReducedField located at (⋅, ⋅, ⋅)
├── architecture: CPU
└── grid: RegularRectilinearGrid{Float64, Bounded, Flat, Bounded}(Nx=128, Ny=1, Nz=64)

We recover the time from the saved FieldTimeSeries and construct two empty arrays to store the volume-averaged kinetic energy and the instantaneous Nusselt number,

t = b_timeseries.times

kinetic_energy, Nu = zeros(length(t)), zeros(length(t))

Now we can loop over the fields in the FieldTimeSeries, compute KineticEnergy and Nu, and plot.

for i = 1:length(t)
s = s_timeseries[i]
sum!(∫ⱽ_s², s^2 * volume)
kinetic_energy[i] = 0.5 * ∫ⱽ_s²[1, 1, 1]  / (Lx * H)

b = b_timeseries[i]
sum!(∫ⱽ_mod²_∇b, (∂x(b)^2 + ∂z(b)^2) * volume)
Nu[i] = (κ *  ∫ⱽ_mod²_∇b[1, 1, 1]) / χ_diff
end

p1 = plot(t, kinetic_energy,
xlabel = "time",
ylabel = "KE / (b⋆H)",
linewidth = 3,
legend = :none)

p2 = plot(t, Nu,
xlabel = "time",
ylabel = "Nu",
linewidth = 3,
legend = :none)

plot(p1, p2, layout = Plots.grid(2, 1))