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.