In [694]:
from __future__ import division

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

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

In [695]:
data_path = 'outputs/beta-eady_C_kappa8_256.npz'

In [696]:
kappa = beta_eady['kappa']
qn = beta_eady['efunc'][1:-1]    # modal q
vtht = beta_eady['efunc'][0]     # theta +
vthb = beta_eady['efunc'][-1]    # theta -

In [697]:
# vertical coordinate
z = np.linspace(-1.,0.,N+3)
dz = z[2]-z[1]

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

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

In [710]:
def alpha_n(n,kappa):
""" Compute the inverse Helmholtz operator in
Fourier space """
return -1./( kappa**2 + (n*pi)**2 )

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

In [712]:
qabs = np.abs(q)
qphase = np.arctan2(q.imag,q.real)
psiabs = np.abs(psi)
psiphase = np.arctan2(psi.imag,psi.real)

vthtabs = np.abs(vtht)
vthbabs = np.abs(vthb)
vthtphase = np.arctan2(vtht.imag,vtht.real)
vthbphase = np.arctan2(vthb.imag,vthb.real)

In [713]:
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 [714]:
PV = Qabs*np.cos(kappa*X + Qphase)
PSI = Psiabs*np.cos(kappa*X + Psiphase)

In [715]:
plt.figure(figsize=(12,8))
plt.contour(X,Z,PSI,20,colors='k')
plt.contourf(X,Z,PV,np.linspace(-6.,6.,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.7, -.1, r"Eady Problem, $\kappa = 2$, Approximation C", size=25, rotation=0.,\
#            ha="center", va="center",\
#            bbox = dict(boxstyle="round",ec='k',fc='w'))
#plt.ylim(-.25,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_2_C_N64.eps',bbox_inches = 'tight')

In [ ]:



### Computing the buoyancy¶

1st method finite differences

In [705]:
def D_matrices(N):
''' create an N x N difference matrices

D1 :: first-order, centered difference
D2 :: second-order, centered difference

'''
D2 = np.zeros((N,N))
D1 = np.zeros((N,N))
for i in range(N):
D2[i,i] = -2.
if i<N-1:
D2[i,i+1],D1[i,i+1] = 1.,-1
if i>0:
D2[i,i-1],D1[i,i-1] = 1.,1.
return D1,D2

In [706]:
# vertical coordinate
z_fd = np.linspace(-1.,0.,N+3)

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

# grid
X,Z_fd= np.meshgrid(x,z_fd)

D1,_ = D_matrices(z_fd.size)
D1 = -D1/(2*dz)                # correct order for this grid

In [668]:
Bfd = np.empty(PSI.shape)
for i in range(x.size):
Bfd[:,i] = np.dot(D1,PSI[:,i])


Append boundary values

In [707]:
Bfd[0,:] = vthbabs*cos(kappa*x + vthbphase)
Bfd[-1,:] = vthtabs*cos(kappa*x + vthtphase)

In [ ]:


In [708]:
ct = np.linspace(-.7,.7,20)

In [709]:
plt.figure(figsize=(12,8))
plt.contour(X,Z_fd,Bfd,ct,colors='k')
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_b_psi_kappa_2_C_N128_fd.eps',bbox_inches = 'tight')


Method 2: use the Galerkin series

In [672]:
z2 = z_fd[1:-1]

for i in range(N+1):
if i == 0 :
b = np.zeros(z2.size) + 0.j
else:
b+= -i*pi*sqrt(2)*(qn[i] - sqrt(2)*vtht + sqrt(2)*((-1)**i * vthb))\
*sin(i*pi*z2)*alpha_n(i,kappa)

In [673]:
b = np.append(vthb,b)
b = np.append(b,vtht)
z2 = np.append(-1,z2)
z2 = np.append(z2,0)
babs = np.abs(b)
bphase = np.arctan2(b.imag,b.real)

# grid
X2,Z2 = np.meshgrid(x,z2)

In [674]:
babs = babs.repeat(x.size).reshape(z.size,x.size)
bphase = bphase.repeat(x.size).reshape(z.size,x.size)

In [688]:
B2 = babs*np.cos(kappa*X + bphase)

In [689]:
plt.figure(figsize=(12,8))
plt.contour(X2,Z2,B2,ct,colors='k')
#plt.contour(X,Z_fd,Bfd,ct,colors='m')
#plt.ylim(-.25,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_b_psi_kappa_2_C_N128.eps',bbox_inches = 'tight')


### Trying to remove gibbs oscillations: Laczos $\sigma$ factor¶

In [677]:
z3 = z_fd[1:-1]

for i in range(N+1):
if i == 0 :
b3 = np.zeros(z3.size) + 0.j
else:
b3+= -i*pi*sqrt(2)*sinc(i/((N+1)))*(qn[i] - sqrt(2)*vtht + sqrt(2)*((-1)**i * vthb))\
*sin(i*pi*z3)*alpha_n(i,kappa)

In [678]:
b3 = np.append(vthb,b3)
b3 = np.append(b3,vtht)
z3 = np.append(-1,z3)
z3 = np.append(z3,0)
babs = np.abs(b3)
bphase = np.arctan2(b3.imag,b3.real)

# grid
X3,Z3 = np.meshgrid(x,z3)

In [679]:
babs = babs.repeat(x.size).reshape(z3.size,x.size)
bphase= bphase.repeat(x.size).reshape(z3.size,x.size)

In [680]:
B3 = babs*np.cos(kappa*X3 + bphase)

In [681]:
plt.figure(figsize=(12,8))
plt.contour(X3,Z3,B3,ct,colors='k')
plt.contour(X,Z_fd,Bfd,ct,colors='m')
#plt.ylim(-.25,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_b_psi_kappa_2_C_N128.eps',bbox_inches = 'tight')

In [416]:
D1,D2 = D_matrices(z.size)
D1 = np.matrix(D1)/(2*dz)
D1[0,0],D1[0,1] = 1./dz,-1/dz
D1[-1,-1],D1[-1,-2] = -1./dz,1/dz

In [417]:
psiabs = np.abs(psi)
psiphase = np.arctan2(psi.imag,psi.real)
psiabs_prime_fd = np.array( D1*np.matrix(psiabs).T )
psiphase_prime_fd = np.array( D1*np.matrix(psiphase).T )

In [418]:
vthtabs = np.abs(vtht)
vthbabs = np.abs(vthb)
vthtphase = np.arctan2(vtht.imag,vtht.real)
vthbphase = np.arctan2(vthb.imag,vthb.real)

In [419]:
vthtabs = vthtabs.repeat(x.size).reshape(1,x.size)
vthbabs = vthbabs.repeat(x.size).reshape(1,x.size)

vthtphase = vthtphase.repeat(x.size).reshape(1,x.size)
vthbphase = vthbphase.repeat(x.size).reshape(1,x.size)

In [420]:
psiabs = psiabs.repeat(x.size).reshape(z.size,x.size)
psiabs_prime_fd = psiabs_prime_fd.repeat(x.size).reshape(z.size,x.size)
psiphase_prime_fd = psiphase_prime_fd.repeat(x.size).reshape(z.size,x.size)

In [421]:
Bfd = psiabs_prime_fd*cos(kappa*X + psiphase) - psiabs*psiphase_prime_fd*sin(kappa*X + psiphase)

In [422]:
Bfd = np.concatenate((vthbabs*cos(kappa*x+vthbphase),Bfd),axis=0)

In [424]:
# vertical coordinate
dz = .01   # vertical resolution
z = np.linspace(-1.+dz,-dz,100)

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

z = np.append(-1.,z)

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

In [ ]:


In [368]:
plt.figure(figsize=(12,8))
plt.contour(X,Z,Bfd,20,colors='k')
#plt.contourf(X,Z,PV,np.linspace(-2.,2.,9),cmap='RdBu_r',extend='both')
#cb = plt.colorbar(extend='both',shrink=.9)
#cb.ax.text(.0,1.075,'PV',rotation=0,fontsize=25)
plt.text(1.7, -.1, r"Eady Problem, $\kappa = 2$, Approximation C", size=25, rotation=0.,\
ha="center", va="center",\
bbox = dict(boxstyle="round",ec='k',fc='w'))
#plt.ylim(-.25,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_b_psi_kappa_2_C_N128_fd.eps',bbox_inches = 'tight')

In [ ]:


In [ ]: