Edit on GitHub

Internal wave example

In this example, we initialize an internal wave packet in two-dimensions and watch it propagate.

using Oceananigans, Oceananigans.Grids, Plots, Printf
┌ Debug: Precompiling Plots [91a5bcdd-55d7-5caf-9e0b-520d859cae80]
└ @ Base loading.jl:1273
┌ Debug: Precompiling GeometryTypes [4d00f742-c7ba-57c2-abde-4428a4b178cb]
└ @ Base loading.jl:1273
┌ Debug: Precompiling ColorTypes [3da002f7-5984-5a60-b8a6-cbb66c0b333f]
└ @ Base loading.jl:1273
┌ Debug: Precompiling FixedPointNumbers [53c48c17-4a7d-5ca2-90c5-79b7896eea93]
└ @ Base loading.jl:1273
┌ Debug: Precompiling FFMPEG [c87230d0-a227-11e9-1b43-d7ebe4e7570a]
└ @ Base loading.jl:1273
┌ Debug: Precompiling FFMPEG_jll [b22a6f82-2f65-5046-a5b2-351ab43fb4e5]
└ @ Base loading.jl:1273
┌ Debug: Precompiling libass_jll [0ac62f75-1d6f-5e53-bd7c-93b484bb37c0]
└ @ Base loading.jl:1273
┌ Debug: Precompiling FreeType2_jll [d7e528f0-a631-5988-bf34-fe36492bcfd7]
└ @ Base loading.jl:1273
┌ Debug: Precompiling Bzip2_jll [6e34b625-4abd-537c-b88f-471c36dfa7a0]
└ @ Base loading.jl:1273
┌ Debug: Precompiling FriBidi_jll [559328eb-81f9-559d-9380-de523a88c83c]
└ @ Base loading.jl:1273
┌ Debug: Precompiling libfdk_aac_jll [f638f0a6-7fb0-5443-88ba-1cc74229b280]
└ @ Base loading.jl:1273
┌ Debug: Precompiling LAME_jll [c1c5ebd0-6772-5130-a774-d5fcae4a789d]
└ @ Base loading.jl:1273
┌ Debug: Precompiling libvorbis_jll [f27f6e37-5d2b-51aa-960f-b287f2bc3b7a]
└ @ Base loading.jl:1273
┌ Debug: Precompiling Ogg_jll [e7412a2a-1a6e-54c0-be00-318e2571c051]
└ @ Base loading.jl:1273
┌ Debug: Precompiling LibVPX_jll [dd192d2f-8180-539f-9fb4-cc70b1dcf69a]
└ @ Base loading.jl:1273
┌ Debug: Precompiling x264_jll [1270edf5-f2f9-52d2-97e9-ab00b5d0237a]
└ @ Base loading.jl:1273
┌ Debug: Precompiling x265_jll [dfaa095f-4041-5dcd-9319-2fabd8486b76]
└ @ Base loading.jl:1273
┌ Debug: Precompiling OpenSSL_jll [458c3c95-2e84-50aa-8efc-19380b2a3a95]
└ @ Base loading.jl:1273
┌ Debug: Precompiling Opus_jll [91d4177d-7536-5919-b921-800302f37372]
└ @ Base loading.jl:1273
┌ Debug: Precompiling RecipesBase [3cdcf5f2-1ef4-517c-9805-6587b60abb01]
└ @ Base loading.jl:1273
┌ Debug: Precompiling PlotUtils [995b91a9-d308-5afd-9ec6-746e21dbc043]
└ @ Base loading.jl:1273
┌ Debug: Precompiling Colors [5ae59095-9a9b-59fe-a467-6f913c188581]
└ @ Base loading.jl:1273
.
┌ Debug: Precompiling PlotThemes [ccf2f8ad-2431-5c83-bf29-c5338b663b6a]
└ @ Base loading.jl:1273
┌ Debug: Precompiling Showoff [992d4aef-0814-514b-bc4d-f2e9a6c4116f]
└ @ Base loading.jl:1273
┌ Debug: Precompiling StatsBase [2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91]
└ @ Base loading.jl:1273
┌ Debug: Precompiling DataAPI [9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a]
└ @ Base loading.jl:1273
┌ Debug: Precompiling SortingAlgorithms [a2af1166-a08f-5f64-846c-94a0d3cef48c]
└ @ Base loading.jl:1273
┌ Debug: Precompiling Missings [e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28]
└ @ Base loading.jl:1273
┌ Debug: Precompiling NaNMath [77ba4419-2d1f-58cd-9bb1-8ffee604a2e3]
└ @ Base loading.jl:1273
┌ Debug: Precompiling Measures [442fdcdd-2543-5da2-b0f3-8c86c306513e]
└ @ Base loading.jl:1273
┌ Debug: Precompiling RecipesPipeline [01d81517-befc-4cb6-b9ec-a95719d0359c]
└ @ Base loading.jl:1273
┌ Debug: Precompiling GR [28b8d3ca-fb5f-59d9-8090-bfdbd6d07a71]
└ @ Base loading.jl:1273

Numerical, 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 extent

We 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 frequency

A 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(ω²)

The 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)

Finally, 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 )

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 * z

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 = IncompressibleModel(    grid = RegularCartesianGrid(size=(Nx, 1, Nx), extent=(Lx, Lx, Lx)),
                             closure = ConstantIsotropicDiffusivity(ν=1e-6, κ=1e-6),
                            coriolis = FPlane(f=f),
                             tracers = :b,
                            buoyancy = BuoyancyTracer())
┌ Debug: Planning transforms for PressureSolver{HorizontallyPeriodic, CPU}...
└ @ Oceananigans.Solvers ~/build/climate-machine/Oceananigans.jl/src/Solvers/horizontally_periodic_pressure_solver.jl:14
┌ Debug: Planning transforms for PressureSolver{HorizontallyPeriodic, CPU} done!
└ @ Oceananigans.Solvers ~/build/climate-machine/Oceananigans.jl/src/Solvers/horizontally_periodic_pressure_solver.jl:20
.

We initialize the velocity and buoyancy fields with our internal wave initial condition.

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

A wave packet on the loose

Finally, we release the packet and watch it go!

simulation = Simulation(model, Δt = 0.001 * 2π/ω, stop_iteration = 0,
                        progress_frequency = 20)

anim = @animate for i=1:100
    simulation.stop_iteration += 20
    run!(simulation)

    x, z = xnodes(Cell, model.grid)[:], znodes(Face, model.grid)[:]
    w = model.velocities.w
    heatmap(x, z, w.data[1:Nx, 1, 1:Nx+1]', title=@sprintf("t = %.2f", model.clock.time),
            xlabel="x", ylabel="z", c=:balance, clims=(-1e-8, 1e-8))
end

This page was generated using Literate.jl.