# An unstable Bickley jet in Shallow Water model

This example uses Oceananigans.jl's ShallowWaterModel to simulate the evolution of an unstable, geostrophically balanced, Bickley jet The example is periodic in $x$ with flat bathymetry and uses the conservative formulation of the shallow water equations. The initial conditions superpose the Bickley jet with small-amplitude perturbations. See "The nonlinear evolution of barotropically unstable jets," J. Phys. Oceanogr. (2003) for more details on this problem.

The mass transport $(uh, vh)$ is the prognostic momentum variable in the conservative formulation of the shallow water equations, where $(u, v)$ are the horizontal velocity components and $h$ is the layer height.

## Install dependencies

First we make sure that we have all of the packages that are required to run the simulation.

using Pkg
pkg"add Oceananigans, NCDatasets, Plots, Printf, Polynomials"
using Oceananigans
using Oceananigans.Models: ShallowWaterModel

## Two-dimensional domain

The shallow water model is a two-dimensional model and thus the number of vertical points Nz must be set to one. Note that $L_z$ is the mean depth of the fluid.

Lx, Ly, Lz = 2π, 20, 1
Nx, Ny = 128, 128

grid = RegularRectilinearGrid(size = (Nx, Ny),
x = (0, Lx), y = (-Ly/2, Ly/2),
topology = (Periodic, Bounded, Flat))
RegularRectilinearGrid{Float64, Periodic, Bounded, Flat}
domain: x ∈ [0.0, 6.283185307179586], y ∈ [-10.0, 10.0], z ∈ [0.0, 0.0]
topology: (Periodic, Bounded, Flat)
resolution (Nx, Ny, Nz): (128, 128, 1)
halo size (Hx, Hy, Hz): (1, 1, 0)
grid spacing (Δx, Δy, Δz): (0.04908738521234052, 0.15625, 0.0)

## Physical parameters

We choose non-dimensional parameters

const U = 1.0         # Maximum jet velocity

f = 1           # Coriolis parameter
g = 9.8         # Gravitational acceleration
Δη = f * U / g  # Maximum free-surface deformation as dictated by geostrophy
0.1020408163265306

## Building a ShallowWaterModel

We build a ShallowWaterModel with the WENO5 advection scheme and 3rd-order Runge-Kutta time-stepping,

model = ShallowWaterModel(architecture = CPU(),
timestepper = :RungeKutta3,
grid = grid,
gravitational_acceleration = g,
coriolis = FPlane(f=f))
ShallowWaterModel{typename(CPU), Float64}(time = 0 seconds, iteration = 0)
├── grid: RegularRectilinearGrid{Float64, Periodic, Bounded, Flat}(Nx=128, Ny=128, Nz=1)
├── tracers: ()
└── coriolis: FPlane{Float64}

Use architecture = GPU() to run this problem on a GPU.

## Background state and perturbation

The background velocity $ū$ and free-surface $η̄$ correspond to a geostrophically balanced Bickely jet with maximum speed of $U$ and maximum free-surface deformation of $Δη$,

h̄(x, y, z) = Lz - Δη * tanh(y)
ū(x, y, z) = U * sech(y)^2
ū (generic function with 1 method)

The total height of the fluid is $h = L_z + \eta$. Linear stability theory predicts that for the parameters we consider here, the growth rate for the most unstable mode that fits our domain is approximately $0.139$.

The vorticity of the background state is

ω̄(x, y, z) = 2 * U * sech(y)^2 * tanh(y)
ω̄ (generic function with 1 method)

The initial conditions include a small-amplitude perturbation that decays away from the center of the jet.

small_amplitude = 1e-4

uⁱ(x, y, z) = ū(x, y, z) + small_amplitude * exp(-y^2) * randn()
uhⁱ(x, y, z) = uⁱ(x, y, z) * h̄(x, y, z)
uhⁱ (generic function with 1 method)

We first set a "clean" initial condition without noise for the purpose of discretely calculating the initial 'mean' vorticity,

ū̄h(x, y, z) = ū(x, y, z) * h̄(x, y, z)

set!(model, uh = ū̄h, h = h̄)

We next compute the initial vorticity and perturbation vorticity,

uh, vh, h = model.solution

# Build velocities
u = uh / h
v = vh / h

# Build and compute mean vorticity discretely
ω = ComputedField(∂x(v) - ∂y(u))
compute!(ω)

# Copy mean vorticity to a new field
ωⁱ = Field(Face, Face, Nothing, model.architecture, model.grid)
ωⁱ .= ω

# Use this new field to compute the perturbation vorticity
ω′ = ComputedField(ω - ωⁱ)
ComputedField located at (Face, Face, Center) of BinaryOperation at (Face, Face, Center)
├── data: OffsetArrays.OffsetArray{Float64, 3, Array{Float64, 3}}, size: (128, 129, 1)
├── grid: RegularRectilinearGrid{Float64, Periodic, Bounded, Flat}(Nx=128, Ny=128, Nz=1)
├── operand: BinaryOperation at (Face, Face, Center)
└── status: time=0.0

and finally set the "true" initial condition with noise,

set!(model, uh = uhⁱ)

## Running a Simulation

We pick the time-step so that we make sure we resolve the surface gravity waves, which propagate with speed of the order $\sqrt{g L_z}$. That is, with Δt = 1e-2 we ensure that $\sqrt{g L_z} Δt / Δx, \sqrt{g L_z} Δt / Δy < 0.7$.

simulation = Simulation(model, Δt = 1e-2, stop_time = 150)
Simulation{typename(ShallowWaterModel){typename(CPU), Float64}}
├── Model clock: time = 0 seconds, iteration = 0
├── Next time step (Float64): 10 ms
├── Iteration interval: 1
├── 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: 2.500 minutes, stop iteration: Inf
├── Diagnostics: typename(OrderedCollections.OrderedDict) with 1 entry:
│   └── nan_checker => typename(NaNChecker)
└── Output writers: typename(OrderedCollections.OrderedDict) with no entries

## Prepare output files

Define a function to compute the norm of the perturbation on the cross channel velocity. We obtain the norm function from LinearAlgebra.

using LinearAlgebra: norm

perturbation_norm(args...) = norm(v)
perturbation_norm (generic function with 1 method)

Build the output_writer for the two-dimensional fields to be output. Output every t = 1.0.

simulation.output_writers[:fields] = NetCDFOutputWriter(model, (; ω, ω′),
filepath = joinpath(@__DIR__, "shallow_water_Bickley_jet_fields.nc"),
schedule = TimeInterval(1),
mode = "c")
NetCDFOutputWriter scheduled on TimeInterval(1 second):
├── filepath: /var/lib/buildkite-agent/builds/tartarus-12/clima/oceananigans/docs/build/generated/shallow_water_Bickley_jet_fields.nc
├── dimensions: zC(1), zF(1), xC(128), yF(129), xF(128), yC(128), time(0)
├── 2 outputs: ["ω", "ω′"]
├── field slicer: FieldSlicer(:, :, :, with_halos=false)
└── array type: Array{Float32}

Build the output_writer for the growth rate, which is a scalar field. Output every time step.

simulation.output_writers[:growth] = NetCDFOutputWriter(model, (; perturbation_norm),
filepath = joinpath(@__DIR__, "shallow_water_Bickley_jet_perturbation_norm.nc"),
schedule = IterationInterval(1),
dimensions = (; perturbation_norm = ()),
mode = "c")
NetCDFOutputWriter scheduled on IterationInterval(1):
├── filepath: /var/lib/buildkite-agent/builds/tartarus-12/clima/oceananigans/docs/build/generated/shallow_water_Bickley_jet_perturbation_norm.nc
├── dimensions: zC(1), zF(1), xC(128), yF(129), xF(128), yC(128), time(0)
├── 1 outputs: ["perturbation_norm"]
├── field slicer: FieldSlicer(:, :, :, with_halos=false)
└── array type: Array{Float32}

And finally run the simulation.

run!(simulation)
[ Info: Updating model auxiliary state before the first time step...
[ Info:     ... updated in 452.855 μs.
[ Info: Executing first time step...
[ Info: Simulation is stopping. Model time 2.500 minutes has hit or exceeded simulation stop time 2.500 minutes.

## Visualize the results

using NCDatasets, Plots, Printf
WARNING: using Plots.grid in module Main conflicts with an existing identifier.

Define the coordinates for plotting.

x, y = xnodes(ω), ynodes(ω)

Define keyword arguments for plotting the contours.

kwargs = (
xlabel = "x",
ylabel = "y",
aspect = 1,
fill = true,
levels = 20,
linewidth = 0,
color = :balance,
colorbar = true,
ylim = (-Ly/2, Ly/2),
xlim = (0, Lx)
)

Read in the output_writer for the two-dimensional fields and then create an animation showing both the total and perturbation vorticities.

ds = NCDataset(simulation.output_writers[:fields].filepath, "r")

iterations = keys(ds["time"])

anim = @animate for (iter, t) in enumerate(ds["time"])
ω = ds["ω"][:, :, 1, iter]
ω′ = ds["ω′"][:, :, 1, iter]

ω′_max = maximum(abs, ω′)

plot_ω = contour(x, y, ω',
clim = (-1, 1),
title = @sprintf("Total vorticity, ω, at t = %.1f", t); kwargs...)

plot_ω′ = contour(x, y, ω′',
clim = (-ω′_max, ω′_max),
title = @sprintf("Perturbation vorticity, ω - ω̄, at t = %.1f", t); kwargs...)

plot(plot_ω, plot_ω′, layout = (1, 2), size = (800, 440))
end

close(ds)

mp4(anim, "shallow_water_Bickley_jet.mp4", fps=15)

Read in the output_writer for the scalar field (the norm of $v$-velocity).

ds2 = NCDataset(simulation.output_writers[:growth].filepath, "r")

iterations = keys(ds2["time"])

t = ds2["time"][:]
norm_v = ds2["perturbation_norm"][:]

close(ds2)

We import the fit function from Polynomials.jl to compute the best-fit slope of the perturbation norm on a logarithmic plot. This slope corresponds to the growth rate.

using Polynomials: fit

I = 6000:7000

degree = 1
linear_fit_polynomial = fit(t[I], log.(norm_v[I]), degree, var = :t)
-9.624761392416577 + 0.13758522902817455∙t

We can get the coefficient of the $n$-th power from the fitted polynomial by using n as an index, e.g.,

constant, slope = linear_fit_polynomial[0], linear_fit_polynomial[1]
(-9.624761392416577, 0.13758522902817455)

We then use the computed linear fit coefficients to construct the best fit and plot it together with the time-series for the perturbation norm for comparison.

best_fit = @. exp(constant + slope * t)

plot(t, norm_v,
yaxis = :log,
ylims = (1e-3, 30),
lw = 4,
label = "norm(v)",
xlabel = "time",
ylabel = "norm(v)",
title = "growth of perturbation norm",
legend = :bottomright)

plot!(t[I], 2 * best_fit[I], # factor 2 offsets fit from curve for better visualization
lw = 4,
label = "best fit")

The slope of the best-fit curve on a logarithmic scale approximates the rate at which instability grows in the simulation. Let's see how this compares with the theoretical growth rate.

println("Numerical growth rate is approximated to be ", round(slope, digits=3), ",\n",
"which is very close to the theoretical value of 0.139.")
Numerical growth rate is approximated to be 0.138,
which is very close to the theoretical value of 0.139.