added option to generate pseudo-spectrum plot and convergence plot using the high(er) order upwind finite difference stencils investigated in the paper by DeSterck et al., 2021, NLAA

Author: danielru

scripts/pseudo-spectrum/figure_13-16.py

@@ -12,7 +12,7 @@
 from solution_dedalus import solution_dedalus
 from solution_linear import solution_linear
 
-from get_matrix import get_upwind, get_centered, get_diffusion
+from get_matrix import get_upwind, get_centered, get_diffusion, get_desterck
 from impeuler import impeuler
 from intexact import intexact
 from trapezoidal import trapezoidal
@@ -29,9 +29,9 @@
 except:
   print("No or wrong command line argument provided, creating figure 13. Use 13, 14, 15 or 16 as command line argument.")
   figure = 5
-assert 13<= figure <= 16 or figure==0, "Figure should be 13, 14, 15 or 16"
+assert 13<= figure <= 16 or figure==0 or figure==-1, "Figure should be 13, 14, 15 or 16"
 
-if figure==13 or figure==14 or figure==0:
+if figure==13 or figure==14 or figure==0 or figure==-1:
   par = parameter(dedalus = False)
   ndof_c   = 24
 elif figure==15: 
@@ -91,7 +91,19 @@
   para     = parareal(0.0, Tend, nslices, impeuler, impeuler, nfine, ncoarse, tol, maxiter, u0fine, u0coarse)  
   filename = 'figure_00.pdf'    
   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()))   
+  print("Normality number of the system matrix (this should be zero): %5.3f" % np.linalg.norm(D.todense()))  
+elif figure==-1:
+  p = 5
+  xaxis_f = np.linspace(0.0, 1.0, ndof_f+1)[0:ndof_f]
+  dx_f    = xaxis_f[1] - xaxis_f[0]
+  xaxis_c = np.linspace(0.0, 1.0, ndof_c+1)[0:ndof_c]
+  dx_c = xaxis_c[1] - xaxis_c[0]  
+  A_f = get_desterck(ndof_f, dx_f, p)
+  A_c = get_desterck(ndof_c, dx_c, p)
+  u0fine   = solution_linear(np.zeros(ndof_f), A_f)
+  u0coarse = solution_linear(np.zeros(ndof_c), A_c)  
+  para     = parareal(0.0, Tend, nslices, impeuler, impeuler, nfine, ncoarse, tol, maxiter, u0fine, u0coarse)  
+  filename = ('figure_conv_desterck_%i.pdf' % p)  
 else:
   sys.exit("Wrong value for figure")
 

scripts/pseudo-spectrum/figure_9-12.py

@@ -9,7 +9,7 @@
 from solution_linear import solution_linear
 from integrator_dedalus import integrator_dedalus
 from solution_dedalus import solution_dedalus
-from get_matrix import get_upwind, get_centered, get_diffusion
+from get_matrix import get_upwind, get_centered, get_diffusion, get_desterck
 import numpy as np
 import scipy.sparse as sparse
 from scipy.linalg import svdvals
@@ -27,9 +27,9 @@
 except:
   print("No or wrong command line argument provided, creating figure 9. Use 9, 10, 11 or 12 as command line argument.")
   figure = 5
-assert (9<= figure <= 12) or figure==0, "Figure should be 9, 10, 11 or 12"
+assert (9<= figure <= 12) or figure==0 or figure==-1, "Figure should be 9, 10, 11 or 12"
 
-if figure==9 or figure==10 or figure==0:
+if figure==9 or figure==10 or figure==0 or figure==-1:
   par = parameter(dedalus = False)
   ndof_c   = 24
 elif figure==11: 
@@ -88,6 +88,18 @@
   u0coarse = solution_linear(np.zeros(ndof_c), A_c)  
   para     = parareal(0.0, Tend, nslices, impeuler, impeuler, nfine, ncoarse, tol, maxiter, u0fine, u0coarse)  
   filename = 'figure_00.pdf'
+elif figure==-1:
+  p = 5
+  xaxis_f = np.linspace(0.0, 1.0, ndof_f+1)[0:ndof_f]
+  dx_f    = xaxis_f[1] - xaxis_f[0]
+  xaxis_c = np.linspace(0.0, 1.0, ndof_c+1)[0:ndof_c]
+  dx_c = xaxis_c[1] - xaxis_c[0]  
+  A_f = get_desterck(ndof_f, dx_f, p)
+  A_c = get_desterck(ndof_c, dx_c, p)
+  u0fine   = solution_linear(np.zeros(ndof_f), A_f)
+  u0coarse = solution_linear(np.zeros(ndof_c), A_c)  
+  para     = parareal(0.0, Tend, nslices, impeuler, impeuler, nfine, ncoarse, tol, maxiter, u0fine, u0coarse)  
+  filename = ('figure_desterck_%i.pdf' % p)
 else:
   sys.exit("Wrong value for figure")