That's Geophysical Fluid Dynamics.

Barotropic Model

Turbulent barotropic flow can result in the emergence of coherent vortices. We will try to recreate this experiment from this James McWilliams paper.

McWilliams performed freely-evolving 2D turbulence experiments on a \(2\pi \times 2\pi\) periodic box:

# create the model object
kwargs = {'L': 2.*np.pi, 'nx': 256, 'beta': 0., 'H': 1., 'rek': 0., 'rd': None, 'tmax': 40., 'dt': 0.001, 'taveint': 1., 'ntd': 4}
m = pyqg.BTModel(**kwargs)
# in this example we used ntd=4, four threads
# if your machine has more (or fewer) cores available, you could try changing it

The initial condition is random, but defined by the following spectrum:

$$ |\hat{\psi}|^2 = A\kappa^{-1}\left[ 1+\left( \frac{\kappa}{6} \right)^4 \right]^{-1} $$

where \(\kappa\) is the wavenumber magnitude. \(A\) is chosen such that the initial kinetic energy is 0.5. This code generates the initial condition:

# generate McWilliams 84 IC condition

fk = (m.wv != 0)
ckappa = np.zeros_like(m.wv2)
ckappa[fk] = np.sqrt( m.wv2[fk]*(1. + (m.wv2[fk]/36.)**2) )**-1

nhx,nhy = m.wv2.shape

Pi_hat = np.random.randn(nhx,nhy)*ckappa +1j*np.random.randn(nhx,nhy)*ckappa

Pi = m.ifft( Pi_hat[np.newaxis,:,:] )
Pi = Pi - Pi.mean()
Pi_hat = m.fft( Pi )
KEaux = m.spec_var( m.wv*Pi_hat )

pih = ( Pi_hat/np.sqrt(KEaux) )
qih = -m.wv2*pih
qi = m.ifft(qih)

# initialize the model with that initial condition
m.set_q(qi)

# define a quick function for plotting and visualize the initial condition
def plot_q(m, qmax=40, show_plot=True):
    extent = [0, 2*np.pi, 2*np.pi, 0]
    fig, ax = plt.subplots()
    pc = ax.imshow(m.q.squeeze(), cmap='RdBu_r', extent=extent)
    pc.set_clim([-qmax, qmax])
    ax.set_aspect(1)
    plt.colorbar(pc)
    plt.title('Time = %g' % m.t)
    
    if show_plot:
        plt.show()
    return fig, ax, pc

plot_q(m)

pyqg comes with a run_with_snapshots feature to enable actions to be performed at select time steps. We can create a movie:

import matplotlib.animation as manimation

FFMpegWriter = manimation.writers['ffmpeg']
writer = FFMpegWriter(fps=5)
vid_title = 'barotropic_movie'

fig, ax, l = plot_q(m, show_plot=False)
with writer.saving(fig, "{}.mp4".format(vid_title), 100):
    writer.grab_frame() # initial data
    for _ in m.run_with_snapshots(tsnapstart=0, tsnapint=1):
        l.set_data(m.q.squeeze())
        plt.title('Time = %g' % m.t)
        writer.grab_frame()

And we get a movie that looks like this:

We can see that small scale variability disappears as large scale vortices form. Let’s calculate the energy and enstrophy spectra. First we set up:

energy = m.get_diagnostic('KEspec')
enstrophy = m.get_diagnostic('Ensspec')
from pyqg import diagnostic_tools as tools
kr, energy_iso = tools.calc_ispec(m,energy.squeeze())
_, enstrophy_iso = tools.calc_ispec(m,enstrophy.squeeze())

This gives the energy spectrum:

slope_1 = -3.6
slope_2 = -2.6
ks1 = np.array([10., 80.])
es1 = 5*ks1**slope_1
ks2 = np.array([4., 11.])
es2 = 4 * ks2**slope_2

plt.loglog(kr, energy_iso, 'k-')

plt.text(30, 1.0e-6, r'$k^{-3.6}$',fontsize=12, color='red')
plt.loglog(ks1, es1, 'r--')

plt.text(12, 1.0e-2, r'$k^{-2.6}$', fontsize=12, color='blue')
plt.loglog(ks2, es2, 'b--')

plt.ylim(1.e-10, 1.e0)
plt.xlabel('wavenumber')
plt.title('Energy Spectrum')
plt.show()

And this gives the enstrophy spectrum:

ks = np.array([3.,80])
es = 5*ks**(-5./3)
plt.loglog(kr,enstrophy_iso)
plt.loglog(ks,es,'k--')
plt.text(5.5,.01,r'$k^{-5/3}$',fontsize=20)
plt.ylim(1.e-3,1.e0)
plt.xlabel('wavenumber')
plt.title('Enstrophy Spectrum')
plt.show()

Next we’ll move into Surface Quasigeostrophy, or SQG, to see how fully 3D flows can be induced via buoyancy advection at the boundary.