script to plot maximum singular value versus coarse time step

Author: Daniel Ruprecht

scripts/plot_dispersion.py

@@ -4,6 +4,7 @@
 from parareal import parareal
 from impeuler import impeuler
 from intexact import intexact
+from trapezoidal import trapezoidal
 from special_integrator import special_integrator
 from solution_linear import solution_linear
 import numpy as np
@@ -22,7 +23,7 @@ def solve_omega(R):
 
 def findroots(R, n):
   assert abs(n - float(int(n)))<1e-14, "n must be an integer or a float equal to an integer"
-  p = np.zeros(n+1, dtype='complex')
+  p = np.zeros(int(n)+1, dtype='complex')
   p[-1] = -R
   p[0]  = 1.0
   return np.roots(p)
@@ -60,8 +61,8 @@ def normalise(R, T, target):
     for i in range(0,np.size(k_vec)):
 
       symb = -(1j*U_speed*k_vec[i] + nu*k_vec[i]**2)
-      symb_coarse = symb
-#      symb_coarse = -(1.0/dx)*(1.0 - np.exp(-1j*k_vec[i]*dx))
+#      symb_coarse = symb
+      symb_coarse = -(1.0/dx)*(1.0 - np.exp(-1j*k_vec[i]*dx))
 
       # Solution objects define the problem
       u0      = solution_linear(u0_val, np.array([[symb]],dtype='complex'))
@@ -94,24 +95,24 @@ def normalise(R, T, target):
       #
 
       # for stab = r*exp(i*theta), r defines the amplitude factor and theta the phase speed
-      stab_tailor = abs(stab_coarse[0,0])*np.exp(1j*np.angle(stab_ex)) # exact phase speed
+      # stab_tailor = abs(stab_coarse[0,0])*np.exp(1j*np.angle(stab_ex)) # exact phase speed
       # stab_tailor = abs(stab_ex)*np.exp(1j*np.angle(stab_coarse[0,0])) # exact amplification factor
 
       # coarse method with exact phase but amplification factor corresponding to a single large backward Euler step
-      stab_coarse_limit = 1.0/(1.0 - Tend*symb_coarse)
-      stab_tailor = abs(stab_coarse_limit)*np.exp(1j*np.angle(stab_ex)) # exact phase speed, massively diffusive amplification factor
+      #stab_coarse_limit = 1.0/(1.0 - Tend*symb_coarse)
+      #stab_tailor = abs(stab_coarse_limit)*np.exp(1j*np.angle(stab_ex)) # exact phase speed, massively diffusive amplification factor
 
       # stab_tailor = abs(stab_ex)*np.exp(1j*np.angle(stab_ex)) ## for testing
       # stab_tailor = abs(stab_coarse[0,0])*np.exp(1j*np.angle(stab_coarse[0,0])) ## for testing
 
       # Re-Create the parareal object to be used in the remainder
-      stab_tailor = sparse.csc_matrix(np.array([stab_tailor], dtype='complex'))
-      para = parareal(0.0, Tend, nslices, intexact, stab_tailor, nfine, ncoarse, 0.0, niter_v[0], u0)
+      # stab_tailor = sparse.csc_matrix(np.array([stab_tailor], dtype='complex'))
+      # para = parareal(0.0, Tend, nslices, intexact, stab_tailor, nfine, ncoarse, 0.0, niter_v[0], u0)
 
       #################################################
 
       # Compute Parareal phase velocity and amplification factor
-      #svds[0,i]         = para.get_max_svd(ucoarse=ucoarse)        
+      svds[0,i]         = para.get_max_svd(ucoarse=ucoarse)        
       for jj in range(0,3):
         stab_para = para.get_parareal_stab_function(k=niter_v[jj], ucoarse=ucoarse)
 

scripts/plot_sine.py

@@ -32,7 +32,7 @@
 
     err = np.zeros(nslices)
     para_show = np.zeros((3,Nx))
-    niter_show = [5, 10, 15]
+    niter_show = [2, 4, 6]
 
     symb      = -(1j*U_speed*k + nu*k**2)
     u0_val    = np.array([[1.0]], dtype='complex')
@@ -41,7 +41,7 @@
 
     stab_ex   = np.exp(-1j*U_speed*k*Tend)*np.exp(-nu*k**2*Tend)
 
-    stab_coarse   = para.timemesh.slices[0].get_coarse_update_matrix(u0)
+    stab_coarse = para.timemesh.slices[0].get_coarse_update_matrix(u0)
     stab_coarse = stab_coarse**nslices
 
     stab_fine = para.timemesh.slices[0].get_fine_update_matrix(u0)
@@ -67,7 +67,7 @@
     plt.plot(x, para_show[0,:], 'r--+', label='Parareal k='+str(niter_show[0]), markevery=(5, 20), markersize=fs/2, mew=1)
     plt.plot(x, para_show[1,:], 'r:s', label='Parareal k='+str(niter_show[1]), markevery=(10,20),  markersize=fs/2, mew=1)
     plt.plot(x, para_show[2,:], 'r-o', label='Parareal k='+str(niter_show[2]), markevery=(15,20),  markersize=fs/2, mew=1)
-    plt.plot(x, y_ex.real,      'g--', label='Fine')
+    #plt.plot(x, y_ex.real,      'g--', label='Fine')
     plt.legend(loc='lower left', fontsize=fs, prop={'size':fs-2}, handlelength=3)
     plt.title((r'$\kappa$ = %4.2f' % k), fontsize=fs)
     plt.ylim([-2.5, 1.5])

scripts/plot_svd_vs_dt.py

@@ -0,0 +1,58 @@
+import sys
+sys.path.append('../src')
+
+from parareal import parareal
+from impeuler import impeuler
+from intexact import intexact
+from trapezoidal import trapezoidal
+from special_integrator import special_integrator
+from solution_linear import solution_linear
+import numpy as np
+import scipy.sparse as sparse
+import math
+
+from pylab import rcParams
+import matplotlib.pyplot as plt
+from matplotlib.patches import Polygon
+from subprocess import call
+import sympy
+from pylab import rcParams
+
+if __name__ == "__main__":
+
+    Tend     = 16.0
+    nslices  = int(Tend) # Make sure each time slice has length 1
+    U_speed  = 1.0
+    nu       = 0.0
+    ncoarse_v = range(1,11)
+    nfine    = 10
+    niter_v  = [5, 10, 15]
+    dx       = 1.0
+    Nsamples = 60
+    u0_val     = np.array([[1.0]], dtype='complex')
+
+    k_vec = np.linspace(0.0, np.pi, Nsamples+1, endpoint=False)
+    k_vec = k_vec[1:]
+    waveno = k_vec[-1]
+
+    svds = np.zeros((1,np.size(ncoarse_v)))
+    dt_v = np.zeros((1,np.size(ncoarse_v)))
+
+    symb = -(1j*U_speed*waveno + nu*waveno**2)
+    symb_coarse = symb
+#    symb_coarse = -(1.0/dx)*(1.0 - np.exp(-1j*waveno*dx))
+
+    # Solution objects define the problem
+    u0      = solution_linear(u0_val, np.array([[symb]],dtype='complex'))
+    ucoarse = solution_linear(u0_val, np.array([[symb_coarse]],dtype='complex'))
+
+    for i in range(0,np.size(ncoarse_v)):
+      para = parareal(0.0, Tend, nslices, intexact, impeuler, nfine, ncoarse_v[i], 0.0, niter_v[2], u0)
+      dt_v[0,i] = Tend/float(ncoarse_v[i]*nslices)
+      svds[0,i]         = para.get_max_svd(ucoarse=ucoarse)        
+
+    rcParams['figure.figsize'] = 7.5, 7.5
+    fs = 8
+    fig  = plt.figure()
+    plt.plot(dt_v[0,:], svds[0,:])
+    plt.show()

scripts/plot_wave_time.py

@@ -60,7 +60,7 @@ def normalise(R, T, target):
     u0      = solution_linear(u0_val, np.array([[symb]],dtype='complex'))
     ucoarse = solution_linear(u0_val, np.array([[symb_coarse]],dtype='complex'))
 
-    para = parareal(0.0, Tend, nslices, intexact, trapezoidal, nfine, ncoarse, 0.0, niter_v[0], u0)
+    para = parareal(0.0, Tend, nslices, intexact, impeuler, nfine, ncoarse, 0.0, niter_v[0], u0)
 
     # get update matrix for imp Euler over one slice
     stab_fine   = para.timemesh.slices[0].get_fine_update_matrix(u0)