Internal wave example

In this example, we initialize an internal wave packet in two-dimensions and watch it propagate. This example illustrates how to set up a two-dimensional model, set initial conditions, and how to use BackgroundFields.

Numerical, domain, and internal wave parameters

First, we pick a resolution and domain size. We use a two-dimensional domain that's periodic in $(x, y, z)$:

using Oceananigans

grid = RegularCartesianGrid(size=(128, 1, 128), x=(-π, π), y=(-π, π), z=(-π, π),
                            topology=(Periodic, Periodic, Periodic))

RegularCartesianGrid{Float64, Periodic, Periodic, Periodic}
                   domain: x ∈ [-3.141592653589793, 3.141592653589793], y ∈ [-3.141592653589793, 3.141592653589793], z ∈ [-3.141592653589793, 3.141592653589793]
                 topology: (Periodic, Periodic, Periodic)
  resolution (Nx, Ny, Nz): (128, 1, 128)
   halo size (Hx, Hy, Hz): (1, 1, 1)
grid spacing (Δx, Δy, Δz): (0.04908738521234052, 6.283185307179586, 0.04908738521234052)

Inertia-gravity waves propagate in fluids that are both (i) rotating, and (ii) density-stratified. We use Oceananigans' coriolis abstraction to implement the background rotation rate:

coriolis = FPlane(f=0.2)

FPlane{Float64}: f = 2.00e-01

On an FPlane, the domain is idealized as rotating at a constant rate with rotation period 2π/f. coriolis is passed to IncompressibleModel below. Our units are arbitrary.

We use Oceananigans' background_fields abstraction to define a background buoyancy field background_b(z) = N^2 * z, where z is the vertical coordinate and N is the "buoyancy frequency". This means that the modeled buoyany field in Oceananigans will be a perturbation away from the basic state background_b. We choose

N = 1 ## buoyancy frequency

and then construct the background buoyancy,

using Oceananigans.Fields: BackgroundField

# Background fields are functions of `x, y, z, t`, and optional parameters.
# Here we have one parameter, the buoyancy frequency
B_func(x, y, z, t, N) = N^2 * z

B = BackgroundField(B_func, parameters=N)

BackgroundField{typeof(Main.ex-internal_wave.B_func), Int64}
├── func: B_func
└── parameters: 1

We are now ready to instantiate our model. We pass grid, coriolis, and background_b to the IncompressibleModel constructor. In addition, we add a small amount of IsotropicDiffusivity to keep the model stable, during time-stepping, and specify our model to use a single tracer called b that we identify as buoyancy by setting buoyancy=BuoyancyTracer().

model = IncompressibleModel(
                 grid = grid,
          timestepper = :RungeKutta3,
              closure = IsotropicDiffusivity(ν=1e-6, κ=1e-6),
             coriolis = coriolis,
              tracers = :b,
    background_fields = (b=B,), # `background_fields` is a `NamedTuple`
             buoyancy = BuoyancyTracer()
)

IncompressibleModel{CPU, Float64}(time = 0 seconds, iteration = 0) 
├── grid: RegularCartesianGrid{Float64, Periodic, Periodic, Periodic}(Nx=128, Ny=1, Nz=128)
├── tracers: (:b,)
├── closure: IsotropicDiffusivity{Float64,NamedTuple{(:b,),Tuple{Float64}}}
├── buoyancy: BuoyancyTracer
└── coriolis: FPlane{Float64}

A Gaussian wavepacket

Next, we set up an initial condition that excites an internal wave that propates through our rotating, stratified fluid. This internal wave has the pressure field

$ p(x, y, z, t) = a(x, z) \, \cos(kx + mz - ω t) $.

where $m$ is the vertical wavenumber, $k$ is the horizontal wavenumber, $ω$ is the wave frequncy, and $a(x, z)$ is a Gaussian envelope. The internal wave dispersion relation links the wave numbers $k$ and $m$, the Coriolis parameter $f$, and the buoyancy frequency N:

# Non-dimensional internal wave parameters
m = 16      # vertical wavenumber
k = 8       # horizontal wavenumber
f = coriolis.f

# Dispersion relation for inertia-gravity waves
ω² = (N^2 * k^2 + f^2 * m^2) / (k^2 + m^2)

ω = sqrt(ω²)

We define a Gaussian envelope for the wave packet so that we can observe wave propagation.

# Some Gaussian parameters
A = 1e-9
δ = grid.Lx / 15

# A Gaussian envelope centered at $(x, z) = (0, 0)$.
a(x, z) = A * exp( -( x^2 + z^2 ) / 2δ^2 )

An inertia-gravity wave is a linear solution to the Boussinesq equations. In order that our initial condition excites an inertia-gravity wave, we initialize the velocity and buoyancy perturbation fields to be consistent with the pressure field $p = a \, \cos(kx + mx - ωt)$ at $t=0$. These relations are sometimes called the "polarization relations". At $t=0$, the polarization relations yield

u₀(x, y, z) = a(x, z) * k * ω   / (ω^2 - f^2) * cos(k*x + m*z)
v₀(x, y, z) = a(x, z) * k * f   / (ω^2 - f^2) * sin(k*x + m*z)
w₀(x, y, z) = a(x, z) * m * ω   / (ω^2 - N^2) * cos(k*x + m*z)
b₀(x, y, z) = a(x, z) * m * N^2 / (ω^2 - N^2) * sin(k*x + m*z)

set!(model, u=u₀, v=v₀, w=w₀, b=b₀)

Recall that the buoyancy b is a perturbation, so that the total buoyancy field is $N^2 z + b$.

A wave packet on the loose

We're ready to release the packet. We build a simulation with a constant time-step,

simulation = Simulation(model, Δt = 0.02 * 2π/ω, stop_iteration = 100)

Simulation{IncompressibleModel{CPU, Float64}}
├── Model clock: time = 0 seconds, iteration = 0 
├── Next time step (Float64): 260.895 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: Inf years, stop iteration: 100
├── Diagnostics: OrderedCollections.OrderedDict with no entries
└── Output writers: OrderedCollections.OrderedDict with no entries

and add an output writer that saves the vertical velocity field every two iterations:

using Oceananigans.OutputWriters: JLD2OutputWriter, IterationInterval

simulation.output_writers[:velocities] = JLD2OutputWriter(model, model.velocities,
                                                          schedule = IterationInterval(2),
                                                            prefix = "internal_wave",
                                                             force = true)

JLD2OutputWriter scheduled on IterationInterval(2):
├── filepath: ./internal_wave.jld2
├── 3 outputs: (:u, :v, :w)
├── field slicer: FieldSlicer(:, :, :, with_halos=false)
├── array type: Array{Float32}
├── including: [:grid, :coriolis, :buoyancy, :closure]
└── max filesize: Inf YiB

With initial conditions set and an output writer at the ready, we run the simulation

run!(simulation)

[ Info: Simulation is stopping. Model iteration 100 has hit or exceeded simulation stop iteration 100.

Animating a propagating packet

To visualize the solution, we load snapshots of the data and use it to make contour plots of vertical velocity.

using JLD2, Plots, Printf, Oceananigans.Grids

WARNING: using Plots.grid in module ex-internal_wave conflicts with an existing identifier.

We use coordinate arrays appropriate for the vertical velocity field,

x, y, z = nodes(model.velocities.w)

open the jld2 file with the data,

file = jldopen(simulation.output_writers[:velocities].filepath)

# Extracts a vector of `iterations` at which data was saved.
iterations = parse.(Int, keys(file["timeseries/t"]))

51-element Array{Int64,1}:
   0
   2
   4
   6
   8
  10
  12
  14
  16
  18
   ⋮
  84
  86
  88
  90
  92
  94
  96
  98
 100

and makes an animation with Plots.jl:

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

    w = file["timeseries/w/$iter"][:, 1, :]
    t = file["timeseries/t/$iter"]

    contourf(x, z, w', title = @sprintf("ωt = %.2f", ω * t),
                      levels = range(-1e-8, stop=1e-8, length=10),
                       clims = (-1e-8, 1e-8),
                      xlabel = "x",
                      ylabel = "z",
                       xlims = (-π, π),
                       ylims = (-π, π),
                   linewidth = 0,
                       color = :balance,
                      legend = false,
                 aspectratio = :equal)
end

This page was generated using Literate.jl.