Here will will use pyqg to reproduce the results of the paper:
J. C. Mcwilliams (1984). The emergence of isolated coherent vortices in turbulent flow. Journal of Fluid Mechanics, 146, pp 21-43 doi:10.1017/S0022112084001750
In [23]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import pyqg
McWilliams performed freely-evolving 2D turbulence ($R_d = \infty$, $\beta =0$) experiments on a $2\pi\times 2\pi$ periodic box.
In [24]:
# create the model object
m = pyqg.BTModel(L=2.*np.pi, nx=256,
beta=0., H=1., rek=0., rd=None,
tmax=40, dt=0.001, taveint=1,
ntd=4)
# in this example we used ntd=4, four threads
# if your machine has more (or fewer) cores available, you could try changing it
In [3]:
# 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)
In [4]:
# initialize the model with that initial condition
m.set_q(qi)
In [5]:
# define a quick function for plotting and visualize the initial condition
def plot_q(m, qmax=40):
fig, ax = plt.subplots()
pc = ax.pcolormesh(m.x,m.y,m.q.squeeze(), cmap='RdBu_r')
pc.set_clim([-qmax, qmax])
ax.set_xlim([0, 2*np.pi])
ax.set_ylim([0, 2*np.pi]);
ax.set_aspect(1)
plt.colorbar(pc)
plt.title('Time = %g' % m.t)
plt.show()
plot_q(m)
In [6]:
for _ in m.run_with_snapshots(tsnapstart=0, tsnapint=10):
plot_q(m)
The genius of McWilliams (1984) was that he has showed that the initial random vorticity field organizes itself into strong coherent vortices. This is true in significant part of the parameter space. This was previously suspected but unproven, mainly because people did not have computer resources to run the simulation long enough. Thirty years later we can perform such simulations in a couple of minutes on a laptop!
Also, note that the energy is nearly conserved, as it should be, and this is a nice test of the model.
In [7]:
energy = m.get_diagnostic('KEspec')
enstrophy = m.get_diagnostic('Ensspec')
In [11]:
# this makes it easy to calculate an isotropic spectrum
from pyqg import diagnostic_tools as tools
kr, energy_iso = tools.calc_ispec(m,energy.squeeze())
_, enstrophy_iso = tools.calc_ispec(m,enstrophy.squeeze())
In [17]:
ks = np.array([3.,80])
es = 5*ks**-4
plt.loglog(kr,energy_iso)
plt.loglog(ks,es,'k--')
plt.text(2.5,.0001,r'$k^{-4}$',fontsize=20)
plt.ylim(1.e-10,1.e0)
plt.xlabel('wavenumber')
plt.title('Energy Spectrum')
Out[17]:
In [20]:
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')
Out[20]:
In [ ]: