Skip to content

Commit

Permalink
order memmap
Browse files Browse the repository at this point in the history
  • Loading branch information
Sebastian Schubert committed Sep 22, 2017
1 parent 72311f9 commit 6e9c253
Show file tree
Hide file tree
Showing 6 changed files with 82 additions and 84 deletions.
24 changes: 12 additions & 12 deletions lorenz96_compute_secondary_instabilities.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,14 +32,14 @@
'dimY' : 10,
'RescaledY' : False,
'expname' : 'secondaryinstabilities_1layer',
'time' : np.arange(0,500,0.01),
'time' : np.arange(0,500,0.1),
'spinup' : 100,
'2lay' : False
}


testzeroclv=True
steplengthforsecondorder = np.arange(0,100,1)
steplengthforsecondorder = np.arange(0,200,1)
hs=[1.0]#, 0.5] # , 0.0625, 0.125 , 0.25 , 0.5 , 1. ]
precision='float64'
for paraL96,h in product([paraL96_1lay],hs):
Expand All @@ -62,18 +62,18 @@


t = paraL96['time']
CLV = np.memmap(savename+'/CLV.dat',mode='r',shape=(len(t),dimN,M),dtype='float64')
lyaploc_clv = np.memmap(savename+'/lyaploc_clv',mode='r',shape=(len(t),M),dtype='float64')
invCLV = np.memmap(savename+'/invCLV.dat',mode='r',shape=(len(t),dimN,M),dtype='float64')
CLV = np.memmap(savename+'/CLV.dat',mode='r',shape=(len(t),dimN,M),dtype='float64', order = 'F')
lyaploc_clv = np.memmap(savename+'/lyaploc_clv',mode='r',shape=(len(t),M),dtype='float64', order = 'F')
invCLV = np.memmap(savename+'/invCLV.dat',mode='r',shape=(len(t),dimN,M),dtype='float64', order = 'F')



lensteparray = len(steplengthforsecondorder)
contracted_CLVs = np.memmap(savename+'/contracted_clvs.dat',mode='w+',shape=(len(t),dimN,M),dtype=precision) #(final_clv,init_clv)
solution = np.memmap(savename+'/solution.dat',mode='w+',shape=(len(t),lensteparray,dimN,M),dtype=precision)
full_solution = np.memmap(savename+'/full_solution.dat',mode='w+',shape=(len(t),lensteparray,dimN,M),dtype=precision)
#normalized_solution = np.memmap(savename+'/normalized_solution.dat',mode='w+',shape=(len(t),lensteparray,dimN,M),dtype=precision)
growth = np.memmap(savename+'/growth.dat',mode='w+',shape=(len(t),lensteparray,M),dtype=precision)
contracted_CLVs = np.memmap(savename+'/contracted_clvs.dat',mode='w+',shape=(len(t),dimN,M),dtype=precision, order = 'F') #(final_clv,init_clv)
solution = np.memmap(savename+'/solution.dat',mode='w+',shape=(len(t),lensteparray,dimN,M),dtype=precision, order = 'F')
full_solution = np.memmap(savename+'/full_solution.dat',mode='w+',shape=(len(t),lensteparray,dimN,M),dtype=precision, order = 'F')
#normalized_solution = np.memmap(savename+'/normalized_solution.dat',mode='w+',shape=(len(t),lensteparray,dimN,M),dtype=precision, order = 'F')
growth = np.memmap(savename+'/growth.dat',mode='w+',shape=(len(t),lensteparray,M),dtype=precision, order = 'F')
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,:,:])
Expand All @@ -92,13 +92,13 @@
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)))
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 % 400 == 0:
if tn % 1 == 0:
np.memmap.flush(solution)
np.memmap.flush(contracted_CLVs)
#np.memmap.flush(normalized_solution)
Expand Down
28 changes: 14 additions & 14 deletions lorenz96_secondary_instabilities.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -34,7 +34,7 @@
'dimY' : 10,
'RescaledY' : False,
'expname' : 'secondaryinstabilities_1layer',
'time' : np.arange(0,500,0.01),
'time' : np.arange(0,500,0.1),
'spinup' : 100,
'2lay' : False
}
Expand All @@ -56,7 +56,7 @@
dimN = paraL96['dimX']

t = paraL96['time']
dt = np.mean(np.diff(t))
dt = 0.01#np.mean(np.diff(t))

savename=paraL96['expname']+"_h_"+str(h)
spinup = paraL96['spinup']
Expand All @@ -66,17 +66,17 @@


if not os.path.exists(savename): os.mkdir(savename)
CLV = np.memmap(savename+'\CLV.dat',mode='w+',shape=(len(t),dimN,M),dtype='float64')
BLV = np.memmap(savename+'\BLV.dat',mode='w+',shape=(len(t),dimN,M),dtype='float64')
R = np.memmap(savename+'\R.dat',mode='w+',shape=(len(t),dimN,M),dtype='float64')
lyapmean_blv = np.memmap(savename+'\lyapmean_blv.dat',mode='w+',shape=(M),dtype='float64')
lyapmean_clv = np.memmap(savename+'\lyapmean_clv.dat',mode='w+',shape=(M),dtype='float64')
lyaploc_clv = np.memmap(savename+'\lyaploc_clv',mode='w+',shape=(len(t),M),dtype='float64')
lyaploc_blv = np.memmap(savename+'\lyaploc_blv',mode='w+',shape=(len(t)-1,M),dtype='float64')
np.save(savename+'\t',t)
trajectory = np.memmap(savename+'\trajectory.dat',mode='w+',shape=(len(t),dimN),dtype='float64')
if testzeroclv: tendency = np.memmap(savename+'\tendency.dat',mode='w+',shape=(len(t),dimN),dtype='float64')
if testzeroclv: tendcorr = np.memmap(savename+'\tendcorr.dat',mode='w+',shape=(len(t)),dtype='float64')
CLV = np.memmap(savename+'/CLV.dat',mode='w+',shape=(len(t),dimN,M),dtype='float64', order = 'F')
BLV = np.memmap(savename+'/BLV.dat',mode='w+',shape=(len(t),dimN,M),dtype='float64', order = 'F')
R = np.memmap(savename+'/R.dat',mode='w+',shape=(len(t),dimN,M),dtype='float64', order = 'F')
lyapmean_blv = np.memmap(savename+'/lyapmean_blv.dat',mode='w+',shape=(M),dtype='float64', order = 'F')
lyapmean_clv = np.memmap(savename+'/lyapmean_clv.dat',mode='w+',shape=(M),dtype='float64', order = 'F')
lyaploc_clv = np.memmap(savename+'/lyaploc_clv',mode='w+',shape=(len(t),M),dtype='float64', order = 'F')
lyaploc_blv = np.memmap(savename+'/lyaploc_blv',mode='w+',shape=(len(t)-1,M),dtype='float64', order = 'F')
np.save(savename+'/t',t)
trajectory = np.memmap(savename+'/trajectory.dat',mode='w+',shape=(len(t),dimN),dtype='float64', order = 'F')
if testzeroclv: tendency = np.memmap(savename+'/tendency.dat',mode='w+',shape=(len(t),dimN),dtype='float64', order = 'F')
if testzeroclv: tendcorr = np.memmap(savename+'/tendcorr.dat',mode='w+',shape=(len(t)),dtype='float64', order = 'F')


# Compute the exponents
Expand Down Expand Up @@ -164,7 +164,7 @@
np.save(savename+"/paraL96",paraL96)
np.save(savename+"/h",h)

invCLV = np.memmap(savename+'\invCLV.dat',mode='w+',shape=(len(t),dimN,M),dtype='float64')
invCLV = np.memmap(savename+'/invCLV.dat',mode='w+',shape=(len(t),dimN,M),dtype='float64')

for tn, (ts,te) in enumerate(zip(t[0:-1],t[1:])):
invCLV[tn,:,:]=np.linalg.inv(CLV[tn,:,:])
Expand Down
13 changes: 8 additions & 5 deletions lorenz96_test_secondary_instabilities.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -55,19 +55,19 @@
integrator = 'rk4'

# first test clv
epsilons=10.0**np.arange(-8,1.5,0.1,dtype = precision)
intsteps=np.arange(0,100,1)
epsilons=10.0**np.arange(-2,1.6,0.1,dtype = precision)
intsteps=np.arange(0,200,1)


CLVs=[1,8,14,30] # ,13, 30]#np.arange(1,2,1)
timeintervall = range(2000,4000,500)
timeintervall = range(2000,4000+1,300)
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[()]
dt = np.mean(np.diff(paraL96['time']))
dt = 0.01#np.mean(np.diff(paraL96['time']))
# M number exponents
if paraL96['2lay']:
M = paraL96['dimX'] + paraL96['dimX']*paraL96['dimY'] # -1 full spectrum
Expand Down Expand Up @@ -123,12 +123,14 @@
print(' eps: '+str(epsilon))
for nc, clv in enumerate(CLVs):
field.x['back']=trajectory[tn,:]+epsilon*CLV[tn,:,clv-1]
#print(0,np.log10(epsilon),np.log10(np.linalg.norm(epsilon*CLV[tn,:,clv-1])))

for step, (a,b) in enumerate(zip(intsteps[0:-1],intsteps[1:])):
stepsize=b
#print('g')
for counting in np.arange(0,b-a):
field.integrate_back(dtau,dt =dt)
#print('h')
#print('counting')
nonlinear_prediction = field.x['back'] - trajectory[tn+stepsize,:]
firstorder_prediction = CLV[tn+stepsize,:,clv-1]*np.exp(dtau*np.memmap.sum(lyaploc_clv[tn:tn+stepsize,clv-1]))*epsilon
secondorder_prediction = full_solution[tn,stepsize,:,clv-1]*epsilon**2.0
Expand All @@ -141,6 +143,7 @@
+CLV[tn+stepsize,:,clv-1]
*np.exp(dtau*np.memmap.sum(lyaploc_clv[tn:tn+stepsize,clv-1])))

#print(step+1,np.log10(epsilon),np.log10(np.linalg.norm(nonlinear_prediction - firstorder_prediction)),np.log10(np.linalg.norm(nonlinear_prediction - firstandsecondorder_prediction)))

correlation[nc][i,en,step] = corr(nonlinear_prediction,firstorder_prediction)

Expand Down
24 changes: 12 additions & 12 deletions plot_results_secondary_instabilities.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,20 +11,20 @@



CLV = np.memmap(savename+'/CLV.dat',mode='r',shape=(len(t),dimN,M,1),dtype='float64')
BLV = np.memmap(savename+'/BLV.dat',mode='r',shape=(len(t),dimN,M,1),dtype='float64')
R = np.memmap(savename+'/R.dat',mode='r',shape=(len(t),dimN,M,1),dtype='float64')
lyapmean_blv = np.memmap(savename+'/lyapmean_blv.dat',mode='r',shape=(M,1),dtype='float64')
lyapmean_clv = np.memmap(savename+'/lyapmean_clv.dat',mode='r',shape=(M,1),dtype='float64')
lyaploc_clv = np.memmap(savename+'/lyaploc_clv',mode='r',shape=(len(t),M,1),dtype='float64')
lyaploc_blv = np.memmap(savename+'/lyaploc_blv',mode='r',shape=(len(t)-1,M,1),dtype='float64')
CLV = np.memmap(savename+'/CLV.dat',mode='r',shape=(len(t),dimN,M,1),dtype='float64', order = 'F')
BLV = np.memmap(savename+'/BLV.dat',mode='r',shape=(len(t),dimN,M,1),dtype='float64', order = 'F')
R = np.memmap(savename+'/R.dat',mode='r',shape=(len(t),dimN,M,1),dtype='float64', order = 'F')
lyapmean_blv = np.memmap(savename+'/lyapmean_blv.dat',mode='r',shape=(M,1),dtype='float64', order = 'F')
lyapmean_clv = np.memmap(savename+'/lyapmean_clv.dat',mode='r',shape=(M,1),dtype='float64', order = 'F')
lyaploc_clv = np.memmap(savename+'/lyaploc_clv',mode='r',shape=(len(t),M,1),dtype='float64', order = 'F')
lyaploc_blv = np.memmap(savename+'/lyaploc_blv',mode='r',shape=(len(t)-1,M,1),dtype='float64', order = 'F')

precision='float32'
contracted_CLVs = np.memmap(savename+'/contracted_clvs.dat',mode='r',shape=(len(t),dimN,M),dtype=precision) #(final_clv,init_clv)
solution = np.memmap(savename+'/solution.dat',mode='r',shape=(len(steplengthforsecondorder),len(t),dimN,M),dtype=precision)
normalized_solution = np.memmap(savename+'/solution.dat',mode='r',shape=(len(steplengthforsecondorder),len(t),dimN,M),dtype=precision)
growth = np.memmap(savename+'/growth.dat',mode='r',shape=(len(steplengthforsecondorder),len(t),M),dtype=precision)
invCLV = np.memmap(savename+'/invCLV.dat',mode='r',shape=(len(t),dimN,M,1),dtype='float64')
contracted_CLVs = np.memmap(savename+'/contracted_clvs.dat',mode='r',shape=(len(t),dimN,M),dtype='float64', order = 'F') #(final_clv,init_clv)
solution = np.memmap(savename+'/solution.dat',mode='r',shape=(len(steplengthforsecondorder),len(t),dimN,M),dtype='float64', order = 'F')
normalized_solution = np.memmap(savename+'/solution.dat',mode='r',shape=(len(steplengthforsecondorder),len(t),dimN,M),dtype='float64', order = 'F')
growth = np.memmap(savename+'/growth.dat',mode='r',shape=(len(steplengthforsecondorder),len(t),M),dtype='float64', order = 'F')
invCLV = np.memmap(savename+'/invCLV.dat',mode='r',shape=(len(t),dimN,M,1),dtype='float64', order = 'F')



Expand Down
17 changes: 9 additions & 8 deletions plot_results_secondary_instabilities_projections.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,7 +31,7 @@
'dimY' : 10,
'RescaledY' : False,
'expname' : 'secondaryinstabilities_1layer',
'time' : np.arange(0,1000,0.1),
'time' : np.arange(0,100,0.1),
'spinup' : 100,
'2lay' : False
}
Expand All @@ -52,7 +52,7 @@

compute = True #compute array projections

averageintervall=np.arange(2000,8000)
averageintervall=np.arange(1000,4000)

for paraL96,h in product(experiments,hs):

Expand All @@ -74,8 +74,8 @@
CLV = np.memmap(savename+'/CLV.dat',mode='r',shape=(len(paraL96['time']),dimN,M),dtype='float64')

if compute:
solution = np.memmap(savename+'/solution.dat',mode='r',shape=(len(paraL96['time']),len(steplengthforsecondorder),dimN,M),dtype=precision)
full_solution = np.memmap(savename+'/full_solution.dat',mode='r',shape=(len(paraL96['time']),len(steplengthforsecondorder),dimN,M),dtype=precision)
solution = np.memmap(savename+'/solution.dat',mode='r',shape=(len(paraL96['time']),len(steplengthforsecondorder),dimN,M),dtype=precision,order = 'F')
full_solution = np.memmap(savename+'/full_solution.dat',mode='r',shape=(len(paraL96['time']),len(steplengthforsecondorder),dimN,M),dtype=precision, order = 'F')
if compute:
projections=np.zeros((len(steplengthforsecondorder),dimN,M))
correlations=np.zeros((len(steplengthforsecondorder),dimN,M))
Expand All @@ -84,10 +84,11 @@
correlations = np.load(savename+"/correlations.npy")
if compute:
dummy = np.zeros((len(steplengthforsecondorder),solution.shape[2],solution.shape[3]))
for tn in averageintervall:
dummy = dummy + np.abs(solution[tn,:,:,:])
dummy = dummy/len(averageintervall)

# for tn in averageintervall:
# print(tn)
# dummy = dummy + np.abs(solution[tn,:,:,:])
# dummy = dummy/len(averageintervall)
dummy = np.memmap.mean(np.abs(solution[averageintervall,:,:,:]))
projections[1:,:,:]=np.divide(dummy[1:,:,:], np.linalg.norm(dummy[1:,:,:], axis = 1, keepdims=True))


Expand Down
60 changes: 27 additions & 33 deletions plot_test_secondary_instabilities.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -52,7 +52,7 @@
integrator = 'classic'

# first test clv
min_epsilon=10**0
min_epsilon=10**-5
max_epsilon=101


Expand Down Expand Up @@ -96,8 +96,8 @@
normerror = np.memmap(savename+'/normerror_clv'+str(clv)+'.dat',mode='r',shape=(len(timeintervall),len(epsilons),len(intsteps)-1),dtype=precisioncorr,order = 'F')
normerror_1st = np.memmap(savename+'/normerror_1st_clv'+str(clv)+'.dat',mode='r',shape=(len(timeintervall),len(epsilons),len(intsteps)-1),dtype=precisioncorr,order = 'F')
normerrorrel = np.memmap(savename+'/normerrorrel_clv'+str(clv)+'.dat',mode='r',shape=(len(timeintervall),len(epsilons),len(intsteps)-1),dtype=precisioncorr,order = 'F')
normerrorrel_1st = np.memmap(savename+'/normerrorrel_clv'+str(clv)+'.dat',mode='r',shape=(len(timeintervall),len(epsilons),len(intsteps)-1),dtype=precisioncorr,order = 'F')
normerrorrel_2nd = np.memmap(savename+'/normerrorrel_clv'+str(clv)+'.dat',mode='r',shape=(len(timeintervall),len(epsilons),len(intsteps)-1),dtype=precisioncorr,order = 'F')
normerrorrel_1st = np.memmap(savename+'/normerrorrel_1st_clv'+str(clv)+'.dat',mode='r',shape=(len(timeintervall),len(epsilons),len(intsteps)-1),dtype=precisioncorr,order = 'F')
normerrorrel_2nd = np.memmap(savename+'/normerrorrel_2nd_clv'+str(clv)+'.dat',mode='r',shape=(len(timeintervall),len(epsilons),len(intsteps)-1),dtype=precisioncorr,order = 'F')

normnonlin = np.memmap(savename+'/normnonlin_clv'+str(clv)+'.dat',mode='r',shape=(len(timeintervall),len(epsilons),len(intsteps)-1),dtype=precisioncorr,order = 'F')

Expand Down Expand Up @@ -210,37 +210,31 @@



# fig=plt.figure()
# for eps,epsval in enumerate(epsilons[imin:imax-1]):
# plt.plot(np.mean(np.log10(np.abs(normerrorrel[:,imin+eps,:])), axis = 0)/np.log10(epsval),np.mean(np.abs(correlationv2[:,imin+eps,:]),axis =0),label=r'$\log_{10}(\epsilon)='+str(np.log10(epsval))+r'$')
# plt.legend(loc = 'upper left',fontsize=6)
# plt.ylabel('Correlation with both orders vs rel error of both orders')
# plt.xlabel(r'$\log_{\epsilon}\left(||error||\right)$')
# plt.ylim(0,1.0)
# plt.xticks(np.arange(-8,4))
# plt.xlim(-8,4)
# ax= plt.gca()
# ax.spines['right'].set_visible(False)
# ax.spines['top'].set_visible(False)
# fig.savefig(resultsfolder+"/erroranalysis/CLV_"+str(clv)+"_errorgrowthvscorrelation.pdf")
# fig.savefig(resultsfolder+"/erroranalysis/CLV_"+str(clv)+"_errorgrowthvscorrelation.png", dpi =400)
# plt.close(fig)
fig=plt.figure()
for ist,istval in enumerate(intsteps[1::20]):
#print(ist,istval)
plt.plot(epsilons[imin:imax],np.mean(normerrorrel_1st[:,imin:imax,ist], axis = 0),label=r'$\tau = '+str(istval*dtau)+'$')
plt.legend(loc = 'upper left',fontsize=6)
plt.title('Relative error scaling with 1st order; CLV '+str(clv),fontsize=8)
plt.xscale('log')
plt.yscale('log')
plt.grid(True)
fig.savefig(resultsfolder+"/erroranalysis/CLV_"+str(clv)+"_errorscaling_rel_1st.pdf")
fig.savefig(resultsfolder+"/erroranalysis/CLV_"+str(clv)+"_errorscaling_rel_1st.png", dpi =400)
plt.close(fig)

# fig=plt.figure()
# for eps,epsval in enumerate(epsilons[imin:imax-1]):
# plt.plot(np.mean(np.log10(np.abs(normerrorrel_1st[:,imin+eps,:])), axis = 0)/np.log10(epsval),np.mean(np.abs(correlation[:,imin+eps,:]),axis =0),label=r'$\log_{10}(\epsilon)='+str(np.log10(epsval))+r'$')
# plt.legend(loc = 'upper left',fontsize=6)
# plt.ylabel('Correlation with 1st order vs rel error of 1st order')
# plt.xlabel(r'$\log_{\epsilon}\left(||error||\right)$')
# plt.ylim(0,1.0)
# plt.xticks(np.arange(-8,4))
# plt.xlim(-8,4)
# ax= plt.gca()
# ax.spines['right'].set_visible(False)
# ax.spines['top'].set_visible(False)
# fig.savefig(resultsfolder+"/erroranalysis/CLV_"+str(clv)+"_errorgrowthvscorrelation_1st.pdf")
# fig.savefig(resultsfolder+"/erroranalysis/CLV_"+str(clv)+"_errorgrowthvscorrelation_1st.png", dpi =400)
# plt.close(fig)
fig=plt.figure()
for ist,istval in enumerate(intsteps[1::20]):
#print(ist,istval)
plt.plot(epsilons[imin:imax],np.mean(normerrorrel[:,imin:imax,ist], axis = 0),label=r'$\tau = '+str(istval*dtau)+'$')
plt.legend(loc = 'upper left',fontsize=6)
plt.title('Relative error scaling with 1st and 2nd order; CLV '+str(clv),fontsize=8)
plt.xscale('log')
plt.yscale('log')
plt.grid(True)
fig.savefig(resultsfolder+"/erroranalysis/CLV_"+str(clv)+"_errorscaling_rel_1stand2nd.pdf")
fig.savefig(resultsfolder+"/erroranalysis/CLV_"+str(clv)+"_errorscaling_rel_1stand2nd.png", dpi =400)
plt.close(fig)
#
# fig=plt.figure()
# for eps,epsval in enumerate(epsilons[imin:imax-1]):
Expand Down

0 comments on commit 6e9c253

Please sign in to comment.