# On Galerkin approximations for the QG equations

## Supplementary material for subsection on the Eady model

</h3>Cesar B. Rocha*</h3> </h3>, William R. Young, and Ian Grooms </h3>

</h4>Winter 2015 </h4>

*Scripps Institution of Oceanography, University of California, San Diego, 9500 Gilman Dr. MC 0213, La Jolla, CA/USA, [email protected]

In [104]:
from __future__ import division

import numpy as np
from numpy import pi, sqrt,cos

import matplotlib.pyplot as plt
plt.rcParams.update({'font.size': 25, 'legend.handlelength'  : 1.25})
%matplotlib inline

import seaborn as sns
#sns.set(style="darkgrid")
sns.set_context("paper", font_scale=5, rc={"lines.linewidth": 1.5})


We computed and saved the eigenvalue and eigenfunction using Approximation A with $\mathrm{N}=64$ using the linear_eady notebook. We now load the data and plot the wavestructure in the $(x,z)$ plane.

In [105]:
data_path = 'outputs/eady_A_k8_64_efunc.npz'

In [106]:
kappa = eady['kappa']


Seting up the domain grid in $(x,z)$ space

In [107]:
# vertical coordinate
dz = 1./N   # vertical resolution
z = np.arange(-dz/2,-1.-dz/2.,-dz)  # level array

# horizontal coordinate
x = np.linspace(0,np.pi,100)

# grid
X,Z = np.meshgrid(x,z)


The eigenfunctions are the modal coefficients $\,\breve{q}_n$. We use the standard modes $\mathsf{p}_0 = 1$ and $\mathsf{p}_n = \sqrt{2}\,\cos(n \pi z)$, $n>1$ to compute the eigenstructure in physical space:

$$\hat{q}(z) = \sum_{n=0}^{\mathrm{N}} \breve{q}_n\, \mathsf{p}_n \,\,\,.$$

In approximation A, the modal streamfunction $\breve{\psi}_n$ is related to the PV by

$$\breve{\psi}_n = \alpha_n\, \breve{q}_n = - (\kappa^2 + (n \pi)^2)^{-1}\,\,\breve{q}_n\,\,.$$

The streamfunction $\hat{\psi}$ is then

$$\hat{q}(z) = \sum_{n=0}^{\mathrm{N}} \breve{q}_n\, \mathsf{p}_n \,.$$$$\hat{\psi}(z) = \sum_{n=0}^{\mathrm{N}} \breve{\psi}_n\, \mathsf{p}_n \,.$$
In [108]:
def alpha_n(n,kappa):
""" Compute the inverse Helmholtz operator in
Fourier space """
return -1./( kappa**2 + (n*pi)**2 )

In [109]:
for i in range(N):
if i == 0 :
q = qn[0]*np.ones(z.size)
psi = qn[0]*np.ones(z.size)*alpha_n(0,kappa)
else:
q+= sqrt(2)*qn[i]*cos(i*pi*z)
psi+= sqrt(2)*qn[i]*cos(i*pi*z)*alpha_n(i,kappa)


Note that $q$ and $\psi$ are complex. The wavestructure is then

$$\psi(x,z) = |\psi(z)| \exp({ik\,x + P_{\psi}(z)})\,, \qquad \text{and}\qquad q(x,z) = |q(z)|\, \exp({ik\,x + P_{q}(z)})$$

where vertical bars denote the magnitude and the phases are

$$P_{\psi}(z) = \text{tan^{-1}}\frac{\text{Im}(\hat{\psi})}{\text{Re}(\hat{\psi})}\,, \qquad \text{and}\qquad P_{q}(z) = \text{tan^{-1}}\frac{\text{Im}(\hat{q})}{\text{Re}(\hat{q})}$$
In [100]:
qabs = np.abs(q)
qphase = np.arctan2(q.imag,q.real)
psiabs = np.abs(psi)
psiphase = np.arctan2(psi.imag,psi.real)

In [101]:
qabs = qabs.repeat(x.size).reshape(z.size,x.size)
qphase = qphase.repeat(x.size).reshape(z.size,x.size)
psiabs = psiabs.repeat(x.size).reshape(z.size,x.size)
psiphase = psiphase.repeat(x.size).reshape(z.size,x.size)

In [102]:
PV = qabs*np.cos(kappa*X + qphase)
PSI = psiabs*np.cos(kappa*X + psiphase)

In [103]:
plt.figure(figsize=(11,8))
plt.contour(X,Z,PSI,colors='k')
plt.contourf(X,Z,PV,np.linspace(-2.,2.,9),cmap='RdBu_r',extend='both')
#plt.text(-0.325,zc,r' $z_c \rightarrow$',fontsize=25)
cb = plt.colorbar(extend='both',shrink=.9)
cb.ax.text(.0,1.075,'PV',rotation=0,fontsize=25)
plt.text(1.725, -.9, r"Eady Problem, $\kappa = 8$, Approximation A", size=25, rotation=0.,\
ha="center", va="center",\
bbox = dict(boxstyle="round",ec='k',fc='w'))
plt.ylim(-1.,0)
plt.xticks([0.,pi/4,pi/2,3*pi/4,pi],[r'$0$',r'$\pi/4$',r'$\pi/2$',\
r'$3\,\pi/4$',r'$\pi$'])
plt.ylabel('$z/H$')
plt.xlabel(r'$x/L_d$')
#plt.savefig('figs/wave-structure_pv_psi_kappa_8_A_N64.eps')

Out[103]:
<matplotlib.text.Text at 0x112324b10>
In [ ]: