Eady turbulence example

In this example, we initialize a random velocity field and observe its viscous, turbulent decay in a two-dimensional domain. This example demonstrates:

  • How to use a tuple of turbulence closures
  • How to use biharmonic diffusivity
  • How to implement background velocity and tracer distributions
  • How to use ComputedFields for output

The Eady problem

The "Eady problem" simulates the baroclinic instability problem proposed by Eric Eady in the classic paper "Long waves and cyclone waves," Tellus (1949). The Eady problem is a simple, canonical model for the generation of mid-latitude atmospheric storms and the ocean eddies that permeate the world sea.

In the Eady problem, baroclinic motion and turublence is generated by the interaction between turbulent motions and a stationary, geostrophically-balanced basic state that is unstable to baroclinic instability. In this example, the baroclinic generation of turbulence due to extraction of energy from the geostrophic basic state is balanced by a bottom boundary condition that extracts momentum from turbulent motions and serves as a crude model for the drag associated with an unresolved and small-scale turbulent bottom boundary layer.

The geostrophic basic state

The geostrophic basic state in the Eady problem is represented by the streamfunction,

$ ψ(y, z) = - α y (z + L_z) \, ,$

where $α$ is the geostrophic shear and $L_z$ is the depth of the domain. The background buoyancy includes both the geostrophic flow component, $f ∂_z ψ$, where $f$ is the Coriolis parameter, and a background stable stratification component, $N^2 z$, where $N$ is the buoyancy frequency:

$ B(y, z) = f ∂_z ψ + N^2 z = - α f y + N^2 z \, .$

The background velocity field is related to the geostrophic streamfunction via $ U = - ∂_y ψ$ such that

$ U(z) = α (z + L_z) \, .$

Boundary conditions

All fields are periodic in the horizontal directions. We use "insulating", or zero-flux boundary conditions on the buoyancy perturbation at the top and bottom. We thus implicitly assume that the background vertical density gradient, $N^2 z$, is maintained by a process external to our simulation. We use free-slip, or zero-flux boundary conditions on $u$ and $v$ at the surface where $z=0$. At the bottom, we impose a momentum flux that extracts momentum and energy from the flow.

Bottom boundary condition: quadratic bottom drag

We model the effects of a turbulent bottom boundary layer on the eddy momentum budget with quadratic bottom drag. A quadratic cottom drag is introduced by imposing a vertical flux of horizontal momentum that removes momentum from the layer immediately above: in other words, the flux is negative (downwards) when the velocity at the bottom boundary is positive, and positive (upwards) with the velocity at the bottom boundary is negative. This drag term is "quadratic" because the rate at which momentum is removed is proportional to $\bm{u}_h \sqrt{u^2 + v^2}$, where $\bm{u}_h = u \bm{\hat x} + v \bm{\hat{v}$ is the horizontal velocity.

The $x$-component of the quadratic bottom drag is thus

$ τ{xz}(z=Lz) = - cᴰ u \sqrt{u^2 + v^2} \, , $

while the $y$-component is

$ τ{yz}(z=Lz) = - cᴰ v \sqrt{u^2 + v^2} \, , $

where $c^D$ is a dimensionless drag coefficient and $τ_{xz}(z=L_z)$ and $τ_{yz}(z=L_z)$ denote the flux of $u$ and $v$ momentum at $z = L_z$, the bottom of the domain.

Vertical and horizontal viscosity and diffusivity

Vertical and horizontal viscosties and diffusivities are required to stabilize the Eady problem and can be idealized as modeling the effect of turbulent mixing below the grid scale. For both tracers and velocities we use a Laplacian vertical diffusivity $κ_z ∂_z^2 c$ and a biharmonic horizontal diffusivity $ϰ_h (∂_x^4 + ∂_y^4) c$.

Eady problem summary and parameters

To summarize, the Eady problem parameters along with the values we use in this example are

Parameter nameDescriptionValueUnits
$ f $Coriolis parameter$ 10^{-4} $$ \mathrm{s^{-1}} $
$ N $Buoyancy frequency (square root of $\partial_z B$)$ 10^{-3} $$ \mathrm{s^{-1}} $
$ \alpha $Background vertical shear $\partial_z U$$ 10^{-3} $$ \mathrm{s^{-1}} $
$ c^D $Bottom quadratic drag coefficient$ 10^{-4} $none
$ κ_z $Laplacian vertical diffusivity$ 10^{-2} $$ \mathrm{m^2 s^{-1}} $
$ \varkappa_h $Biharmonic horizontal diffusivity$ 10^{-2} \times \Delta x^4 / \mathrm{day} $$ \mathrm{m^4 s^{-1}} $

We start off by importing Oceananigans, some convenient aliases for dimensions, and a function that generates a pretty string from a number that represents 'time' in seconds:

using Oceananigans, Printf
using Oceananigans.Utils: hour, day, prettytime

The grid

We use a three-dimensional grid with a depth of 4000 m and a horizontal extent of 1000 km, appropriate for mesoscale ocean dynamics with characteristic scales of 50-200 km.

grid = RegularCartesianGrid(size=(48, 48, 16), x=(0, 1e6), y=(0, 1e6), z=(-4e3, 0))
RegularCartesianGrid{Float64, Periodic, Periodic, Bounded}
                   domain: x ∈ [0.0, 1.0e6], y ∈ [0.0, 1.0e6], z ∈ [-4000.0, 0.0]
                 topology: (Periodic, Periodic, Bounded)
  resolution (Nx, Ny, Nz): (48, 48, 16)
   halo size (Hx, Hy, Hz): (1, 1, 1)
grid spacing (Δx, Δy, Δz): (20833.333333333332, 20833.333333333332, 250.0)

Rotation

The classical Eady problem is posed on an $f$-plane. We use a Coriolis parameter typical to mid-latitudes on Earth,

coriolis = FPlane(f=1e-4) # [s⁻¹]
FPlane{Float64}: f = 1.00e-04

The background flow

We build a NamedTuple of parameters that describe the background flow,

background_parameters = ( α = 10 * coriolis.f, # s⁻¹, geostrophic shear
                          f = coriolis.f,      # s⁻¹, Coriolis parameter
                          N = 1e-3,            # s⁻¹, buoyancy frequency
                         Lz = grid.Lz)         # m, ocean depth
(α = 0.001, f = 0.0001, N = 0.001, Lz = 4000.0)

and then construct the background fields $U$ and $B$

using Oceananigans.Fields: BackgroundField

# Background fields are defined via functions of x, y, z, t, and optional parameters
U(x, y, z, t, p) = + p.α * (z + p.Lz)
B(x, y, z, t, p) = - p.α * p.f * y + p.N^2 * z

U_field = BackgroundField(U, parameters=background_parameters)
B_field = BackgroundField(B, parameters=background_parameters)
BackgroundField{typeof(Main.ex-eady_turbulence.B), NamedTuple{(:α, :f, :N, :Lz),NTuple{4,Float64}}}
├── func: B
└── parameters: (α = 0.001, f = 0.0001, N = 0.001, Lz = 4000.0)

Boundary conditions

These boundary conditions prescribe a quadratic drag at the bottom as a flux condition. We also fix the surface and bottom buoyancy to enforce a buoyancy gradient N^2.

drag_coefficient = 1e-4

@inline drag_u(u, v, cᴰ) = - cᴰ * sqrt(u^2 + v^2) * u
@inline drag_v(u, v, cᴰ) = - cᴰ * sqrt(u^2 + v^2) * v

@inline bottom_drag_u(i, j, grid, clock, f, cᴰ) = @inbounds drag_u(f.u[i, j, 1], f.v[i, j, 1], cᴰ)
@inline bottom_drag_v(i, j, grid, clock, f, cᴰ) = @inbounds drag_v(f.u[i, j, 1], f.v[i, j, 1], cᴰ)

drag_bc_u = BoundaryCondition(Flux, bottom_drag_u, discrete_form=true, parameters=drag_coefficient)
drag_bc_v = BoundaryCondition(Flux, bottom_drag_v, discrete_form=true, parameters=drag_coefficient)

u_bcs = UVelocityBoundaryConditions(grid, bottom = drag_bc_u)
v_bcs = VVelocityBoundaryConditions(grid, bottom = drag_bc_v)
Oceananigans.FieldBoundaryConditions (NamedTuple{(:x, :y, :z)}), with boundary conditions
├── x: CoordinateBoundaryConditions{BoundaryCondition{Oceananigans.BoundaryConditions.Periodic,Nothing},BoundaryCondition{Oceananigans.BoundaryConditions.Periodic,Nothing}}
├── y: CoordinateBoundaryConditions{BoundaryCondition{Oceananigans.BoundaryConditions.Periodic,Nothing},BoundaryCondition{Oceananigans.BoundaryConditions.Periodic,Nothing}}
└── z: CoordinateBoundaryConditions{BoundaryCondition{Flux,Oceananigans.BoundaryConditions.ParameterizedDiscreteBoundaryFunction{typeof(Main.ex-eady_turbulence.bottom_drag_v),Float64}},BoundaryCondition{Flux,Nothing}}

Turbulence closures

We use a horizontal biharmonic diffusivity and a Laplacian vertical diffusivity to dissipate energy in the Eady problem. To use both of these closures at the same time, we set the keyword argument closure to a tuple of two closures.

κ₂z = 1e-2 # [m² s⁻¹] Laplacian vertical viscosity and diffusivity
κ₄h = 1e-1 / day * grid.Δx^4 # [m⁴ s⁻¹] biharmonic horizontal viscosity and diffusivity

Laplacian_vertical_diffusivity = AnisotropicDiffusivity(νh=0, κh=0, νz=κ₂z, κz=κ₂z)
biharmonic_horizontal_diffusivity = AnisotropicBiharmonicDiffusivity(νh=κ₄h, κh=κ₄h)
AnisotropicBiharmonicDiffusivity{Float64,Float64,Float64,Float64}(2.1803253690129166e11, 2.1803253690129166e11, 0.0, 2.1803253690129166e11, 2.1803253690129166e11, 0.0)

Model instantiation

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

using Oceananigans.Advection: WENO5

model = IncompressibleModel(
           architecture = CPU(),
                   grid = grid,
              advection = WENO5(),
            timestepper = :RungeKutta3,
               coriolis = coriolis,
                tracers = :b,
               buoyancy = BuoyancyTracer(),
      background_fields = (b=B_field, u=U_field),
                closure = (Laplacian_vertical_diffusivity, biharmonic_horizontal_diffusivity),
    boundary_conditions = (u=u_bcs, v=v_bcs)
)
IncompressibleModel{CPU, Float64}(time = 0 seconds, iteration = 0) 
├── grid: RegularCartesianGrid{Float64, Periodic, Periodic, Bounded}(Nx=48, Ny=48, Nz=16)
├── tracers: (:b,)
├── closure: Tuple{AnisotropicDiffusivity{Float64,Float64,Float64,NamedTuple{(:b,),Tuple{Float64}},NamedTuple{(:b,),Tuple{Float64}},NamedTuple{(:b,),Tuple{Float64}}},AnisotropicBiharmonicDiffusivity{Float64,NamedTuple{(:b,),Tuple{Float64}},NamedTuple{(:b,),Tuple{Float64}},NamedTuple{(:b,),Tuple{Float64}}}}
├── buoyancy: BuoyancyTracer
└── coriolis: FPlane{Float64}

Initial conditions

We use non-trivial initial conditions consisting of an array of vortices superposed with large-amplitude noise to (hopefully) stimulate the rapid growth of baroclinic instability.

# A noise function, damped at the top and bottom
Ξ(z) = randn() * z/grid.Lz * (z/grid.Lz + 1)

# Scales for the initial velocity and buoyancy
Ũ = 1e-1 * background_parameters.α * grid.Lz
B̃ = 1e-2 * background_parameters.α * coriolis.f

uᵢ(x, y, z) = Ũ * Ξ(z)
vᵢ(x, y, z) = Ũ * Ξ(z)
bᵢ(x, y, z) = B̃ * Ξ(z)

set!(model, u=uᵢ, v=vᵢ, b=bᵢ)

We subtract off any residual mean velocity to avoid exciting domain-scale inertial oscillations. We use a sum over the entire parent arrays or data to ensure this operation is efficient on the GPU (set architecture = GPU() in IncompressibleModel constructor to run this problem on the GPU if one is available).

Ū = sum(model.velocities.u.data.parent) / (grid.Nx * grid.Ny * grid.Nz)
V̄ = sum(model.velocities.v.data.parent) / (grid.Nx * grid.Ny * grid.Nz)

model.velocities.u.data.parent .-= Ū
model.velocities.v.data.parent .-= V̄

Simulation set-up

We set up a simulation that runs for 10 days with a JLD2OutputWriter that saves the vertical vorticity and divergence every 2 hours.

The TimeStepWizard

The TimeStepWizard manages the time-step adaptively, keeping the CFL close to a desired value.

# Calculate absolute limit on time-step using diffusivities and
# background velocity.
Ū = background_parameters.α * grid.Lz

max_Δt = min(grid.Δx / Ū, grid.Δx^4 / κ₄h, grid.Δz^2 / κ₂z)

wizard = TimeStepWizard(cfl=1.0, Δt=0.1*max_Δt, max_change=1.1, max_Δt=max_Δt)
TimeStepWizard{Float64}(1.0, Inf, 1.1, 0.5, 5208.333333333333, 520.8333333333334)

A progress messenger

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

using Oceananigans.Diagnostics: AdvectiveCFL

CFL = AdvectiveCFL(wizard)

start_time = time_ns()

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

Build the simulation

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

simulation = Simulation(model, Δt = wizard, iteration_interval = 20,
                                                     stop_time = 10day,
                                                      progress = progress)
Simulation{IncompressibleModel{CPU, Float64}}
├── Model clock: time = 0 seconds, iteration = 0 
├── Next time step (TimeStepWizard{Float64}): 8.681 minutes 
├── Iteration interval: 20
├── 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: 10 days, stop iteration: Inf
├── Diagnostics: OrderedCollections.OrderedDict with no entries
└── Output writers: OrderedCollections.OrderedDict with no entries

Output

To visualize the baroclinic turbulence ensuing in the Eady problem, we use ComputedFields to diagnose and output vertical vorticity and divergence. Note that ComputedFields take "AbstractOperations" on Fields as input:

using Oceananigans.AbstractOperations
using Oceananigans.Fields: ComputedField

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

# Vertical vorticity [s⁻¹]
ζ = ComputedField(∂x(v) - ∂y(u))

# Horizontal divergence, or ∂x(u) + ∂y(v) [s⁻¹]
δ = ComputedField(-∂z(w))

With the vertical vorticity, ζ, and the horizontal divergence, δ in hand, we create a JLD2OutputWriter that saves ζ and δ and add it to simulation.

using Oceananigans.OutputWriters: JLD2OutputWriter, TimeInterval

simulation.output_writers[:fields] = JLD2OutputWriter(model, (ζ=ζ, δ=δ),
                                                      schedule = TimeInterval(2hour),
                                                        prefix = "eady_turbulence",
                                                         force = true)

All that's left is to press the big red button:

run!(simulation)
i:     20, sim time: 2.894 hours, wall time: 1.292 minutes, Δt: 8.681 minutes, CFL: 2.13e-02
i:     40, sim time: 6.076 hours, wall time: 1.343 minutes, Δt: 9.549 minutes, CFL: 1.99e-02
i:     60, sim time: 9.578 hours, wall time: 1.393 minutes, Δt: 10.503 minutes, CFL: 1.59e-02
i:     80, sim time: 13.429 hours, wall time: 1.443 minutes, Δt: 11.554 minutes, CFL: 1.78e-02
i:    100, sim time: 17.665 hours, wall time: 1.493 minutes, Δt: 12.709 minutes, CFL: 2.28e-02
i:    120, sim time: 22.325 hours, wall time: 1.544 minutes, Δt: 13.980 minutes, CFL: 2.69e-02
i:    140, sim time: 1.144 days, wall time: 1.593 minutes, Δt: 15.378 minutes, CFL: 3.80e-02
i:    160, sim time: 1.379 days, wall time: 1.643 minutes, Δt: 16.916 minutes, CFL: 5.85e-02
i:    180, sim time: 1.637 days, wall time: 1.693 minutes, Δt: 18.608 minutes, CFL: 8.25e-02
i:    200, sim time: 1.921 days, wall time: 1.743 minutes, Δt: 20.468 minutes, CFL: 9.11e-02
i:    220, sim time: 2.234 days, wall time: 1.793 minutes, Δt: 22.515 minutes, CFL: 1.75e-01
i:    240, sim time: 2.578 days, wall time: 1.843 minutes, Δt: 24.767 minutes, CFL: 2.62e-01
i:    260, sim time: 2.957 days, wall time: 1.893 minutes, Δt: 27.243 minutes, CFL: 4.44e-01
i:    280, sim time: 3.373 days, wall time: 1.943 minutes, Δt: 29.968 minutes, CFL: 6.97e-01
i:    300, sim time: 3.831 days, wall time: 1.993 minutes, Δt: 32.964 minutes, CFL: 1.16e+00
i:    320, sim time: 4.225 days, wall time: 2.043 minutes, Δt: 28.412 minutes, CFL: 1.31e+00
i:    340, sim time: 4.526 days, wall time: 2.093 minutes, Δt: 21.654 minutes, CFL: 7.86e-01
i:    360, sim time: 4.857 days, wall time: 2.143 minutes, Δt: 23.820 minutes, CFL: 8.64e-01
i:    380, sim time: 5.221 days, wall time: 2.193 minutes, Δt: 26.202 minutes, CFL: 1.17e+00
i:    400, sim time: 5.531 days, wall time: 2.243 minutes, Δt: 22.349 minutes, CFL: 7.28e-01
i:    420, sim time: 5.873 days, wall time: 2.293 minutes, Δt: 24.584 minutes, CFL: 1.03e+00
i:    440, sim time: 6.205 days, wall time: 2.343 minutes, Δt: 23.923 minutes, CFL: 1.05e+00
i:    460, sim time: 6.521 days, wall time: 2.394 minutes, Δt: 22.753 minutes, CFL: 1.03e+00
i:    480, sim time: 6.829 days, wall time: 2.444 minutes, Δt: 22.166 minutes, CFL: 8.01e-01
i:    500, sim time: 7.167 days, wall time: 2.494 minutes, Δt: 24.383 minutes, CFL: 8.36e-01
i:    520, sim time: 7.540 days, wall time: 2.544 minutes, Δt: 26.821 minutes, CFL: 1.13e+00
i:    540, sim time: 7.869 days, wall time: 2.594 minutes, Δt: 23.669 minutes, CFL: 9.42e-01
i:    560, sim time: 8.218 days, wall time: 2.644 minutes, Δt: 25.134 minutes, CFL: 1.10e+00
i:    580, sim time: 8.536 days, wall time: 2.694 minutes, Δt: 22.924 minutes, CFL: 1.02e+00
i:    600, sim time: 8.847 days, wall time: 2.744 minutes, Δt: 22.365 minutes, CFL: 1.07e+00
i:    620, sim time: 9.138 days, wall time: 2.793 minutes, Δt: 20.953 minutes, CFL: 1.09e+00
i:    640, sim time: 9.405 days, wall time: 2.843 minutes, Δt: 19.209 minutes, CFL: 1.29e+00
i:    660, sim time: 9.612 days, wall time: 2.893 minutes, Δt: 14.941 minutes, CFL: 9.47e-01
i:    680, sim time: 9.831 days, wall time: 2.943 minutes, Δt: 15.778 minutes, CFL: 9.44e-01
i:    700, sim time: 10.063 days, wall time: 2.993 minutes, Δt: 16.719 minutes, CFL: 8.92e-01
[ Info: Simulation is stopping. Model time 10.063 days has hit or exceeded simulation stop time 10 days.

Visualizing Eady turbulence

We animate the results by opening the JLD2 file, extracting data for the iterations we ended up saving at, and ploting slices of the saved fields. We prepare for animating the flow by creating coordinate arrays, opening the file, building a vector of the iterations that we saved data at, and defining a function for computing colorbar limits:

using JLD2, Plots

using Oceananigans.Grids: nodes, x_domain, y_domain, z_domain # for nice domain limits

# Coordinate arrays
xζ, yζ, zζ = nodes(ζ)
xδ, yδ, zδ = nodes(δ)

# Open the file with our data
file = jldopen(simulation.output_writers[:fields].filepath)

# Extract a vector of iterations
iterations = parse.(Int, keys(file["timeseries/t"]))
121-element Array{Int64,1}:
   0
  14
  27
  40
  51
  63
  73
  83
  93
 102
   ⋮
 635
 642
 650
 658
 665
 673
 681
 688
 695

This utility is handy for calculating nice contour intervals:

function nice_divergent_levels(c, clim, nlevels=30)
    levels = range(-clim, stop=clim, length=nlevels)

    cmax = maximum(abs, c)
    if clim < cmax # add levels on either end
        levels = vcat([-cmax], range(-clim, stop=clim, length=nlevels), [cmax])
    end

    return levels
end
nice_divergent_levels (generic function with 2 methods)

Now we're ready to animate.

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

anim = @animate for (i, iter) in enumerate(iterations)

    # Load 3D fields from file
    t = file["timeseries/t/$iter"]
    R = file["timeseries/ζ/$iter"] ./ coriolis.f
    δ = file["timeseries/δ/$iter"]

    surface_R = R[:, :, grid.Nz]
    surface_δ = δ[:, :, grid.Nz]

    slice_R = R[:, 1, :]
    slice_δ = δ[:, 1, :]

    Rlim = 0.5 * maximum(abs, R) + 1e-9
    δlim = 0.5 * maximum(abs, δ) + 1e-9

    Rlevels = nice_divergent_levels(R, Rlim)
    δlevels = nice_divergent_levels(δ, δlim)

    @info @sprintf("Drawing frame %d from iteration %d: max(ζ̃ / f) = %.3f \n",
                   i, iter, maximum(abs, surface_R))

    R_xy = contourf(xζ, yζ, surface_R';
                    aspectratio = 1,
                      linewidth = 0,
                          color = :balance,
                         legend = false,
                          clims = (-Rlim, Rlim),
                         levels = Rlevels,
                          xlims = (0, grid.Lx),
                          ylims = (0, grid.Lx),
                         xlabel = "x (m)",
                         ylabel = "y (m)")

    δ_xy = contourf(xδ, yδ, surface_δ';
                    aspectratio = 1,
                      linewidth = 0,
                          color = :balance,
                         legend = false,
                          clims = (-δlim, δlim),
                         levels = δlevels,
                          xlims = (0, grid.Lx),
                          ylims = (0, grid.Lx),
                         xlabel = "x (m)",
                         ylabel = "y (m)")

    R_xz = contourf(xζ, zζ, slice_R';
                    aspectratio = grid.Lx / grid.Lz * 0.5,
                      linewidth = 0,
                          color = :balance,
                         legend = false,
                          clims = (-Rlim, Rlim),
                         levels = Rlevels,
                          xlims = (0, grid.Lx),
                          ylims = (-grid.Lz, 0),
                         xlabel = "x (m)",
                         ylabel = "z (m)")

    δ_xz = contourf(xδ, zδ, slice_δ';
                    aspectratio = grid.Lx / grid.Lz * 0.5,
                      linewidth = 0,
                          color = :balance,
                         legend = false,
                          clims = (-δlim, δlim),
                         levels = δlevels,
                          xlims = (0, grid.Lx),
                          ylims = (-grid.Lz, 0),
                         xlabel = "x (m)",
                         ylabel = "z (m)")

    plot(R_xy, δ_xy, R_xz, δ_xz,
           size = (1000, 800),
           link = :x,
         layout = Plots.grid(2, 2, heights=[0.5, 0.5, 0.2, 0.2]),
          title = [@sprintf("ζ/f(t=%s) (s⁻¹)", prettytime(t)) @sprintf("δ(t=%s) (s⁻¹)", prettytime(t)) "" ""])

    iter == iterations[end] && close(file)
end

This page was generated using Literate.jl.