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.

using Pkg
pkg"add Oceananigans, JLD2, CairoMakie"

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 = 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

\[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: 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, $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(; 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
├── tracers: b
├── closure: ScalarDiffusivity{ExplicitTimeDiscretization}(ν=0.000282843, κ=(b=0.000282843,))
├── buoyancy: BuoyancyTracer with -ĝ = ZDirection
└── 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$.

simulation = Simulation(model, Δt=1e-2, stop_time=40.0)
Simulation of NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── Next time step: 10 ms
├── Elapsed wall time: 0 seconds
├── Wall time per iteration: NaN years
├── Stop time: 40 seconds
├── Stop iteration : Inf
├── Wall time limit: Inf
├── 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
└── Diagnostics: OrderedDict with no entries

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.

wizard = TimeStepWizard(cfl=0.75, max_change=1.2, max_Δt=1e-1)

simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(50))
Callback of TimeStepWizard(cfl=0.75, max_Δt=0.1, min_Δt=0.0) on IterationInterval(50)

A progress messenger

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

progress(sim) = @printf("i: % 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, $\zeta$, and the buoyancy, $b$. Note that computed Fields take "AbstractOperations" on Fields as input:

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

# total flow speed
s = @at (Center, Center, Center) sqrt(u^2 + w^2)

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

We create a JLD2OutputWriter that saves the speed, and the vorticity. Because we want to post-process buoyancy and compute the buoyancy variance dissipation (which is proportional to `|\boldsymbol{\nabla} b|^2) we use thewith_halos = true`. This way, the halos for the fields are saved and thus when we load them as fields the will come with the proper boundary conditions.

We then add the JLD2OutputWriter to the simulation.

saved_output_filename = "horizontal_convection.jld2"

simulation.output_writers[:fields] = JLD2OutputWriter(model, (; s, b, ζ),
                                                      schedule = TimeInterval(0.5),
                                                      filename = saved_output_filename,
                                                      with_halos = true,
                                                      overwrite_existing = true)

Ready to press the big red button:

run!(simulation)
[ Info: Initializing simulation...
i:      0, sim time:  0.000, wall time:  0 seconds, Δt:  0.0120, advective CFL: 0.00e+00, diffusive CFL: 1.39e-02
[ Info:     ... simulation initialization complete (10.713 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (2.928 minutes).
i:     50, sim time:  0.596, wall time: 3.313 minutes, Δt:  0.0144, advective CFL: 6.08e-03, diffusive CFL: 1.67e-02
i:    100, sim time:  1.302, wall time: 3.535 minutes, Δt:  0.0173, advective CFL: 2.64e-02, diffusive CFL: 2.00e-02
i:    150, sim time:  2.156, wall time: 3.750 minutes, Δt:  0.0207, advective CFL: 6.95e-02, diffusive CFL: 2.40e-02
i:    200, sim time:  3.166, wall time: 3.910 minutes, Δt:  0.0249, advective CFL: 1.52e-01, diffusive CFL: 2.88e-02
i:    250, sim time:  4.373, wall time: 4.042 minutes, Δt:  0.0299, advective CFL: 3.62e-01, diffusive CFL: 3.46e-02
i:    300, sim time:  5.828, wall time: 4.195 minutes, Δt:  0.0190, advective CFL: 7.50e-01, diffusive CFL: 2.20e-02
i:    350, sim time:  6.746, wall time: 4.320 minutes, Δt:  0.0172, advective CFL: 7.50e-01, diffusive CFL: 1.99e-02
i:    400, sim time:  7.586, wall time: 4.448 minutes, Δt:  0.0169, advective CFL: 7.50e-01, diffusive CFL: 1.96e-02
i:    450, sim time:  8.423, wall time: 4.571 minutes, Δt:  0.0165, advective CFL: 7.50e-01, diffusive CFL: 1.91e-02
i:    500, sim time:  9.231, wall time: 4.715 minutes, Δt:  0.0177, advective CFL: 7.50e-01, diffusive CFL: 2.05e-02
i:    550, sim time:  10.088, wall time: 4.847 minutes, Δt:  0.0165, advective CFL: 7.50e-01, diffusive CFL: 1.91e-02
i:    600, sim time:  10.895, wall time: 4.973 minutes, Δt:  0.0169, advective CFL: 7.50e-01, diffusive CFL: 1.96e-02
i:    650, sim time:  11.703, wall time: 5.104 minutes, Δt:  0.0188, advective CFL: 7.50e-01, diffusive CFL: 2.17e-02
i:    700, sim time:  12.631, wall time: 5.228 minutes, Δt:  0.0190, advective CFL: 7.50e-01, diffusive CFL: 2.20e-02
i:    750, sim time:  13.538, wall time: 5.357 minutes, Δt:  0.0194, advective CFL: 7.50e-01, diffusive CFL: 2.25e-02
i:    800, sim time:  14.500, wall time: 5.488 minutes, Δt:  0.0202, advective CFL: 7.50e-01, diffusive CFL: 2.34e-02
i:    850, sim time:  15.500, wall time: 5.612 minutes, Δt:  0.0212, advective CFL: 7.50e-01, diffusive CFL: 2.45e-02
i:    900, sim time:  16.542, wall time: 5.756 minutes, Δt:  0.0211, advective CFL: 7.50e-01, diffusive CFL: 2.44e-02
i:    950, sim time:  17.563, wall time: 5.959 minutes, Δt:  0.0213, advective CFL: 7.50e-01, diffusive CFL: 2.47e-02
i:   1000, sim time:  18.606, wall time: 6.148 minutes, Δt:  0.0220, advective CFL: 7.50e-01, diffusive CFL: 2.55e-02
i:   1050, sim time:  19.698, wall time: 6.295 minutes, Δt:  0.0232, advective CFL: 7.50e-01, diffusive CFL: 2.68e-02
i:   1100, sim time:  20.824, wall time: 6.464 minutes, Δt:  0.0242, advective CFL: 7.50e-01, diffusive CFL: 2.80e-02
i:   1150, sim time:  22.000, wall time: 6.635 minutes, Δt:  0.0209, advective CFL: 7.50e-01, diffusive CFL: 2.42e-02
i:   1200, sim time:  23.042, wall time: 6.806 minutes, Δt:  0.0212, advective CFL: 7.50e-01, diffusive CFL: 2.46e-02
i:   1250, sim time:  24.064, wall time: 6.935 minutes, Δt:  0.0234, advective CFL: 7.50e-01, diffusive CFL: 2.71e-02
i:   1300, sim time:  25.210, wall time: 7.070 minutes, Δt:  0.0256, advective CFL: 7.50e-01, diffusive CFL: 2.97e-02
i:   1350, sim time:  26.461, wall time: 7.200 minutes, Δt:  0.0274, advective CFL: 7.50e-01, diffusive CFL: 3.18e-02
i:   1400, sim time:  27.774, wall time: 7.332 minutes, Δt:  0.0217, advective CFL: 7.50e-01, diffusive CFL: 2.52e-02
i:   1450, sim time:  28.848, wall time: 7.475 minutes, Δt:  0.0170, advective CFL: 7.50e-01, diffusive CFL: 1.97e-02
i:   1500, sim time:  29.687, wall time: 7.607 minutes, Δt:  0.0158, advective CFL: 7.50e-01, diffusive CFL: 1.83e-02
i:   1550, sim time:  30.473, wall time: 7.738 minutes, Δt:  0.0157, advective CFL: 7.50e-01, diffusive CFL: 1.82e-02
i:   1600, sim time:  31.220, wall time: 7.864 minutes, Δt:  0.0188, advective CFL: 7.07e-01, diffusive CFL: 2.18e-02
i:   1650, sim time:  32.151, wall time: 8.085 minutes, Δt:  0.0210, advective CFL: 7.50e-01, diffusive CFL: 2.43e-02
i:   1700, sim time:  33.189, wall time: 8.300 minutes, Δt:  0.0222, advective CFL: 7.50e-01, diffusive CFL: 2.57e-02
i:   1750, sim time:  34.267, wall time: 8.512 minutes, Δt:  0.0231, advective CFL: 7.50e-01, diffusive CFL: 2.68e-02
i:   1800, sim time:  35.393, wall time: 8.720 minutes, Δt:  0.0237, advective CFL: 7.50e-01, diffusive CFL: 2.74e-02
i:   1850, sim time:  36.524, wall time: 8.855 minutes, Δt:  0.0241, advective CFL: 7.50e-01, diffusive CFL: 2.79e-02
i:   1900, sim time:  37.717, wall time: 8.987 minutes, Δt:  0.0247, advective CFL: 7.50e-01, diffusive CFL: 2.86e-02
i:   1950, sim time:  38.920, wall time: 9.113 minutes, Δt:  0.0256, advective CFL: 7.50e-01, diffusive CFL: 2.97e-02
[ 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 CairoMakie
using Oceananigans
using Oceananigans.Fields
using Oceananigans.AbstractOperations: volume

saved_output_filename = "horizontal_convection.jld2"

# Open the file with our data
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
xc, yc, zc = nodes(b_timeseries[1])
xζ, yζ, zζ = nodes(ζ_timeseries[1])

χ_timeseries = deepcopy(b_timeseries)

for i in 1:length(times)
  bᵢ = b_timeseries[i]
  χ_timeseries[i] .= @at (Center, Center, Center) κ * (∂x(bᵢ)^2 + ∂z(bᵢ)^2)
end

Now we're ready to animate using Makie.

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

n = Observable(1)

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

sₙ = @lift interior(s_timeseries[$n], :, 1, :)
ζₙ = @lift interior(ζ_timeseries[$n], :, 1, :)
bₙ = @lift interior(b_timeseries[$n], :, 1, :)
χₙ = @lift interior(χ_timeseries[$n], :, 1, :)

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(resolution = (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];
            title = L"vorticity, $(∂u/∂z - ∂w/∂x) \, (L_x / b_*)^{1/2}$", axis_kwargs...)

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

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

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

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

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

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

And, finally, we record a movie.

frames = 1:length(times)

record(fig, "horizontal_convection.mp4", frames, framerate=8) do i
    msg = string("Plotting frame ", i, " of ", frames[end])
    print(msg * " \r")
    n[] = i
end
Plotting frame 1 of 81 
Plotting frame 2 of 81 
Plotting frame 3 of 81 
Plotting frame 4 of 81 
Plotting frame 5 of 81 
Plotting frame 6 of 81 
Plotting frame 7 of 81 
Plotting frame 8 of 81 
Plotting frame 9 of 81 
Plotting frame 10 of 81 
Plotting frame 11 of 81 
Plotting frame 12 of 81 
Plotting frame 13 of 81 
Plotting frame 14 of 81 
Plotting frame 15 of 81 
Plotting frame 16 of 81 
Plotting frame 17 of 81 
Plotting frame 18 of 81 
Plotting frame 19 of 81 
Plotting frame 20 of 81 
Plotting frame 21 of 81 
Plotting frame 22 of 81 
Plotting frame 23 of 81 
Plotting frame 24 of 81 
Plotting frame 25 of 81 
Plotting frame 26 of 81 
Plotting frame 27 of 81 
Plotting frame 28 of 81 
Plotting frame 29 of 81 
Plotting frame 30 of 81 
Plotting frame 31 of 81 
Plotting frame 32 of 81 
Plotting frame 33 of 81 
Plotting frame 34 of 81 
Plotting frame 35 of 81 
Plotting frame 36 of 81 
Plotting frame 37 of 81 
Plotting frame 38 of 81 
Plotting frame 39 of 81 
Plotting frame 40 of 81 
Plotting frame 41 of 81 
Plotting frame 42 of 81 
Plotting frame 43 of 81 
Plotting frame 44 of 81 
Plotting frame 45 of 81 
Plotting frame 46 of 81 
Plotting frame 47 of 81 
Plotting frame 48 of 81 
Plotting frame 49 of 81 
Plotting frame 50 of 81 
Plotting frame 51 of 81 
Plotting frame 52 of 81 
Plotting frame 53 of 81 
Plotting frame 54 of 81 
Plotting frame 55 of 81 
Plotting frame 56 of 81 
Plotting frame 57 of 81 
Plotting frame 58 of 81 
Plotting frame 59 of 81 
Plotting frame 60 of 81 
Plotting frame 61 of 81 
Plotting frame 62 of 81 
Plotting frame 63 of 81 
Plotting frame 64 of 81 
Plotting frame 65 of 81 
Plotting frame 66 of 81 
Plotting frame 67 of 81 
Plotting frame 68 of 81 
Plotting frame 69 of 81 
Plotting frame 70 of 81 
Plotting frame 71 of 81 
Plotting frame 72 of 81 
Plotting frame 73 of 81 
Plotting frame 74 of 81 
Plotting frame 75 of 81 
Plotting frame 76 of 81 
Plotting frame 77 of 81 
Plotting frame 78 of 81 
Plotting frame 79 of 81 
Plotting frame 80 of 81 
Plotting frame 81 of 81

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)} \, ,\]

where $b_s(x)$ is the surface boundary condition. The diffusive solution implies $\langle \chi_{\rm diff} \rangle = \kappa b_*^2 \pi \tanh(2 \pi Η /Lx) / (L_x H)$.

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):

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

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 kinetic energy and $Nu$, and plot. We make use of Integral to compute the volume integral of fields over our domain.

for i = 1:length(t)
    ke = Field(Integral(1/2 * s_timeseries[i]^2 / (Lx * H)))
    compute!(ke)
    kinetic_energy[i] = ke[1, 1, 1]

    χ = Field(Integral(χ_timeseries[i] / (Lx * H)))
    compute!(χ)

    Nu[i] = χ[1, 1, 1] / χ_diff
end

fig = Figure(resolution = (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, t, kinetic_energy; linewidth = 3)

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

fig

This page was generated using Literate.jl.