Merge branch 'master' of https://github.com/Parallel-in-Time/pyParareal

Author: Daniel Ruprecht

README.md

@@ -29,12 +29,12 @@ How can I reproduce Figures from the publication?
  - Fig. 7 --> scripts/plot_dispersion.py with artificial_coarse=1 and artifical_coarse=2
  - Fig. 8 --> example/run.py with artificial_coarse=2
  - Fig. 9 --> scripts/plot_svd_vs_waveno_different_G.py
- - Fig. 8 --> scripts/plot_dispersion.py with artifical_fine=1
- - Fig. 9 --> example/run.py with artifical_fine=1
- - Fig. 10 -> scripts/plot_dispersion.py with ncoarse=2
- - Fig. 11 -> scripts/plot_dispersion.py with Tend=64
- - Fig. 12 -> scripts/plot_svd_vs_P.py
- - Fig. 13 -> scripts/plot_conv_vs_waveno.py
+ - Fig. 10 -> scripts/plot_dispersion.py with artifical_fine=1
+ - Fig. 11 -> example/run.py with artifical_fine=1
+ - Fig. 12 -> scripts/plot_dispersion.py with ncoarse=2
+ - Fig. 13 -> scripts/plot_dispersion.py with Tend=64 and Nsamples=120
+ - Fig. 14 -> scripts/plot_svd_vs_P.py
+ - Fig. 15 -> scripts/plot_conv_vs_waveno.py
  - Fig. 14 -> scripts/plot_svd_vs_waveno.py 
 
 

scripts/plot_conv_vs_waveno.py

@@ -14,18 +14,14 @@ def solve_omega(R):
 
 if __name__ == "__main__":
 
-    Nx = 200
-    x = np.linspace(0,20,Nx+1,endpoint=False)
-    x = x[0:Nx]
-
     Nk    = 6
     k_vec = np.linspace(0, np.pi, Nk+1, endpoint=False)
     k_vec = k_vec[1:]
     k_vec = [k_vec[0], k_vec[1], k_vec[-1]]
     Tend    = 16.0    
     nslices = 16
     U_speed = 1.0
-    nu      = 0.1
+    nu      = 0.0
     ncoarse = 1
     nfine   = 1
 
@@ -54,64 +50,34 @@ def solve_omega(R):
       if abs(stab_fine[0,0]-stab_ex)>1e-14:
         print "WARNING: Fine method is not the exact integrator..."    
 
-      y_start = np.exp(1j*k_vec[kk]*x)
-      y_ex    = stab_ex*y_start
-
-      y_coarse= stab_coarse[0,0]*y_start
-      y_fine  = stab_fine[0,0]*y_start
-
       omega_fine = solve_omega(stab_fine[0,0])
       phase_fine = omega_fine.real/k_vec[kk]
       amp_fine   = np.exp(omega_fine.imag)
 
       for n in range(1,nslices+1):
         stab_para = para.get_parareal_stab_function(n)
-        y_para  = stab_para[0,0]*y_start
-        err[kk,n-1] = np.linalg.norm(y_para - y_fine, np.inf)/np.linalg.norm(y_fine, np.inf)
-        omega_para = solve_omega(stab_para[0,0])
-        phase_para = omega_para.real/k_vec[kk]
-        amp_para = np.exp(omega_para.imag)
-        err_phase[kk,n-1] = abs(phase_para - phase_fine)/abs(phase_fine)
-        err_amp[kk,n-1]   = abs(amp_para - amp_fine)/abs(amp_fine)
+        err[kk,n-1] = abs(stab_para[0,0] - stab_fine[0,0])
 
     fs = 8
-    rcParams['figure.figsize'] = 2.5, 2.5
+    rcParams['figure.figsize'] = 3.54, 3.54
     fig = plt.figure()
     iter_v = range(1,nslices)
     assert np.max(err[:,-1])<1e-10, ("For at least one wavenumber, Parareal did not fully converge for niter=nslices. Error: %5.3e" % np.max(err[:,-1]))
     plt.semilogy(iter_v, err[0,0:-1], 'b-o', label=(r"$\kappa$=%4.2f" % k_vec[0]), markersize=fs/2)
-    plt.semilogy(iter_v, err[0,0]*np.power(svds[0], iter_v), 'b--')
+    #plt.semilogy(iter_v, err[0,0]*np.power(svds[0], iter_v), 'b--')
 
     plt.semilogy(iter_v, err[1,0:-1], 'r-s', label=(r"$\kappa$=%4.2f" % k_vec[1]), markersize=fs/2)
-    plt.semilogy(iter_v, err[1,0]*np.power(svds[1], iter_v), 'r--')
+    #plt.semilogy(iter_v, err[1,0]*np.power(svds[1], iter_v), 'r--')
 
     plt.semilogy(iter_v, err[2,0:-1], 'g-x', label=(r"$\kappa$=%4.2f" % k_vec[2]), markersize=fs/2)
-    plt.semilogy(iter_v, err[2,0]*np.power(svds[2], iter_v), 'g--')
-
-    #plt.semilogy(iter_v, err[3,0:-1], 'k-d', label=(r"$\kappa$=%4.2f" % k_vec[3]), markersize=fs/2)
-    #plt.semilogy(iter_v, err[3,0]*np.power(svds[3], iter_v), 'k--')
+    #plt.semilogy(iter_v, err[2,0]*np.power(svds[2], iter_v), 'g--')
 
     plt.legend(loc='lower left', fontsize=fs, prop={'size':fs-2}, handlelength=3)
     plt.xlabel('Parareal Iteration $k$', fontsize=fs)
     plt.ylabel('Parareal defect', fontsize=fs)    
-    plt.xticks(np.arange(iter_v[0], iter_v[-1], 2), fontsize=fs)
+    plt.xticks(np.arange(iter_v[0], nslices, 2), fontsize=fs)
     plt.yticks([1e-14, 1e-12, 1e-10, 1e-8, 1e-6, 1e-4, 1e-2, 1e0, 1e2], fontsize=fs)
     filename = 'parareal-conv-waveno.pdf'
     plt.gcf().savefig(filename, bbox_inches='tight')
     call(["pdfcrop", filename, filename])
-
-    #fig = plt.figure()
-    #plt.semilogy(iter_v, err_phase[0,0:-1], 'b-o', label=(r"$\kappa$=%4.2f" % k_vec[0]), markersize=fs/2)
-    #plt.semilogy(iter_v, err_phase[1,0:-1], 'r-s', label=(r"$\kappa$=%4.2f" % k_vec[1]), markersize=fs/2)
-    #plt.semilogy(iter_v, err_phase[2,0:-1], 'g-x', label=(r"$\kappa$=%4.2f" % k_vec[2]), markersize=fs/2)
-    #plt.semilogy(iter_v, err_phase[3,0:-1], 'k-d', label=(r"$\kappa$=%4.2f" % k_vec[3]), markersize=fs/2)
-    #plt.legend(loc='lower left', fontsize=fs, prop={'size':fs}, handlelength=3)
-
-    #fig = plt.figure()
-    #plt.semilogy(iter_v, err_amp[0,0:-1], 'b-o', label=(r"$\kappa$=%4.2f" % k_vec[0]), markersize=fs/2)
-    #plt.semilogy(iter_v, err_amp[1,0:-1], 'r-s', label=(r"$\kappa$=%4.2f" % k_vec[1]), markersize=fs/2)
-    #plt.semilogy(iter_v, err_amp[2,0:-1], 'g-x', label=(r"$\kappa$=%4.2f" % k_vec[2]), markersize=fs/2)
-    #plt.semilogy(iter_v, err_amp[3,0:-1], 'k-d', label=(r"$\kappa$=%4.2f" % k_vec[3]), markersize=fs/2)
-    #plt.legend(loc='lower left', fontsize=fs, prop={'size':fs}, handlelength=3)
-
     #plt.show()

scripts/plot_svd_vs_P.py

@@ -49,7 +49,7 @@
           para = parareal(0.0, float(nslices_v[i]), nslices_v[i], intexact, impeuler, nfine, ncoarse, 0.0, 1, u0)
           svds[j,i] = para.get_max_svd(ucoarse=ucoarse)
 
-    rcParams['figure.figsize'] = 2.5, 2.5
+    rcParams['figure.figsize'] = 3.54, 3.54
     fs = 8
     fig  = plt.figure()
     plt.plot(nslices_v, svds[0,:], 'b-o', label=(r"$\kappa$=%4.2f" % waveno_v[0]), markersize=fs/2, markevery=(1,6))
@@ -62,4 +62,4 @@
     plt.legend(loc='lower right', fontsize=fs, prop={'size':fs-2}, handlelength=3)
     filename = 'parareal-svd-vs-p.pdf'
     plt.gcf().savefig(filename, bbox_inches='tight')
-    call(["pdfcrop", filename, filename])
+    call(["pdfcrop", filename, filename])