diff --git a/l96.py b/l96.py old mode 100644 new mode 100755 index 7e67cb7..6e1a1b6 --- a/l96.py +++ b/l96.py @@ -23,20 +23,29 @@ def setupL96(para): Im1=(I-1) % N Ip1=(I+1) % N #Ip2=(I+2) % N - + M = N L96 = lambda t,x: (x[Ip1]-x[Im2]) * x[Im1]- x[I] + F - L96Jac = lambda x,t: np.array(np.stack((-x[Im1],x[Ip1]-x[Im2],-1.0*np.ones(N),x[Im1]),axis = 1)) + mydiag = lambda a,k: np.roll(np.diag(a),shift=k,axis=1) + print("hier") + L96Jac = lambda x,t: -np.eye(N) + mydiag(x[Ip1]-x[Im2],-1) + mydiag(x[Im1],1) - mydiag(x[Im1],-2) L96JacV = lambda x,v,t : np.multiply(v[Im2,:],-x[Im1,np.newaxis])+np.multiply(v[Im1,:],x[Ip1,np.newaxis]-x[Im2,np.newaxis])-v[I,:]+np.multiply(v[Ip1,:],x[Im1,np.newaxis]) #L96Full = lambda t,x: np.concatenate((L96(x[0:dimN],t), np.reshape(L96JacV(x[0:dimN],np.reshape(x[dimN:],(dimN,dimN)),t),((dimN)**2,1))) - - #for k in np.arange(0,dimN): - # HessMat[XI[k],XIp1[k],XIm1[k]]=1 - # HessMat[XI[k],XIm1[k],XIp1[k]]=1 - # HessMat[XI[k],XIm2[k],XIm1[k]]=-1 - # HessMat[XI[k],XIm1[k],XIm2[k]]=-1 - HessMat = None - L96JacFull = None + HessMat = np.zeros((N,N,N)) + for k in np.arange(0,dimN): + HessMat[I[k],Ip1[k],Im1[k]]=1 + HessMat[I[k],Im1[k],Ip1[k]]=1 + HessMat[I[k],Im2[k],Im1[k]]=-1 + HessMat[I[k],Im1[k],Im2[k]]=-1 + #HessMat = None + def L96JacFull(x,t): + res = np.zeros((N+N*M,N+N*M)) + res[:N,:N] = L96Jac(x[:N],t) + for ix in range(0,M): + res[N+N*ix:N+N*(ix+1),N+N*ix:N+N*(ix+1)] = L96Jac(x[:N],t) + res[N+N*ix:N+N*(ix+1),:N] = np.sum(np.multiply(HessMat,x[np.newaxis,np.newaxis,N+N*ix:N+N*(ix+1)]),axis=2) + return res + return L96,L96Jac,L96JacV,L96JacFull,dimN def contract_hess_l96_1layer(a,b,c): @@ -433,7 +442,7 @@ def integrate(self,f,y0, delta_t, integrator = None , method = None,jac = None, r.integrate(r.t+delta_t) return r.y,r.t elif integrator == 'classic': - intres = odeint(lambda x,t : f(t,x), y0, [self.step_t,self.step_t+delta_t], mxstep=0, Dfun = jac,rtol=self.rtol,atol=self.atol) + intres = odeint(lambda x,t : f(t,x), y0, [self.step_t,self.step_t+delta_t], Dfun = jac,rtol=self.rtol,atol=self.atol) return intres[-1,:],self.step_t+delta_t elif integrator == 'rk4': tend=self.RK4(f) @@ -447,7 +456,7 @@ def integrate(self,f,y0, delta_t, integrator = None , method = None,jac = None, def integrate_back(self,delta_t,integrator=None,dt= None): if not integrator: integrator=self.defaultintegrator self.xold['back'] = self.x['back'] - self.x['back'], self.step_t = self.integrate(self.tendency, self.x['back'], delta_t,integrator=integrator,dt= dt)# ,max_step = 1,jac=None)#self.jacobian) + self.x['back'], self.step_t = self.integrate(self.tendency, self.x['back'], delta_t,integrator=integrator,dt= dt,jac=self.jacobian) # integrate background and tangent linear at the same time def integrate_all(self, delta_t, normalize = False,integrator=None,dt= None): @@ -462,3 +471,21 @@ def integrate_all(self, delta_t, normalize = False,integrator=None,dt= None): self.step = self.step + 1 if normalize: self.x['lin'] = self.x['lin']/np.linalg.norm(self.x['lin'],axis=0,keepdims=True) + + + + +def test_jac(tend,jac,dim,n=100,time=0,epsrange = 10**np.arange(-8,1.1),tendtimeord=0,jactimeord=1): + correlation = np.zeros((len(epsrange))) + for ieps,eps in enumerate(epsrange): + for i in range(0,n): + testa = np.random.rand(dim) + testb = testa+eps*np.random.rand(dim) + if tendtimeord==0: tenddiff = (tend(time,testb) - tend(time,testa))/eps + else: tenddiff = (tend(testb,time) - tend(testa,time))/eps + if jactimeord==0: jacdiff = np.matmul(jac(time,testa),(testb-testa)/eps) + else: jacdiff = np.matmul(jac(testa,time),(testb-testa)/eps) + correlation[ieps] = correlation[ieps] + np.corrcoef(tenddiff,jacdiff)[0,1]/np.float(n) + return correlation + + \ No newline at end of file diff --git a/lorenz96.py b/lorenz96.py old mode 100644 new mode 100755 diff --git a/lorenz96_anglesclvs.py b/lorenz96_anglesclvs.py old mode 100644 new mode 100755 diff --git a/lorenz96_compute_2_secondary_instabilities.py b/lorenz96_compute_2_secondary_instabilities.py old mode 100644 new mode 100755 diff --git a/lorenz96_compute_secondary_instabilities.py b/lorenz96_compute_secondary_instabilities.py old mode 100644 new mode 100755 index 471d1ef..9eefb2c --- a/lorenz96_compute_secondary_instabilities.py +++ b/lorenz96_compute_secondary_instabilities.py @@ -1,4 +1,4 @@ -#!/home/user/anaconda3/bin/python + #!/home/user/anaconda3/bin/python import numpy as np import os import l96 @@ -20,7 +20,8 @@ 'expname' : 'secondaryinstabilities_2layer', 'time' : np.arange(0,500,0.1), 'spinup' : 100, - '2lay' : True + '2lay' : True, + 'integrator': 'classic' } paraL96_1lay = {'F1' : 10, @@ -32,14 +33,15 @@ 'dimY' : 10, 'RescaledY' : False, 'expname' : 'secondaryinstabilities_1layer', - 'time' : np.arange(0,500,0.1), + 'time' : np.arange(0,500,0.05), 'spinup' : 100, - '2lay' : False + '2lay' : False, + 'integrator': 'classic' } testzeroclv=True -steplengthforsecondorder = np.arange(0,200,1) +steplengthforsecondorder = np.arange(0,400,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): @@ -62,18 +64,18 @@ t = paraL96['time'] - 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') + 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 = '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') + 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,:,:]) diff --git a/lorenz96_condnb.py b/lorenz96_condnb.py old mode 100644 new mode 100755 diff --git a/lorenz96_condnb_max.py b/lorenz96_condnb_max.py old mode 100644 new mode 100755 diff --git a/lorenz96_definterun.py b/lorenz96_definterun.py old mode 100644 new mode 100755 diff --git a/lorenz96_fft.py b/lorenz96_fft.py old mode 100644 new mode 100755 diff --git a/lorenz96_gabrielesetup.py b/lorenz96_gabrielesetup.py old mode 100644 new mode 100755 diff --git a/lorenz96_gabrielesetup_convergence.py b/lorenz96_gabrielesetup_convergence.py old mode 100644 new mode 100755 diff --git a/lorenz96_gabrielesetup_memmap.py b/lorenz96_gabrielesetup_memmap.py old mode 100644 new mode 100755 diff --git a/lorenz96_inttest.py b/lorenz96_inttest.py old mode 100644 new mode 100755 diff --git a/lorenz96_secondary_instabilities.py b/lorenz96_secondary_instabilities.py old mode 100644 new mode 100755 index f1d24fa..c9c340a --- a/lorenz96_secondary_instabilities.py +++ b/lorenz96_secondary_instabilities.py @@ -22,21 +22,23 @@ 'expname' : 'secondaryinstabilities_2layer', 'time' : np.arange(0,500,0.1), 'spinup' : 100, - '2lay' : True + '2lay' : True, + 'integrator': 'classic' } -paraL96_1lay = {'F1' : 10, +paraL96_1lay = {'F1' : 8, 'F2' : 0, 'b' : 10, 'c' : 10, 'h' : 1, - 'dimX': 44, + 'dimX': 5, 'dimY' : 10, 'RescaledY' : False, 'expname' : 'secondaryinstabilities_1layer', - 'time' : np.arange(0,500,0.1), + 'time' : np.arange(0,100,0.1), 'spinup' : 100, - '2lay' : False + '2lay' : False, + 'integrator': 'classic' } @@ -66,17 +68,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', 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') + CLV = np.memmap(savename+'/CLV.dat',mode='w+',shape=(len(t),dimN,M),dtype='float64', order = 'C') + BLV = np.memmap(savename+'/BLV.dat',mode='w+',shape=(len(t),dimN,M),dtype='float64', order = 'C') + R = np.memmap(savename+'/R.dat',mode='w+',shape=(len(t),dimN,M),dtype='float64', order = 'C') + lyapmean_blv = np.memmap(savename+'/lyapmean_blv.dat',mode='w+',shape=(M),dtype='float64', order = 'C') + lyapmean_clv = np.memmap(savename+'/lyapmean_clv.dat',mode='w+',shape=(M),dtype='float64', order = 'C') + lyaploc_clv = np.memmap(savename+'/lyaploc_clv',mode='w+',shape=(len(t),M),dtype='float64', order = 'C') + lyaploc_blv = np.memmap(savename+'/lyaploc_blv',mode='w+',shape=(len(t)-1,M),dtype='float64', order = 'C') 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') + trajectory = np.memmap(savename+'/trajectory.dat',mode='w+',shape=(len(t),dimN),dtype='float64', order = 'C') + if testzeroclv: tendency = np.memmap(savename+'/tendency.dat',mode='w+',shape=(len(t),dimN),dtype='float64', order = 'C') + if testzeroclv: tendcorr = np.memmap(savename+'/tendcorr.dat',mode='w+',shape=(len(t)),dtype='float64', order = 'C') # Compute the exponents @@ -87,11 +89,10 @@ if paraL96['2lay']: L96,L96Jac,L96JacV,L96JacFull,dimN = l96.setupL96_2layer(paraL96) else: L96,L96Jac,L96JacV,L96JacFull,dimN = l96.setupL96(paraL96) - integrator = 'rk4' - field = l96.GinelliForward(dimN,M,tendfunc = L96, jacfunc = None, jacVfunc = L96JacV,jacfull=None, integrator=integrator) - field.rtol = 1e-12 - field.atol = 1e-12 + field = l96.GinelliForward(dimN,M,tendfunc = L96, jacfunc = L96Jac, jacVfunc = L96JacV,jacfull=L96JacFull, integrator=paraL96['integrator']) + field.rtol = 1e-8 + field.atol = 1e-8 # initialize fields print("\nInitialize ...") @@ -101,7 +102,10 @@ # spinup print("\nSpinup ...") - field.integrate_back(spinup,dt=dt) + for i in range(0,spinup): + field.integrate_back(1,dt=dt) + field.integrate_back(spinup % 1,dt=dt) + field.step_t = 0.0 BLV[0,:,:]=field.x['lin'] diff --git a/lorenz96_simplerun.py b/lorenz96_simplerun.py old mode 100644 new mode 100755 diff --git a/lorenz96_test_secondary_instabilities.py b/lorenz96_test_secondary_instabilities.py old mode 100644 new mode 100755 index 96b121f..45e4b5a --- a/lorenz96_test_secondary_instabilities.py +++ b/lorenz96_test_secondary_instabilities.py @@ -30,7 +30,8 @@ 'expname' : 'secondaryinstabilities_2layer', 'time' : np.arange(0,500,0.1), 'spinup' : 100, - '2lay' : True + '2lay' : True, + 'integrator': 'classic' } paraL96_1lay = {'F1' : 10, @@ -42,9 +43,10 @@ 'dimY' : 10, 'RescaledY' : False, 'expname' : 'secondaryinstabilities_1layer', - 'time' : np.arange(0,500,0.1), + 'time' : np.arange(0,500,0.05), 'spinup' : 100, - '2lay' : False + '2lay' : False, + 'integrator': 'classic' } @@ -52,22 +54,22 @@ hs=[1.0] #, 0.5] # , 0.0625, 0.125 , 0.25 , 0.5 , 1. ] experiments = [paraL96_1lay]#,paraL96_2lay]#,paraL96_1lay] -integrator = 'rk4' +integrator = 'classic' # first test clv -epsilons=10.0**np.arange(-2,1.6,0.1,dtype = precision) -intsteps=np.arange(0,200,1) +epsilons=10.0**np.arange(0,2.1,0.2,dtype = precision) +intsteps=np.arange(0,100,1) -CLVs=[1,8,14,30] # ,13, 30]#np.arange(1,2,1) -timeintervall = range(2000,4000+1,300) +CLVs=[1,2,3,4,5] # ,13, 30]#np.arange(1,2,1) +timeintervall = range(200,600,100) 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 = 0.01#np.mean(np.diff(paraL96['time'])) + dt = 0.005#np.mean(np.diff(paraL96['time'])) # M number exponents if paraL96['2lay']: M = paraL96['dimX'] + paraL96['dimX']*paraL96['dimY'] # -1 full spectrum @@ -94,28 +96,28 @@ normnonlin =[] for clv in CLVs: - correlation.append(np.memmap(savename+'/correlation_only_1_clv'+str(clv)+'.dat',mode='w+',shape=(len(timeintervall),len(epsilons),len(intsteps)-1),dtype=precision,order = 'F')) - correlationv2.append(np.memmap(savename+'/correlation_1and2_clv'+str(clv)+'.dat',mode='w+',shape=(len(timeintervall),len(epsilons),len(intsteps)-1),dtype=precision,order = 'F')) - correlationv3.append(np.memmap(savename+'/correlation_only_2_clv'+str(clv)+'.dat',mode='w+',shape=(len(timeintervall),len(epsilons),len(intsteps)-1),dtype=precision,order = 'F')) - normerror.append(np.memmap(savename+'/normerror_clv'+str(clv)+'.dat',mode='w+',shape=(len(timeintervall),len(epsilons),len(intsteps)-1),dtype=precision,order = 'F')) - normerror_1st.append(np.memmap(savename+'/normerror_1st_clv'+str(clv)+'.dat',mode='w+',shape=(len(timeintervall),len(epsilons),len(intsteps)-1),dtype=precision,order = 'F')) - normerrorrel_1st.append(np.memmap(savename+'/normerrorrel_1st_clv'+str(clv)+'.dat',mode='w+',shape=(len(timeintervall),len(epsilons),len(intsteps)-1),dtype=precision,order = 'F')) - normerrorrel_2nd.append(np.memmap(savename+'/normerrorrel_2nd_clv'+str(clv)+'.dat',mode='w+',shape=(len(timeintervall),len(epsilons),len(intsteps)-1),dtype=precision,order = 'F')) - normerrorrel.append(np.memmap(savename+'/normerrorrel_clv'+str(clv)+'.dat',mode='w+',shape=(len(timeintervall),len(epsilons),len(intsteps)-1),dtype=precision,order = 'F')) - normnonlin.append(np.memmap(savename+'/normnonlin_clv'+str(clv)+'.dat',mode='w+',shape=(len(timeintervall),len(epsilons),len(intsteps)-1),dtype=precision,order = 'F')) - realgrowth.append(np.memmap(savename+'/realgrowth_clv'+str(clv)+'.dat',mode='w+',shape=(len(timeintervall),len(epsilons),len(intsteps)-1),dtype=precision,order = 'F')) + correlation.append(np.memmap(savename+'/correlation_only_1_clv'+str(clv)+'.dat',mode='w+',shape=(len(timeintervall),len(epsilons),len(intsteps)-1),dtype=precision,order = 'C')) + correlationv2.append(np.memmap(savename+'/correlation_1and2_clv'+str(clv)+'.dat',mode='w+',shape=(len(timeintervall),len(epsilons),len(intsteps)-1),dtype=precision,order = 'C')) + correlationv3.append(np.memmap(savename+'/correlation_only_2_clv'+str(clv)+'.dat',mode='w+',shape=(len(timeintervall),len(epsilons),len(intsteps)-1),dtype=precision,order = 'C')) + normerror.append(np.memmap(savename+'/normerror_clv'+str(clv)+'.dat',mode='w+',shape=(len(timeintervall),len(epsilons),len(intsteps)-1),dtype=precision,order = 'C')) + normerror_1st.append(np.memmap(savename+'/normerror_1st_clv'+str(clv)+'.dat',mode='w+',shape=(len(timeintervall),len(epsilons),len(intsteps)-1),dtype=precision,order = 'C')) + normerrorrel_1st.append(np.memmap(savename+'/normerrorrel_1st_clv'+str(clv)+'.dat',mode='w+',shape=(len(timeintervall),len(epsilons),len(intsteps)-1),dtype=precision,order = 'C')) + normerrorrel_2nd.append(np.memmap(savename+'/normerrorrel_2nd_clv'+str(clv)+'.dat',mode='w+',shape=(len(timeintervall),len(epsilons),len(intsteps)-1),dtype=precision,order = 'C')) + normerrorrel.append(np.memmap(savename+'/normerrorrel_clv'+str(clv)+'.dat',mode='w+',shape=(len(timeintervall),len(epsilons),len(intsteps)-1),dtype=precision,order = 'C')) + normnonlin.append(np.memmap(savename+'/normnonlin_clv'+str(clv)+'.dat',mode='w+',shape=(len(timeintervall),len(epsilons),len(intsteps)-1),dtype=precision,order = 'C')) + realgrowth.append(np.memmap(savename+'/realgrowth_clv'+str(clv)+'.dat',mode='w+',shape=(len(timeintervall),len(epsilons),len(intsteps)-1),dtype=precision,order = 'C')) - CLV = np.memmap(savename+'/CLV.dat',mode='r',shape=(len(paraL96['time']),dimN,M),dtype='float64') - lyaploc_clv = np.memmap(savename+'/lyaploc_clv',mode='r',shape=(len(paraL96['time']),M),dtype='float64') - trajectory = np.memmap(savename+'/trajectory.dat',mode='r',shape=(len(paraL96['time']),dimN),dtype='float64') - full_solution = np.memmap(savename+'/full_solution.dat',mode='r',shape=(len(paraL96['time']),len(steplengthforsecondorder),dimN,M),dtype='float64') + CLV = np.memmap(savename+'/CLV.dat',mode='r',shape=(len(paraL96['time']),dimN,M),dtype='float64',order = 'C') + lyaploc_clv = np.memmap(savename+'/lyaploc_clv',mode='r',shape=(len(paraL96['time']),M),dtype='float64',order = 'C') + trajectory = np.memmap(savename+'/trajectory.dat',mode='r',shape=(len(paraL96['time']),dimN),dtype='float64',order = 'C') + full_solution = np.memmap(savename+'/full_solution.dat',mode='r',shape=(len(paraL96['time']),len(steplengthforsecondorder),dimN,M),dtype='float64',order = 'C') if paraL96['2lay']: L96,L96Jac,L96JacV,L96JacFull,dimN = l96.setupL96_2layer(paraL96) else: L96,L96Jac,L96JacV,L96JacFull,dimN = l96.setupL96(paraL96) - field = l96.GinelliForward(dimN,M,tendfunc = L96, jacfunc = L96Jac, jacVfunc = L96JacV,jacfull=L96JacFull, integrator=integrator) - field.rtol = 1e-12 - field.atol = 1e-12 + field = l96.GinelliForward(dimN,M,tendfunc = L96, jacfunc = L96Jac, jacVfunc = L96JacV,jacfull=L96JacFull, integrator=paraL96['integrator']) + field.rtol = 1e-8 + field.atol = 1e-8 for i,tn in enumerate(timeintervall): print(tn) diff --git a/lorenz96_traditionalrun.py b/lorenz96_traditionalrun.py old mode 100644 new mode 100755 diff --git a/plot_condnb.py b/plot_condnb.py old mode 100644 new mode 100755 diff --git a/plot_results.py b/plot_results.py old mode 100644 new mode 100755 diff --git a/plot_results_definiterun.py b/plot_results_definiterun.py old mode 100644 new mode 100755 diff --git a/plot_results_gabrielesetup.py b/plot_results_gabrielesetup.py old mode 100644 new mode 100755 diff --git a/plot_results_gabrielesetup_convergence.py b/plot_results_gabrielesetup_convergence.py old mode 100644 new mode 100755 diff --git a/plot_results_gabrielesetup_memmap.py b/plot_results_gabrielesetup_memmap.py old mode 100644 new mode 100755 diff --git a/plot_results_secondary_instabilities.py b/plot_results_secondary_instabilities.py old mode 100644 new mode 100755 diff --git a/plot_results_secondary_instabilities_convergence.py b/plot_results_secondary_instabilities_convergence.py old mode 100644 new mode 100755 diff --git a/plot_results_secondary_instabilities_lyapunovexponents.py b/plot_results_secondary_instabilities_lyapunovexponents.py old mode 100644 new mode 100755 diff --git a/plot_results_secondary_instabilities_projections.py b/plot_results_secondary_instabilities_projections.py old mode 100644 new mode 100755 diff --git a/plot_results_simplerun.py b/plot_results_simplerun.py old mode 100644 new mode 100755 diff --git a/plot_results_v2.py b/plot_results_v2.py old mode 100644 new mode 100755 diff --git a/plot_test_secondary_instabilities.py b/plot_test_secondary_instabilities.py old mode 100644 new mode 100755 index 9f85eca..d6e355b --- a/plot_test_secondary_instabilities.py +++ b/plot_test_secondary_instabilities.py @@ -25,7 +25,8 @@ 'expname' : 'secondaryinstabilities_2layer', 'time' : np.arange(0,500,0.1), 'spinup' : 100, - '2lay' : True + '2lay' : True, + 'integrator': 'classic' } paraL96_1lay = {'F1' : 10, @@ -39,7 +40,8 @@ 'expname' : 'secondaryinstabilities_1layer', 'time' : np.arange(0,500,0.1), 'spinup' : 100, - '2lay' : False + '2lay' : False, + 'integrator': 'classic' } @@ -52,7 +54,7 @@ integrator = 'classic' # first test clv -min_epsilon=10**-5 +min_epsilon=10**-2 max_epsilon=101 @@ -89,21 +91,21 @@ precisioncorr = 'float64' for clv in CLVs: - correlation = np.memmap(savename+'/correlation_only_1_clv'+str(clv)+'.dat',mode='r',shape=(len(timeintervall),len(epsilons),len(intsteps)-1),dtype=precisioncorr,order = 'F') - correlationv2 = np.memmap(savename+'/correlation_1and2_clv'+str(clv)+'.dat',mode='r',shape=(len(timeintervall),len(epsilons),len(intsteps)-1),dtype=precisioncorr,order = 'F') - correlationv3 = np.memmap(savename+'/correlation_only_2_clv'+str(clv)+'.dat',mode='r',shape=(len(timeintervall),len(epsilons),len(intsteps)-1),dtype=precisioncorr,order = 'F') - realgrowth = np.memmap(savename+'/realgrowth_clv'+str(clv)+'.dat',mode='r',shape=(len(timeintervall),len(epsilons),len(intsteps)-1),dtype=precisioncorr,order = 'F') - 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_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') + correlation = np.memmap(savename+'/correlation_only_1_clv'+str(clv)+'.dat',mode='r',shape=(len(timeintervall),len(epsilons),len(intsteps)-1),dtype=precisioncorr,order = 'C') + correlationv2 = np.memmap(savename+'/correlation_1and2_clv'+str(clv)+'.dat',mode='r',shape=(len(timeintervall),len(epsilons),len(intsteps)-1),dtype=precisioncorr,order = 'C') + correlationv3 = np.memmap(savename+'/correlation_only_2_clv'+str(clv)+'.dat',mode='r',shape=(len(timeintervall),len(epsilons),len(intsteps)-1),dtype=precisioncorr,order = 'C') + realgrowth = np.memmap(savename+'/realgrowth_clv'+str(clv)+'.dat',mode='r',shape=(len(timeintervall),len(epsilons),len(intsteps)-1),dtype=precisioncorr,order = 'C') + normerror = np.memmap(savename+'/normerror_clv'+str(clv)+'.dat',mode='r',shape=(len(timeintervall),len(epsilons),len(intsteps)-1),dtype=precisioncorr,order = 'C') + normerror_1st = np.memmap(savename+'/normerror_1st_clv'+str(clv)+'.dat',mode='r',shape=(len(timeintervall),len(epsilons),len(intsteps)-1),dtype=precisioncorr,order = 'C') + normerrorrel = np.memmap(savename+'/normerrorrel_clv'+str(clv)+'.dat',mode='r',shape=(len(timeintervall),len(epsilons),len(intsteps)-1),dtype=precisioncorr,order = 'C') + normerrorrel_1st = np.memmap(savename+'/normerrorrel_1st_clv'+str(clv)+'.dat',mode='r',shape=(len(timeintervall),len(epsilons),len(intsteps)-1),dtype=precisioncorr,order = 'C') + normerrorrel_2nd = np.memmap(savename+'/normerrorrel_2nd_clv'+str(clv)+'.dat',mode='r',shape=(len(timeintervall),len(epsilons),len(intsteps)-1),dtype=precisioncorr,order = 'C') - normnonlin = np.memmap(savename+'/normnonlin_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 = 'C') - CLV = np.memmap(savename+'/CLV.dat',mode='r',shape=(len(paraL96['time']),dimN,M),dtype='float64') - lyaploc_clv = np.memmap(savename+'/lyaploc_clv',mode='r',shape=(len(paraL96['time']),M),dtype='float64') - trajectory = np.memmap(savename+'/trajectory.dat',mode='r',shape=(len(paraL96['time']),dimN),dtype='float64') + CLV = np.memmap(savename+'/CLV.dat',mode='r',shape=(len(paraL96['time']),dimN,M),dtype='float64',order = 'C') + lyaploc_clv = np.memmap(savename+'/lyaploc_clv',mode='r',shape=(len(paraL96['time']),M),dtype='float64',order = 'C') + trajectory = np.memmap(savename+'/trajectory.dat',mode='r',shape=(len(paraL96['time']),dimN),dtype='float64',order = 'C') imin = np.abs(epsilons-min_epsilon).argmin() @@ -184,26 +186,42 @@ # plt.close(fig) fig=plt.figure() - for ist,istval in enumerate(intsteps[1::20]): + for ist,istval in enumerate(intsteps[1::40]): #print(ist,istval) plt.plot(epsilons[imin:imax],np.mean(normerror[:,imin:imax,ist], axis = 0),label=r'$\tau = '+str(istval*dtau)+'$') plt.legend(loc = 'upper left',fontsize=6) plt.title('scaling of difference between sum of 1st and 2nd order and nonlinear evolution\n for perturbation along CLV '+str(clv),fontsize=8) plt.xscale('log') plt.yscale('log') - plt.grid(True) + plt.grid(False) + ax=plt.gca() + xlim = ax.get_xlim() + ylim = ax.get_ylim() + for i in np.arange(-10,10,0.5): + plt.plot(epsilons[imin:imax],10**i*epsilons[imin:imax]**2.0*np.min(np.mean(normerror[:,imin:imax,ist], axis = 0)),label='_nolegend_',color='k',linestyle='-',alpha=0.1) + plt.plot(epsilons[imin:imax],10**i*epsilons[imin:imax]**3.0*np.min(np.mean(normerror[:,imin:imax,ist], axis = 0)),label='_nolegend_',color='k',linestyle='-',alpha=0.1) + plt.xlim(xlim[0],xlim[1]) + plt.ylim(ylim[0],ylim[1]) fig.savefig(resultsfolder+"/erroranalysis/CLV_"+str(clv)+"_errorscaling_1stand2nd.pdf") fig.savefig(resultsfolder+"/erroranalysis/CLV_"+str(clv)+"_errorscaling_1stand2nd.png", dpi =400) plt.close(fig) fig=plt.figure() - for ist,istval in enumerate(intsteps[1::20]): + for ist,istval in enumerate(intsteps[1::40]): plt.plot(epsilons[imin:imax],np.mean(normerror_1st[:,imin:imax,ist], axis = 0),label=r'$\tau = '+str(istval*dtau)+'$') plt.legend(loc = 'upper left',fontsize=6) plt.title('scaling of difference between sum of 1st order and nonlinear evolution\n for perturbation along CLV '+str(clv),fontsize=8) plt.xscale('log') plt.yscale('log') - plt.grid(True) + plt.grid(False) + ax=plt.gca() + xlim = ax.get_xlim() + ylim = ax.get_ylim() + for i in np.arange(-10,10,0.5): + plt.plot(epsilons[imin:imax],10**i*epsilons[imin:imax]**2.0*np.min(np.mean(normerror_1st[:,imin:imax,ist], axis = 0)),label='_nolegend_',color='k',linestyle='-',alpha=0.1) + plt.plot(epsilons[imin:imax],10**i*epsilons[imin:imax]**3.0*np.min(np.mean(normerror_1st[:,imin:imax,ist], axis = 0)),label='_nolegend_',color='k',linestyle='-',alpha=0.1) + plt.xlim(xlim[0],xlim[1]) + plt.ylim(ylim[0],ylim[1]) fig.savefig(resultsfolder+"/erroranalysis/CLV_"+str(clv)+"_errorscaling_1st.pdf") fig.savefig(resultsfolder+"/erroranalysis/CLV_"+str(clv)+"_errorscaling_1st.png", dpi =400) plt.close(fig) @@ -211,27 +229,43 @@ fig=plt.figure() - for ist,istval in enumerate(intsteps[1::20]): + for ist,istval in enumerate(intsteps[1::40]): #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) + plt.grid(False) + ax=plt.gca() + xlim = ax.get_xlim() + ylim = ax.get_ylim() + for i in np.arange(-10,10,0.5): + plt.plot(epsilons[imin:imax],10**i*epsilons[imin:imax]**2.0*np.min(np.mean(normerrorrel_1st[:,imin:imax,ist], axis = 0)),label='_nolegend_',color='k',linestyle='-',alpha=0.1) + plt.plot(epsilons[imin:imax],10**i*epsilons[imin:imax]**3.0*np.min(np.mean(normerrorrel_1st[:,imin:imax,ist], axis = 0)),label='_nolegend_',color='k',linestyle='-',alpha=0.1) + plt.xlim(xlim[0],xlim[1]) + plt.ylim(ylim[0],ylim[1]) 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 ist,istval in enumerate(intsteps[1::20]): + for ist,istval in enumerate(intsteps[1::40]): #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) + plt.grid(False) + ax=plt.gca() + xlim = ax.get_xlim() + ylim = ax.get_ylim() + for i in np.arange(-10,10,0.5): + plt.plot(epsilons[imin:imax],10**i*epsilons[imin:imax]**2.0*np.min(np.mean(normerrorrel[:,imin:imax,ist], axis = 0)),label='_nolegend_',color='k',linestyle='-',alpha=0.1) + plt.plot(epsilons[imin:imax],10**i*epsilons[imin:imax]**3.0*np.min(np.mean(normerrorrel[:,imin:imax,ist], axis = 0)),label='_nolegend_',color='k',linestyle='-',alpha=0.1) + plt.xlim(xlim[0],xlim[1]) + plt.ylim(ylim[0],ylim[1]) 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)