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

Author: Daniel Ruprecht

README.md

@@ -1,2 +1,43 @@
-# pyParareal
-Python Parareal code for rapid protoyping
+pyParareal
+============
+
+Python Parareal code to analyse Parareal's wave propagation characteristics.
+
+Attribution
+-----------
+You can freely use and reuse this code in line with the BSD license. 
+If you use it (or parts of it) for a publication, please cite
+
+@unpublished{Ruprecht2017,   
+  author = {Ruprecht, Daniel},    
+  title = {{Wave propagation characteristics of Parareal}},    
+  year = {2017}    
+}
+
+(this will be updated once published).
+
+[![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.228104.svg)](https://doi.org/10.5281/zenodo.228104)
+
+How can I reproduce Figures from the publication?
+-----------------
+
+ - Fig. 1 --> scripts/plot_dispersion.py with different parameters
+ - Fig. 2 --> example/run.py
+ - Fig. 3 --> scripts/plot_ieuler_dispersion.py
+ - Fig. 4 --> scripts/plot_dispersion.py with approx_symbol=True
+ - Fig. 5 --> scripts/plot_dispersion.py with artificial_coarse=1 or 2
+ - Fig. 6 --> example/run.py with artificial_coarse=2
+ - Fig. 7 --> 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. 14 -> scripts/plot_svd_vs_waveno.py 
+
+
+Who do I talk to?
+-----------------
+
+This code is written by [Daniel Ruprecht](http://www.parallelintime.org/groups/leeds.html).

example/run.py

@@ -0,0 +1,163 @@
+# Parts of this code are taken from the PseudoSpectralPython tutorial.
+#
+# Author: D. Ketcheson
+# https://github.com/ketch/PseudoSpectralPython/blob/master/PSPython_01-linear-PDEs.ipynb
+#
+# The MIT License (MIT)
+#
+# Copyright (c) 2015 David Ketcheson
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in all
+# copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+# SOFTWARE.
+
+import sys
+sys.path.append('../src')
+import numpy as np
+import scipy.sparse as sp
+import scipy.linalg as la
+
+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 matplotlib.pyplot as plt
+from subprocess import call
+from pylab import rcParams
+
+from joblib import Parallel, delayed
+
+# advection speed
+uadv = 1.0
+
+# diffusivity parameter
+nu = 0.0
+
+# Spatial grid
+m = 64  # Number of grid points in space
+L = 4.0  # Width of spatial domain
+
+# select coarse integrator:
+# 0 = normal backward Euler method
+# 1 = artificially constructed method with phase error from backward Euler and exact amplification factor
+# 2 = artificially constructed method with exact phase and amplification factor from backward Euler
+artificial_coarse = 0
+
+artificial_fine = 0
+
+# Grid points
+x = np.linspace(0, L, m, endpoint=False)
+# Grid spacing
+dx = x[1]-x[0]
+
+# Temporal grid
+tmax  = 16  # Final time
+#N    = 100  # number grid points in time
+#dt   = tmax/float(N)
+nproc = int(tmax)
+
+Kiter_v = [5, 10, 15]
+
+xi = np.fft.fftfreq(m)*m*2*np.pi/L  # Wavenumber "grid"
+# (this is the order in which numpy's FFT gives the frequencies)
+
+# define spatial operator matrix
+D = -1j*xi*uadv - nu*xi**2
+
+# Initial data
+sig   = 1.0
+u     = np.exp(-(x-0.5*L)**2/sig**2)
+u     = np.exp(2.0*np.pi*1j*x)
+uhat0 = np.fft.fft(u)
+
+yend = np.zeros((3,m), dtype='complex')
+
+### Because the exact integrator can only deal with scalar problems, we will need to solve
+### to initial value problem for each wave number independently
+### ...we use the JobLib module to speed up the computational
+def run_parareal(uhat, D, k):
+  sol = solution_linear(np.asarray([[uhat]]), np.asarray([[D]]))
+  para = parareal(tstart=0.0, tend=tmax, nslices=nproc, fine=intexact, coarse=impeuler, nsteps_fine=10, nsteps_coarse=1, tolerance=0.0, iter_max=k, u0 = sol)
+  stab_coarse = para.timemesh.slices[0].get_coarse_update_matrix(sol)
+  stab_ex     = np.exp(D)
+  
+  if artificial_coarse==2:
+    stab_tailor = abs(stab_coarse[0,0])*np.exp(1j*np.angle(stab_ex)) # exact phase speed
+  elif artificial_coarse==1:
+    stab_tailor = abs(stab_ex)*np.exp(1j*np.angle(stab_coarse[0,0])) # exact amplification factor
+
+  if not artificial_coarse == 0:
+    stab_tailor = sp.csc_matrix(np.array([stab_tailor], dtype='complex'))
+    para = parareal(tstart=0.0, tend=tmax, nslices=nproc, fine=intexact, coarse=stab_tailor, nsteps_fine=10, nsteps_coarse=1, tolerance=0.0, iter_max=k, u0 = sol)
+
+  if artificial_fine==1:
+    assert artificial_coarse==0, "Using artifical coarse and fine propagators together is not implemented and probably not working correctly"
+    stab_fine = abs(stab_ex)*np.exp(1j*np.angle(stab_coarse[0,0]))
+    stab_fine = sp.csc_matrix(np.array([stab_fine], dtype='complex'))
+    # Must use nfine=1 in this case
+    para = parareal(tstart=0.0, tend=tmax, nslices=nproc, fine=stab_fine, coarse=impeuler, nsteps_fine=1, nsteps_coarse=1, tolerance=0.0, iter_max=k, u0 = sol)
+
+  para.run()
+  temp = para.get_last_end_value()
+  return temp.y[0,0]
+
+for k in range(0,3):
+  temp = Parallel(n_jobs=2, verbose=5)(delayed(run_parareal)(uhat0[n],D[n], Kiter_v[k]) for n in range(m))
+  yend[k,:] = temp
+
+rcParams['figure.figsize'] = 2.5, 2.5
+fs = 8
+fig1 = plt.figure()
+plt.plot(x, np.fft.ifft(yend[0,:]).real,  '-s', color='r', linewidth=1.5, label='Parareal k='+str(Kiter_v[0]), markevery=(1,6), mew=1.0, markersize=fs/2)
+plt.plot(x, np.fft.ifft(yend[1,:]).real,  '-d', color='r', linewidth=1.5, label='Parareal k='+str(Kiter_v[1]), markevery=(3,6), mew=1.0, markersize=fs/2)
+plt.plot(x, np.fft.ifft(yend[2,:]).real,  '-x', color='r', linewidth=1.5, label='Parareal k='+str(Kiter_v[2]), markevery=(5,6), mew=1.0, markersize=fs/2)
+plt.plot(x, np.fft.ifft(uhat0).real, '--', color='k', linewidth=1.5, label='Exact')
+plt.xlim([x[0], x[-1]])
+plt.ylim([-.2, 1.4])
+plt.xlabel('x', fontsize=fs, labelpad=0.25)
+plt.ylabel('u', fontsize=fs, labelpad=0.5)
+fig1.gca().tick_params(axis='both', labelsize=fs)
+plt.legend(loc='upper left', fontsize=fs, prop={'size':fs-2})
+#plt.show()
+filename = 'parareal-gauss-peak.pdf'
+plt.gcf().savefig(filename, bbox_inches='tight')
+call(["pdfcrop", filename, filename])
+
+xi = np.fft.fftshift(xi)
+xi = xi[m/2:m]/m
+uhat0     = np.fft.fftshift(uhat0)
+yend[0,:] = np.around(np.fft.fftshift(yend[0,:]), 10)
+yend[1,:] = np.around(np.fft.fftshift(yend[1,:]), 10)
+yend[2,:] = np.around(np.fft.fftshift(yend[2,:]), 10)
+
+fig2 = plt.figure()
+plt.semilogy(xi, np.absolute(uhat0[m/2:m])/m, '--', color='k', linewidth=1.5, label='Exact')
+plt.semilogy(xi, np.absolute(yend[0,m/2:m])/m, '-s', color='r', linewidth=1.5, label='Parareal k='+str(Kiter_v[0]), markevery=(1,3), mew=1.0, markersize=fs/2)
+plt.semilogy(xi, np.absolute(yend[1,m/2:m])/m, '-d', color='r', linewidth=1.5, label='Parareal k='+str(Kiter_v[1]), markevery=(2,3), mew=1.0, markersize=fs/2)
+plt.semilogy(xi, np.absolute(yend[2,m/2:m])/m, '-x', color='r', linewidth=1.5, label='Parareal k='+str(Kiter_v[2]), markevery=(3,3), mew=1.0, markersize=fs/2)
+plt.xlabel('Wave number', fontsize=fs, labelpad=0.25)
+plt.ylabel(r'abs($\hat{u}$)')
+plt.xticks([0.0, 0.1, 0.3], fontsize=fs)
+plt.xlim([0.0, 0.3])
+plt.yticks([1e0, 1e-3, 1e-6, 1e-9, 1e-12], fontsize=fs)
+plt.legend(loc='upper right', fontsize=fs, prop={'size':fs-2})
+filename = 'parareal-gauss-peak-spectrum.pdf'
+plt.gcf().savefig(filename, bbox_inches='tight')
+call(["pdfcrop", filename, filename])

scripts/plot_conv_vs_waveno.py

@@ -25,7 +25,7 @@ def solve_omega(R):
     Tend    = 16.0    
     nslices = 16
     U_speed = 1.0
-    nu      = 0.0
+    nu      = 0.1
     ncoarse = 1
     nfine   = 1
 
@@ -78,7 +78,7 @@ def solve_omega(R):
     rcParams['figure.figsize'] = 2.5, 2.5
     fig = plt.figure()
     iter_v = range(1,nslices)
-    assert np.max(err[:,-1])<1e-14, "For at least one wavenumber, Parareal did not fully converge for niter=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--')