finalizing scripts to generate paper figures

Author: danielru

scripts/pseudo-spectrum/figure_5_6.py renamed to scripts/pseudo-spectrum/figure_1-4.py

@@ -15,32 +15,44 @@
 try:
   figure      =  int(sys.argv[1]) # 5 generates figure_5, 6 generates figure_6
 except:
-  print("No or wrong command line argument provided, creating figure 5. Use 5 or 6 as command line argument.")
-  figure = 5
+  print("No or wrong command line argument provided, creating figure 1. Use 1, 2, 3 or 4 as command line argument.")
+  figure = 1
 
 
-if figure==5:
+if figure==1:
     A    = get_upwind(nx, h)
-    filename='figure_5.pdf'
-elif figure==6:
+    filename='figure_1.pdf'
+elif figure==2:
     A    = get_centered(nx, h)
-    filename='figure_6.pdf'
+    filename='figure_2.pdf'
+elif figure==3:
+    A    = get_upwind(nx, h)
+    filename='figure_3.pdf'
+elif figure==4:
+    A   = get_centered(nx, h)
+    filename='figure_4.pdf'    
 else:
-    sys.exit("Figure needs to be 5 or 6")
+    sys.exit("Figure needs to be 1, 2, 3 or 4")
 
-eigs, dummy = LA.eig(A.todense())
+eigs_space, dummy = LA.eig(A.todense())
 svds        = LA.svdvals(A.todense())
 
-A_exp                = LA.expm(A*dt)
-eigs_exp, dummy      = LA.eig(A_exp.todense())
+if figure==1 or figure==2:
+  A_exp                = LA.expm(A*dt)
+  eigs_time, dummy      = LA.eig(A_exp.todense())
+  mylabel = '$\exp(\Delta t A)$'
+else:
+  A_ie = LA.inv(np.eye(nx) - dt*A)
+  eigs_time, dummy = LA.eig(A_ie)
+  mylabel = '$(I - \Delta t A)^{-1}$'
 
 rcParams['figure.figsize'] = 2.5, 2.5
 fs = 8
 ms = 4
 fig = plt.figure(1)
-plt.plot(0*eigs.real, np.linspace(-2.1,2.1,nx), 'k', linewidth=2)
-plt.plot(eigs.real, eigs.imag, 'bo', markersize=ms, markerfacecolor='b', label='$\lambda(A)$')
-plt.plot(eigs_exp.real, eigs_exp.imag, 'rs', markersize=ms, markerfacecolor='r', label='$\exp(\lambda(A))$')
+plt.plot(0*eigs_space.real, np.linspace(-2.1,2.1,nx), 'k', linewidth=2)
+plt.plot(eigs_space.real, eigs_space.imag, 'bo', markersize=ms, markerfacecolor='b', label='$A$')
+plt.plot(eigs_time.real, eigs_time.imag, 'rs', markersize=ms, markerfacecolor='r', label=mylabel)
 circle = plt.Circle((0.0,0.0), 1.0, color='k', fill=False)
 fig.gca().add_patch(circle)
 plt.xlim([-2.1, 2.1])

scripts/pseudo-spectrum/plot_norm_vs_dt.py renamed to scripts/pseudo-spectrum/figure_5-8.py

@@ -26,16 +26,16 @@
 try:
   figure      =  int(sys.argv[1]) # 1 generates figure_1, 2 generates figure_2
 except:
-  print("No or wrong command line argument provided, creating figure 11. Use 11, 12, 13 or 14 as command line argument.")
-  figure = 11
-assert 11<= figure <= 14, "Figure should be 11, 12, 13 or 14"
+  print("No or wrong command line argument provided, creating figure 5. Use 5, 6, 7 or 8 as command line argument.")
+  figure = 5
+assert 5<= figure <= 8, "Figure should be 5, 6, 7 or 8"
 
-if figure==11 or figure==12:
+if figure==5 or figure==6:
   par = parameter(dedalus = False)
-elif figure==13 or figure==14:
+elif figure==7 or figure==8:
   par = parameter(dedalus = True)
 else:
-  sys.exit("Figure should be 11, 12, 13 or 14")
+  sys.exit("This should have been caught above")
 
 Tend, nslices, maxiter, nfine, ncoarse, tol, epsilon, ndof_f = par.getpar()
 
@@ -44,16 +44,16 @@
 xaxis_f  = np.linspace(0.0, 1.0, ndof_f+1)[0:ndof_f]
 dx_f     = xaxis_f[1] - xaxis_f[0]
 
-if figure==11:
+if figure==5:
   A_f = get_upwind(ndof_f, dx_f)
   u0fine   = solution_linear(np.zeros(ndof_f), A_f)
-elif figure==12:
+elif figure==6 or figure==7:
   A_f = get_centered(ndof_f, dx_f)
   u0fine   = solution_linear(np.zeros(ndof_f), A_f)
-elif figure==13:
+elif figure==8:
   u0fine = solution_dedalus(np.zeros(ndof_f), ndof_f)
 else:
-  sys.exit("Figure can only have values 11, 12 or 13")    
+  sys.exit("This should have been caught above")
 
 norm_l2  = np.zeros((3,np.size(nsteps)))
 norm_inf = np.zeros((3,np.size(nsteps)))
@@ -66,30 +66,34 @@
   xaxis_c = np.linspace(0.0, 1.0, ndof_c+1)[0:ndof_c]
   dx_c    = xaxis_c[1] - xaxis_c[0]
 
-  if figure==11:
+  if figure==5:
     A_c = get_upwind(ndof_c, dx_c)
     u0coarse = solution_linear(np.zeros(ndof_c), A_c)
-    filename = 'figure_11.pdf'
-  elif figure==12:
+    filename = 'figure_5.pdf'
+  elif figure==6:
     A_c = get_centered(ndof_c, dx_c)
     u0coarse = solution_linear(np.zeros(ndof_c), A_c)
-    filename = 'figure_12.pdf'
-  elif figure==13:
+    filename = 'figure_6.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'    
+  elif figure==8:
     u0coarse = solution_dedalus(np.zeros(ndof_c), ndof_c)
-    filename = 'figure_13.pdf'
+    filename = 'figure_8.pdf'
   else:
-    sys.exit("Problem can only have values 1, 2 or 3")    
+    sys.exit("Value of figure should be")    
 
   for mm in range(np.size(nsteps)):
 
-    if figure==11:
+    if figure==5 or figure==7:
       para     = parareal(0.0, Tend, nslices, impeuler, impeuler, nsteps[mm], nsteps[mm], tol, maxiter, u0fine, u0coarse)
-    elif figure==12:
-      para     = parareal(0.0, Tend, nslices, trapezoidal, trapezoidal, nsteps[mm], nsteps[mm], tol, maxiter, u0fine, u0coarse)
-    elif figure==13:
+    elif figure==6:
+      para     = parareal(0.0, Tend, nslices, trapezoidal, trapezoidal, nsteps[mm], nsteps[mm], tol, maxiter, u0fine, u0coarse)    
+    elif figure==8:
      para     = parareal(0.0, Tend, nslices, integrator_dedalus, integrator_dedalus, nsteps[mm], nsteps[mm], tol, maxiter, u0fine, u0coarse)    
     else:
-      quit()
+      sys.exit("Value of figure should be")    
     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