added function to produce the finite difference matrices studied in the paper by De Sterck et al., NLAA 2021

Author: danielru

scripts/pseudo-spectrum/figure_5-8-loglog.py

@@ -43,7 +43,7 @@ def trap(z):
 Tend, nslices, maxiter, nfine, ncoarse, tol, epsilon, ndof_f = par.getpar()
 
 nsteps   = [1, 2, 4, 8, 12, 16, 20]
-nsteps   = [1, 10, 25, 50, 75, 100]
+#nsteps   = [10, 100, 1000, 10000, 100000]
 
 ndof_c_v = [16, 24, 30]
 xaxis_f  = np.linspace(0.0, 1.0, ndof_f+1)[0:ndof_f]
@@ -52,6 +52,8 @@ def trap(z):
 if figure==5:
   A_f = get_upwind(ndof_f, dx_f)
   u0fine   = solution_linear(np.zeros(ndof_f), A_f)
+  D = A_f*A_f.H - A_f.H*A_f
+  print("Normality number of the system matrix (this should be zero): %5.3f" % np.linalg.norm(D.todense()))    
 elif figure==6 or figure==7:
   A_f = get_centered(ndof_f, dx_f)
   u0fine   = solution_linear(np.zeros(ndof_f), A_f)
@@ -78,15 +80,15 @@ def trap(z):
   if figure==5:
     A_c = get_upwind(ndof_c, dx_c)
     u0coarse = solution_linear(np.zeros(ndof_c), A_c)
-    filename = 'figure_5.pdf'
+    filename = 'figure_5_mod.pdf'
   elif figure==6:
     A_c = get_centered(ndof_c, dx_c)
     u0coarse = solution_linear(np.zeros(ndof_c), A_c)
-    filename = 'figure_6.pdf'
+    filename = 'figure_6_mod.pdf'
   elif figure==7:
     A_c = get_centered(ndof_c, dx_c)
     u0coarse = solution_linear(np.zeros(ndof_c), A_c)
-    filename = 'figure_7.pdf'    
+    filename = 'figure_7_mod.pdf'    
   elif figure==8:
     u0coarse = solution_dedalus(np.zeros(ndof_c), ndof_c)
     filename = 'figure_8.pdf'
@@ -109,24 +111,21 @@ def trap(z):
     Pmat, Bmat = para.get_parareal_matrix()
     dt_f_v[0,mm] = para.timemesh.slices[0].int_fine.dt
     dt_c_v[0,mm] = para.timemesh.slices[0].int_coarse.dt
-    
     # Sort according to absolute values of R(z) with z = lambda*dt_f
     sort_index = np.argsort(np.abs(Rz(eig_val*dt_f_v[0,mm])))
 
     # Store eigenvalue ndof_c+1. first "truncated" EV
-    A_f_eig[nn,mm] = np.flip((eig_val[sort_index]))[ndof_c+1]
-  
+    A_f_eig[nn,mm] = np.flip((eig_val[sort_index]))[ndof_c]
     ### Parareal iteration: y^k+1 = Pmat*y^k + Bmat*b
     norm_l2[nn,mm] = np.linalg.norm(Pmat.todense(), 2)
-    
 
 rcParams['figure.figsize'] = 2.5, 2.5
 fs = 8
 ms = 4
 fig = plt.figure(1)
-plt.loglog(dt_f_v[0,:], norm_l2[0,:]-np.abs(np.exp(np.multiply(A_f_eig[0,:],dt_f_v[0,:]))), 'bo-', label='m='+str(ndof_c_v[0]), markersize=ms)
-plt.loglog(dt_f_v[0,:], norm_l2[1,:]-np.abs(np.exp(np.multiply(A_f_eig[1,:],dt_f_v[0,:]))), 'rx-', label='m='+str(ndof_c_v[1]), markersize=ms)
-plt.loglog(dt_f_v[0,:], norm_l2[2,:]-np.abs(np.exp(np.multiply(A_f_eig[2,:],dt_f_v[0,:]))), 'cd-', label='m='+str(ndof_c_v[2]), markersize=ms)
+plt.loglog(dt_f_v[0,:], norm_l2[0,:]-np.abs(np.exp(np.multiply(A_f_eig[0,:],np.multiply(dt_f_v[0,:],nsteps[:])))), 'bo-', label='m='+str(ndof_c_v[0]), markersize=ms)
+plt.loglog(dt_f_v[0,:], norm_l2[1,:]-np.abs(np.exp(np.multiply(A_f_eig[1,:],np.multiply(dt_f_v[0,:],nsteps[:])))), 'rx-', label='m='+str(ndof_c_v[1]), markersize=ms)
+plt.loglog(dt_f_v[0,:], norm_l2[2,:]-np.abs(np.exp(np.multiply(A_f_eig[2,:],np.multiply(dt_f_v[0,:],nsteps[:])))), 'cd-', label='m='+str(ndof_c_v[2]), markersize=ms)
 plt.legend(loc='best', bbox_to_anchor=(0.5, 0.5), fontsize=fs, prop={'size':fs-2}, handlelength=3)
 plt.xlabel(r'$\delta t = \Delta t$', fontsize=fs)
 plt.ylabel(r'$|| \mathbf{E} ||_2 - | \exp(\lambda_{m+1} \delta t) |$', fontsize=fs)

src/get_matrix.py

@@ -28,3 +28,35 @@ def get_diffusion(n,h):
     col[-1] = 1.0
     A       = (1.0/h**2)*spla.circulant(col)
     return sp.csc_matrix(A)
+  
+def get_desterck(n, h, p):
+  col = np.zeros(n)
+  if p==2:
+    col[0] = 3.0
+    col[1] = -4.0
+    col[2] = 1-0
+    A = -1.0/(2.0*h)*spla.circulant(col)
+  elif p==3:
+    col[0] = 3.0
+    col[1] = -6.0
+    col[2] = 1.0
+    col[-1] = 2.0
+    A = -1.0/(6.0*h)*spla.circulant(col)
+  elif p==4:
+    col[0] = 10.0
+    col[1] = -18.0
+    col[2] = 6.0
+    col[3] = -1.0
+    col[-1] = 3.0
+    A = -1.0/(12.0*h)*spla.circulant(col)
+  elif p==5:
+    col[0] = 20.0
+    col[1] = -60.0
+    col[2] = 15.0
+    col[3] = -2.0
+    col[-1] = 30.0
+    col[-2] = -3.0
+    A = -1.0/(60.0*h)*spla.circulant(col)
+  else:
+    sys.exit("Do not have coefficients for requested p")
+  return sp.csc_matrix(A)

tests/test_getmatrix.py

@@ -0,0 +1,64 @@
+import sys
+sys.path.append('./src')
+
+import pytest
+import numpy as np
+from get_matrix import get_upwind, get_centered, get_diffusion, get_desterck
+
+class TestClass:
+    
+    def f(x):
+      return np.sin(2.0*np.pi*x)
+    
+    def minus_fx(x):
+      return -2.0*np.pi*np.cos(2.0*np.pi*x)
+  
+    def fxx(x):
+      return -4.0*np.pi**2*np.sin(2.0*np.pi*x)  
+    
+    def setUp(self):
+      self.ndof = 256
+      self.xaxis = np.linspace(0.0, 1.0, self.ndof+1)[0:self.ndof]
+      self.dx = self.xaxis[1] - self.xaxis[0]
+      self.u  = TestClass.f(self.xaxis)
+      self.ux = TestClass.minus_fx(self.xaxis)
+      self.uxx = TestClass.fxx(self.xaxis)
+      
+    '''
+    The following tests confirm that the finite difference provide an approximation that is within
+    some fixed tolerance to the analytic solution
+    '''
+    
+    def testUpwindIsAccurate(self):
+      self.setUp()
+      A = get_upwind(self.ndof, self.dx)    
+      ux_fd = A*self.u
+      err = np.linalg.norm(self.ux - ux_fd, np.inf)
+      assert err < 1e-1, "Upwind finite difference seems inaccurate"
+      
+    def testCenteredIsAccurate(self):
+      self.setUp()
+      A = get_centered(self.ndof, self.dx)      
+      ux_fd = A*self.u
+      err = np.linalg.norm(self.ux - ux_fd, np.inf)
+      assert err < 1e-2, "Centered finite difference seems inaccurate"
+        
+    def testDiffusionIsAccurate(self):
+      self.setUp()
+      A = get_diffusion(self.ndof, self.dx)      
+      uxx_fd = A*self.u
+      err = np.linalg.norm(self.uxx - uxx_fd, np.inf)
+      assert err < 1e-2, "Centered finite difference for diffusion seems inaccurate"
+        
+    def testDeSterckIsAccurate(self):
+      self.setUp()
+      tols = [2e-2, 8e-6, 2e-7, 1e-9]
+      for p in range(4):
+          A = get_desterck(self.ndof, self.dx, p+2)
+          ux_fd = A*self.u
+          err = np.linalg.norm(self.ux - ux_fd, np.inf)
+          assert err < tols[p], ("Order %i of De Sterck et al finite difference seems inaccurate" % (p+2))
+      
+            
+      
+