From 9cad92c4489256e59728a3c4bd8c577ad8c66e9e Mon Sep 17 00:00:00 2001 From: Daniel Ruprecht Date: Mon, 25 Apr 2016 16:15:16 +0100 Subject: [PATCH] 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 --- scripts/plot_dispersion.py | 9 +++++++-- src/special_integrator.py | 4 ++-- src/timemesh.py | 6 +++++- tests/test_timemesh.py | 34 ++++++++++++++++++++++++++++++++++ 4 files changed, 48 insertions(+), 5 deletions(-) diff --git a/scripts/plot_dispersion.py b/scripts/plot_dispersion.py index 98e5ad1..af40410 100644 --- a/scripts/plot_dispersion.py +++ b/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) diff --git a/src/special_integrator.py b/src/special_integrator.py index f22243f..c00d09a 100644 --- a/src/special_integrator.py +++ b/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 diff --git a/src/timemesh.py b/src/timemesh.py index 561757b..e626e7b 100644 --- a/src/timemesh.py +++ b/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 diff --git a/tests/test_timemesh.py b/tests/test_timemesh.py index 0b8fa0d..d4133a8 100644 --- a/tests/test_timemesh.py +++ b/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)