Internal wave example
In this example, we initialize an internal wave packet in two-dimensions and watch is propagate.
using Oceananigans, PyPlot, Printf[ Info: Installing matplotlib via the Conda matplotlib package...
[ Info: Running `conda install -q -y matplotlib` in root environment
Collecting package metadata (current_repodata.json): ...working... done
Solving environment: ...working... done
## Package Plan ##
environment location: /home/travis/.julia/conda/3
added / updated specs:
- matplotlib
The following packages will be downloaded:
package | build
---------------------------|-----------------
cycler-0.10.0 | py37_0 13 KB
dbus-1.13.12 | h746ee38_0 611 KB
expat-2.2.6 | he6710b0_0 146 KB
fontconfig-2.13.0 | h9420a91_0 227 KB
freetype-2.9.1 | h8a8886c_1 550 KB
glib-2.63.1 | h5a9c865_0 3.4 MB
gst-plugins-base-1.14.0 | hbbd80ab_1 4.8 MB
gstreamer-1.14.0 | hb453b48_1 3.1 MB
icu-58.2 | h9c2bf20_1 10.3 MB
kiwisolver-1.1.0 | py37he6710b0_0 82 KB
libpng-1.6.37 | hbc83047_0 278 KB
libuuid-1.0.3 | h1bed415_2 15 KB
libxcb-1.13 | h1bed415_1 421 KB
libxml2-2.9.9 | hea5a465_1 1.6 MB
matplotlib-3.1.1 | py37h5429711_0 5.0 MB
pcre-8.43 | he6710b0_0 209 KB
pyparsing-2.4.5 | py_0 62 KB
pyqt-5.9.2 | py37h05f1152_2 4.5 MB
python-dateutil-2.8.1 | py_0 224 KB
pytz-2019.3 | py_0 231 KB
qt-5.9.7 | h5867ecd_1 68.5 MB
sip-4.19.8 | py37hf484d3e_0 274 KB
tornado-6.0.3 | py37h7b6447c_0 584 KB
------------------------------------------------------------
Total: 104.9 MB
The following NEW packages will be INSTALLED:
cycler pkgs/main/linux-64::cycler-0.10.0-py37_0
dbus pkgs/main/linux-64::dbus-1.13.12-h746ee38_0
expat pkgs/main/linux-64::expat-2.2.6-he6710b0_0
fontconfig pkgs/main/linux-64::fontconfig-2.13.0-h9420a91_0
freetype pkgs/main/linux-64::freetype-2.9.1-h8a8886c_1
glib pkgs/main/linux-64::glib-2.63.1-h5a9c865_0
gst-plugins-base pkgs/main/linux-64::gst-plugins-base-1.14.0-hbbd80ab_1
gstreamer pkgs/main/linux-64::gstreamer-1.14.0-hb453b48_1
icu pkgs/main/linux-64::icu-58.2-h9c2bf20_1
kiwisolver pkgs/main/linux-64::kiwisolver-1.1.0-py37he6710b0_0
libpng pkgs/main/linux-64::libpng-1.6.37-hbc83047_0
libuuid pkgs/main/linux-64::libuuid-1.0.3-h1bed415_2
libxcb pkgs/main/linux-64::libxcb-1.13-h1bed415_1
libxml2 pkgs/main/linux-64::libxml2-2.9.9-hea5a465_1
matplotlib pkgs/main/linux-64::matplotlib-3.1.1-py37h5429711_0
pcre pkgs/main/linux-64::pcre-8.43-he6710b0_0
pyparsing pkgs/main/noarch::pyparsing-2.4.5-py_0
pyqt pkgs/main/linux-64::pyqt-5.9.2-py37h05f1152_2
python-dateutil pkgs/main/noarch::python-dateutil-2.8.1-py_0
pytz pkgs/main/noarch::pytz-2019.3-py_0
qt pkgs/main/linux-64::qt-5.9.7-h5867ecd_1
sip pkgs/main/linux-64::sip-4.19.8-py37hf484d3e_0
tornado pkgs/main/linux-64::tornado-6.0.3-py37h7b6447c_0
Preparing transaction: ...working... done
Verifying transaction: ...working... done
Executing transaction: ...working... doneNumerical, domain, and internal wave parameters
First, we pick some numerical and physical parameters for our model and its rotation rate.
Nx = 128 # resolution
Lx = 2π # domain extent6.283185307179586We set up an internal wave with 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.
# Non-dimensional internal wave parameters
m = 16 # vertical wavenumber
k = 1 # horizontal wavenumber
N = 1 # buoyancy frequency
f = 0.2 # inertial frequency0.2A Gaussian wavepacket
Next, we set up an initial condition corresponding to a propagating wave packet with a Gaussian envelope. The internal wave dispersion relation yields
ω² = (N^2 * k^2 + f^2 * m^2) / (k^2 + m^2)
# and thus
ω = sqrt(ω²)0.20913012351239907The internal wave polarization relations follow from the linearized Boussinesq equations,
U = k * ω / (ω^2 - f^2)
V = k * f / (ω^2 - f^2)
W = m * ω / (ω^2 - N^2)
B = m * N^2 / (ω^2 - N^2)-16.731770833333336Finally, we set-up a small-amplitude, Gaussian envelope for the wave packet
# Some Gaussian parameters
A, x₀, z₀, δ = 1e-9, Lx/2, -Lx/2, Lx/15
# A Gaussian envelope
a(x, z) = A * exp( -( (x - x₀)^2 + (z - z₀)^2 ) / 2δ^2 )a (generic function with 1 method)Create initial condition functions
u₀(x, y, z) = a(x, z) * U * cos(k*x + m*z)
v₀(x, y, z) = a(x, z) * V * sin(k*x + m*z)
w₀(x, y, z) = a(x, z) * W * cos(k*x + m*z)
b₀(x, y, z) = a(x, z) * B * sin(k*x + m*z) + N^2 * zb₀ (generic function with 1 method)We are now ready to instantiate our model on a uniform grid. We give the model a constant rotation rate with background vorticity f, use temperature as a buoyancy tracer, and use a small constant viscosity and diffusivity to stabilize the model.
model = Model(
grid = RegularCartesianGrid(size=(Nx, 1, Nx), length=(Lx, Lx, Lx)),
closure = ConstantIsotropicDiffusivity(ν=1e-6, κ=1e-6),
coriolis = FPlane(f=f),
tracers = :b,
buoyancy = BuoyancyTracer()
)Oceananigans.Model on a CPU architecture (time = 0.000 ns, iteration = 0)
├── grid: RegularCartesianGrid{Float64,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}
├── tracers: (:b,)
├── closure: ConstantIsotropicDiffusivity{Float64,NamedTuple{(:b,),Tuple{Float64}}}
├── buoyancy: BuoyancyTracer
├── coriolis: FPlane{Float64}
├── output writers: OrderedCollections.OrderedDict{Symbol,Oceananigans.AbstractOutputWriter} with no entries
└── diagnostics: OrderedCollections.OrderedDict{Symbol,Oceananigans.AbstractDiagnostic} with no entriesWe initialize the velocity and buoyancy fields with our internal wave initial condition.
set!(model, u=u₀, v=v₀, w=w₀, b=b₀)Some plotting utilities
To watch the wave packet propagate interactively as the model runs, we build some plotting utilities.
xplot(u) = repeat(dropdims(xnodes(u), dims=2), 1, u.grid.Nz)
zplot(u) = repeat(dropdims(znodes(u), dims=2), u.grid.Nx, 1)
function plot_field!(ax, w, t)
pcolormesh(xplot(w), zplot(w), interior(model.velocities.w)[:, 1, :])
xlabel(L"x")
ylabel(L"z")
title(@sprintf("\$ \\omega t / 2 \\pi = %.2f\$", t*ω/2π))
ax.set_aspect(1)
pause(0.1)
return nothing
end
close("all")
fig, ax = subplots();(PyPlot.Figure(PyObject <Figure size 640x480 with 1 Axes>), PyObject <matplotlib.axes._subplots.AxesSubplot object at 0x7f2c3b3c1590>)A wave packet on the loose
Finally, we release the packet and plot its trajectory:
for i = 1:10
time_step!(model, Nt = 200, Δt = 0.001 * 2π/ω)
plot_field!(ax, model.velocities.w, model.clock.time)
end
gcf()
This page was generated using Literate.jl.