pyqg: Python Quasigeostrophic Model

1
2
3
4
5
6

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.


1
2
3
4
5
6

Leave a Reply

Your email address will not be published. Required fields are marked *

This site uses Akismet to reduce spam. Learn how your comment data is processed.