Allows to provide a matrix as argument for parareal and timemesh instead of the name of a coarse integrator. The matrix is used as stability matrix to initialise an object of type special_integrator

Author: Daniel Ruprecht

scripts/plot_dispersion.py

@@ -4,6 +4,7 @@
 from parareal import parareal
 from impeuler import impeuler
 from intexact import intexact
+from special_integrator import special_integrator
 from solution_linear import solution_linear
 import numpy as np
 #from scipy.sparse import linalg
@@ -62,10 +63,10 @@ def normalise(R, T, target):
 #      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'))
       ucoarse = solution_linear(u0_val, np.array([[symb_coarse]],dtype='complex'))
-      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)
       stab_fine   = stab_fine
@@ -87,6 +88,10 @@ def normalise(R, T, target):
 
       phase[2,i]      = sol_coarse.real/k_vec[i]
       amp_factor[2,i] = np.exp(sol_coarse.imag)
+
+      # Create the parareal object to be used in the remainder
+      para = parareal(0.0, Tend, nslices, intexact, impeuler, nfine, ncoarse, 0.0, niter_v[0], u0)
+
 
       # Compute Parareal phase velocity and amplification factor
       svds[0,i]         = para.get_max_svd(ucoarse=ucoarse)        

src/special_integrator.py

@@ -13,5 +13,5 @@ def run(self, u0):
     for i in range(0,self.nsteps):
       u0.y = self.stab_function.dot(u0.y)
 
-  def get_update_matrix(self):
-    return self.stab_function
+  def get_update_matrix(self, sol):
+    return self.stab_function**self.nsteps

src/timemesh.py

@@ -1,4 +1,5 @@
 from timeslice import timeslice
+from special_integrator import special_integrator
 import numpy as np
 from scipy.sparse import linalg
 from scipy import sparse
@@ -21,7 +22,10 @@ def __init__(self, tstart, tend, nslices, fine, coarse, nsteps_fine, nsteps_coar
     # ... however, this has not yet been tested!!
     for i in range(0,nslices):
       ts_fine   =   fine(self.timemesh[i], self.timemesh[i+1], nsteps_fine)
-      ts_coarse = coarse(self.timemesh[i], self.timemesh[i+1], nsteps_coarse)
+      if sparse.issparse(coarse):
+        ts_coarse = special_integrator(self.timemesh[i], self.timemesh[i+1], nsteps_coarse, coarse)      
+      else:
+        ts_coarse = coarse(self.timemesh[i], self.timemesh[i+1], nsteps_coarse)
       self.slices.append( timeslice(ts_fine, ts_coarse, tolerance, iter_max) )
 
   # Run the coarse method serially over all slices

tests/test_timemesh.py

@@ -26,6 +26,11 @@ def setUp(self):
   def test_caninstantiate(self):
     tm = timemesh(self.tstart, self.tend, self.nslices, impeuler, impeuler, self.nfine, self.ncoarse, 1e-10, 5)
 
+  # timemesh class can be instantiated with matrix as argument for coarse
+  def test_caninstantiatewithmatrix(self):
+    mat = sparse.eye(1)
+    tm = timemesh(self.tstart, self.tend, self.nslices, impeuler, mat, self.nfine, self.ncoarse, 1e-10, 5)
+
   # initial value can be set for first time slice
   def test_cansetinitial(self):
     tm = timemesh(self.tstart, self.tend, self.nslices, impeuler, impeuler, self.nfine, self.ncoarse, 1e-10, 5)
@@ -105,6 +110,20 @@ def test_runcoarse(self):
     err = abs(u[-1] - tm.get_coarse_value(self.nslices-1).y)
     assert err<1e-12, ("run_coarse and successive application of update matrix does not give identical results - error: %5.3e" % err)
 
+  # with coarse method provided as matrix, run_coarse is callable and provides expected output at very end
+  def test_runcoarsewithmatrix(self):
+    dt = (self.tend - self.tstart)/(float(self.nslices)*float(self.ncoarse))
+    mat = sparse.csc_matrix(np.array([1.0/(1.0 - dt)]))
+    tm = timemesh(self.tstart, self.tend, self.nslices, impeuler, mat, self.nfine, self.ncoarse, 1e-10, 5)
+    tm.run_coarse(self.u0)
+    u = np.array([1.0])
+    Mat = tm.get_coarse_matrix(self.u0)
+    b = np.zeros(self.nslices+1)
+    b[0] = u[0]
+    u = linalg.spsolve(Mat, b)
+    err = abs(u[-1] - tm.get_coarse_value(self.nslices-1).y)
+    assert err<1e-12, ("for coarse provided as matrix, run_coarse and successive application of update matrix does not give identical results - error: %5.3e" % err)
+
   # run_fine is callable and provides expected output at very end
   def test_runfine(self):
     tm = timemesh(self.tstart, self.tend, self.nslices, impeuler, impeuler, self.nfine, self.ncoarse, 1e-10, 5)
@@ -130,6 +149,21 @@ def test_runcoarseintermediate(self):
       err = np.linalg.norm( tm.get_coarse_value(i).y - b, np.inf )
       assert err<2e-12, ("Successive application of update matrix does not reproduce intermediata coarse values generated by run_coarse. Error: %5.3e in slice %2i" % (err, i))
 
+  # with coarse provided as matrix, run_coarse provides expected intermediate values
+  def test_runcoarseintermediatewithmatrix(self):
+    dt = (self.tend - self.tstart)/(float(self.nslices)*float(self.ncoarse))
+    mat = sparse.csc_matrix( np.array([1.0/(1.0 - dt)]) )
+    tm = timemesh(self.tstart, self.tend, self.nslices, impeuler, mat, self.nfine, self.ncoarse, 1e-10, 5)
+    tm.run_coarse(self.u0)
+    u = np.array([1.0])
+    Mat = tm.slices[0].int_coarse.get_update_matrix(self.u0)
+    b = self.u0.y
+    for i in range(0,3):
+      # matrix update to end of slice
+      b = Mat.dot(b)
+      err = np.linalg.norm( tm.get_coarse_value(i).y - b, np.inf )
+      assert err<2e-12, ("Successive application of update matrix does not reproduce intermediata coarse values generated by run_coarse. Error: %5.3e in slice %2i" % (err, i))
+
   # run_fine provides expected intermediate values
   def test_runfineintermediate(self):
     tm = timemesh(self.tstart, self.tend, self.nslices, impeuler, impeuler, self.nfine, self.ncoarse, 1e-10, 5)