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.
Install dependencies
First let's make sure we have all required packages installed.
using Pkg
pkg"add Oceananigans, JLD2, Plots"Resolving package versions... No Changes to `/storage7/buildkite-agent/builds/tartarus-mit-edu-6/clima/oceananigans/docs/Project.toml` No Changes to `/storage7/buildkite-agent/builds/tartarus-mit-edu-6/clima/oceananigans/docs/Manifest.toml`
The physical domain
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)Internal wave parameters
Inertia-gravity waves propagate in fluids that are both (i) rotating, and (ii) density-stratified. We use Oceananigans' Coriolis abstraction to implement a background rotation rate:
coriolis = FPlane(f=0.2)FPlane{Float64}: f = 2.00e-01On 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 B(z) = N^2 * z, where z is the vertical coordinate and N is the "buoyancy frequency". This means that the modeled buoyancy field perturbs the basic state B(z).
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
N = 1 ## buoyancy frequency
B = BackgroundField(B_func, parameters=N)BackgroundField{typeof(Main.ex-internal_wave.B_func), Int64}
├── func: B_func
└── parameters: 1We are now ready to instantiate our model. We pass grid, coriolis, and B to the IncompressibleModel constructor. We add a small amount of IsotropicDiffusivity to keep the model stable during time-stepping, and specify that we're using a single tracer called b that we identify as buoyancy by setting buoyancy=BuoyancyTracer().
using Oceananigans.Advection
model = IncompressibleModel(
grid = grid,
advection = CenteredFourthOrder(),
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.1 * 2π/ω, stop_iteration = 15)Simulation{IncompressibleModel{CPU, Float64}}
├── Model clock: time = 0 seconds, iteration = 0
├── Next time step (Float64): 1.304 seconds
├── 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: 15
├── Diagnostics: OrderedCollections.OrderedDict with 1 entry:
│ └── nan_checker => NaNChecker
└── Output writers: OrderedCollections.OrderedDict with no entriesand 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(1),
prefix = "internal_wave",
force = true)JLD2OutputWriter scheduled on IterationInterval(1):
├── 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 YiBWith initial conditions set and an output writer at the ready, we run the simulation
run!(simulation)[ Info: Simulation is stopping. Model iteration 15 has hit or exceeded simulation stop iteration 15.
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, Printf, Plots, Oceananigans.GridsWARNING: 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"]))16-element Array{Int64,1}:
0
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15and makes an animation with Plots.jl:
anim = @animate for (i, iter) in enumerate(iterations)
@info "Drawing frame $i from iteration $iter..."
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.