-
Notifications
You must be signed in to change notification settings - Fork 2
/
qg_1L_channel_stab.py
117 lines (90 loc) · 2.78 KB
/
qg_1L_channel_stab.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
import sys, slepc4py
slepc4py.init(sys.argv)
from petsc4py import PETSc
from slepc4py import SLEPc
import numpy as np
import qg_1L_fxs as fx
import matplotlib.pyplot as plt
Print = PETSc.Sys.Print # For printing with only 1 processor
rank = PETSc.COMM_WORLD.Get_rank()
opts = PETSc.Options()
nEV = opts.getInt('nev', 5)
Ly = 350e03
Lj = 20e03
Hm = 500
Ny = opts.getInt('Ny',400)
y = fx.vec(0, float(Ly)/float(Ny), Ny+1)
hy = float(Ly)/float(Ny)
Dy = fx.Dy(hy, Ny)
Dy2 = fx.Dy(hy, Ny, Dy2=True)
f0 = 1.e-4
bet = 0
g0 = 9.81
Phi = fx.Phi(y, Ly, Lj, Ny)
U = fx.U(Dy, Phi, Ny)
etaB = fx.etaB(y)
F0 = f0**2/(g0*Hm)
dkk = 2e-2
kk = np.arange(dkk,2+dkk,dkk)/Lj
nk = len(kk)
# Temporary vector to store Dy2*Phi
temp = PETSc.Vec().createMPI(Ny-1,comm=PETSc.COMM_WORLD)
Dy2.mult(Phi,temp)
temp.assemble()
temp2 = PETSc.Vec().createMPI(Ny+1, comm=PETSc.COMM_WORLD)
ts,te = temp2.getOwnershipRange()
if ts == 0: temp2[0] = temp[0]; ts+=1
if te == Ny+1: temp2[Ny] = temp[Ny-2]; te -= 1
for i in xrange(ts,te):
temp2[i] = temp[i-1]
temp2.assemble()
Q = temp2 - F0*Phi + bet*y + f0/Hm*etaB # size Ny+1
grow = np.zeros([nEV,nk])
freq = np.zeros([nEV,nk])
mode = np.zeros([Ny+1,nEV,nk], dtype=complex)
# grOut = open('grow.dat', 'wb')
# frOut = open('freq.dat', 'wb')
# mdOut = open('mode.dat', 'wb')
cnt = 0
for kx in kk[0:nk]:
kx2=kx**2
Lap = fx.Lap(Dy2, kx2, Ny)
B = fx.B(Lap, F0, Ny)
A = fx.A(U,Lap, F0,Dy,Q,Ny)
# Set up slepc, generalized eig problem
E = SLEPc.EPS(); E.create(comm=SLEPc.COMM_WORLD)
E.setOperators(A,B); E.setDimensions(nEV, PETSc.DECIDE)
E.setProblemType(SLEPc.EPS.ProblemType.GNHEP); E.setFromOptions()
E.setWhichEigenpairs(SLEPc.EPS.Which.LARGEST_IMAGINARY)
E.setTolerances(1e-8, max_it=50)
E.solve()
nconv = E.getConverged()
vr, wr = A.getVecs()
vi, wi = A.getVecs()
if nconv <= nEV: evals = nconv
else: evals = nEV
for i in xrange(evals):
eigVal = E.getEigenvalue(i)
grow[i,cnt] = eigVal.imag*kx
freq[i,cnt] = eigVal.real*kx
if rank == 0:
# grow[i,cnt].tofile(grOut) # save values to file
# freq[i,cnt].tofile(frOut)
plt.plot(kk*Lj, grow[i]*3600*24, 'o')
eigVec=E.getEigenvector(i,vr,vi)
start,end = vi.getOwnershipRange()
if start == 0: mode[0,i,cnt] = 0; start+=1
if end == Ny: mode[Ny,i,cnt] = 0; end -=1
for j in xrange(start,end):
mode[j,i,cnt] = 1j*vi[j]; mode[j,i,cnt] = vr[j]
# if rank == 0:
# mode[j,i,cnt].tofile(mdOut)
cnt = cnt+1
# grOut.close(); frOut.close(); mdOut.close()
if rank == 0:
ky = np.pi/Ly
plt.ylabel('1/day')
plt.xlabel('k')
plt.title('Growth Rate: 1-Layer QG')
plt.savefig('Grow1L_QG.eps', format='eps', dpi=1000)
#plt.show()