newsest

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
+
+

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,:,:])

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']

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,32 +43,33 @@
            '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
 
 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)