-
Notifications
You must be signed in to change notification settings - Fork 0
/
lorenz96_compute_secondary_instabilities.py
executable file
·118 lines (101 loc) · 4.69 KB
/
lorenz96_compute_secondary_instabilities.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
118
#!/home/user/anaconda3/bin/python
import numpy as np
import os
import l96
from scipy import triu
import scipy.linalg as linalg
from itertools import product
paraL96_2lay = {'F1' : 64,
'F2' : 0,
'b' : 10,
'c' : 10,
'h' : 1,
'dimX': 20,
'dimY' : 10,
'RescaledY' : False,
'expname' : 'secondaryinstabilities_2layer',
'time' : np.arange(0,50,0.01),
'spinup' : 100,
'2lay' : True,
'integrator': 'classic'
}
paraL96_1lay = {'F1' : 8,
'F2' : 0,
'b' : 10,
'c' : 10,
'h' : 1,
'dimX': 44,
'dimY' : 10,
'RescaledY' : False,
'expname' : 'secondaryinstabilities_1layer',
'time' : np.arange(0,200,0.1),
'spinup' : 100,
'2lay' : False,
'integrator': 'classic'
}
testzeroclv=True
hs=[ 0.0 , 0.0625, 0.125 , 0.25 , 0.5 , 1. ]
experiments = [paraL96_2lay,paraL96_1lay]
steplengthforsecondorder = np.arange(0,100,1)
precision='float64'
for paraL96,h in product(experiments ,hs):
if not paraL96['2lay'] and not h == 1.0: print("1 lay only with h = 1.");break
savename=paraL96['expname']+"_h_"+str(h)
spinup = paraL96['spinup']
paraL96=np.load(savename+"/paraL96.npy")
paraL96=paraL96[()]
# M number exponents
if paraL96['2lay']:
contract_func = lambda x,y: l96.contract_hess_l96_2layer_v2(x,y,paraL96)
M = paraL96['dimX'] + paraL96['dimX']*paraL96['dimY'] # -1 full spectrum
dimN = paraL96['dimX'] + paraL96['dimX']*paraL96['dimY'] # -1 full spectrum
else:
contract_func = l96.contract_hess_l96_1layer_v2
M = paraL96['dimX']
dimN = paraL96['dimX']
np.save(savename+'/steplengthforsecondorder',steplengthforsecondorder)
t = paraL96['time']
CLV = np.memmap(savename+'/CLV.dat',mode='r',shape=(len(t),dimN,M),dtype='float64', order = 'C')
lyaploc_clv = np.memmap(savename+'/lyaploc_clv',mode='r',shape=(len(t),M),dtype='float64', order = 'C')
invCLV = np.memmap(savename+'/invCLV.dat',mode='r',shape=(len(t),dimN,M),dtype='float64', order = 'C')
lensteparray = len(steplengthforsecondorder)
contracted_CLVs = np.memmap(savename+'/contracted_clvs.dat',mode='w+',shape=(len(t),dimN,M),dtype=precision, order = 'C') #(final_clv,init_clv)
solution = np.memmap(savename+'/solution.dat',mode='w+',shape=(len(t),lensteparray,dimN,M),dtype=precision, order = 'C')
full_solution = np.memmap(savename+'/full_solution.dat',mode='w+',shape=(len(t),lensteparray,dimN,M),dtype=precision, order = 'C')
#normalized_solution = np.memmap(savename+'/normalized_solution.dat',mode='w+',shape=(len(t),lensteparray,dimN,M),dtype=precision, order = 'C')
growth = np.memmap(savename+'/growth.dat',mode='w+',shape=(len(t),lensteparray,M),dtype=precision, order = 'C')
for tn, (ts,te) in enumerate(zip(t[0:-2],t[1:-1])):
print(tn)
contracted_CLVs[tn,:,:]=1/2*contract_func(invCLV[tn,:,:],CLV[tn,:,:])
if tn % 1000: np.memmap.flush(contracted_CLVs)
np.memmap.flush(contracted_CLVs)
for tn, (ts,te) in enumerate(zip(t[0:-steplengthforsecondorder.max()-1],t[1:-steplengthforsecondorder.max()])):
dtau=te -ts
print(tn,paraL96['h'],paraL96['2lay'])
for n_ind, n_step in enumerate(steplengthforsecondorder):
if n_ind == 0: solution[tn,n_ind,:,:] = 0
else:
tni = tn
solution[tn,n_ind,:,:] = (solution[tn,n_ind-1,:,:] +
dtau * np.multiply(np.exp(-dtau*np.memmap.sum(lyaploc_clv[tni:tn+n_step,:,np.newaxis], axis =0)),np.multiply(contracted_CLVs[tn+n_step,:,:],np.exp(2.0*np.memmap.sum(dtau*lyaploc_clv[tni:tn+n_step,np.newaxis,:], axis =0))))
)
for n_ind, n_step in enumerate(steplengthforsecondorder):
solution[tn,n_ind,:,:] = np.multiply(solution[tn,n_ind,:,:] ,np.exp(np.sum(dtau*lyaploc_clv[tn:tn+n_step,:,np.newaxis], axis = 0)))
if n_ind == 0: continue
else:
full_solution[tn,n_ind,:,:] = np.matmul(CLV[tn+n_step,:,:],solution[tn,n_step,:,:])
growth[tn,n_ind,:]=np.log(np.linalg.norm(full_solution[tn,n_ind,:,:],axis=0))/(n_step*dtau)
if tn % 1 == 0:
np.memmap.flush(solution)
np.memmap.flush(contracted_CLVs)
#np.memmap.flush(normalized_solution)
np.memmap.flush(full_solution)
np.memmap.flush(growth)
print("Saveing results in folder "+savename+".")
del contracted_CLVs
del solution
del full_solution
del growth
del CLV
del lyaploc_clv
del invCLV