From 9564bbb8234f34432a6a3243aa005945617c917e Mon Sep 17 00:00:00 2001 From: Daniel Ruprecht Date: Fri, 22 Apr 2016 11:09:56 +0100 Subject: [PATCH] adds a function to parareal to return the maximum singular value of the error propagatin matrix --- scripts/plot_dispersion.py | 20 +++++++++++++++++--- src/parareal.py | 6 ++++++ 2 files changed, 23 insertions(+), 3 deletions(-) diff --git a/scripts/plot_dispersion.py b/scripts/plot_dispersion.py index 9ab52fe..245d7d3 100644 --- a/scripts/plot_dispersion.py +++ b/scripts/plot_dispersion.py @@ -52,6 +52,7 @@ def normalise(R, T, target): phase = np.zeros((6,Nsamples)) amp_factor = np.zeros((6,Nsamples)) + svds = np.zeros((3,Nsamples)) u0_val = np.array([[1.0]], dtype='complex') targets = np.zeros((3,Nsamples)) @@ -89,13 +90,13 @@ def normalise(R, T, target): amp_factor[2,i] = np.exp(sol_coarse.imag) # Compute Parareal phase velocity and amplification factor - + 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) if i==0: targets[jj,0] = np.angle(stab_ex) - + stab_para_norm = normalise(stab_para[0,0], Tend, targets[jj,i]) # Make sure that stab_norm*dt = stab err = abs(stab_para_norm**Tend - stab_para) @@ -112,7 +113,6 @@ def normalise(R, T, target): phase[3+jj,i] = sol_para.real/k_vec[i] amp_factor[3+jj,i] = np.exp(sol_para.imag) - ### rcParams['figure.figsize'] = 2.5, 2.5 fs = 8 @@ -155,3 +155,17 @@ def normalise(R, T, target): plt.gcf().savefig(filename, bbox_inches='tight') call(["pdfcrop", filename, filename]) + fig = plt.figure() + plt.plot(k_vec, svds[0,:], '-s', color='r', linewidth=1.5, markevery=(1,6), mew=1.0, markersize=fs/2) + plt.xlabel('Wave number', fontsize=fs, labelpad=0.25) + plt.ylabel('Maximal singular value', fontsize=fs, labelpad=0.5) + fig.gca().tick_params(axis='both', labelsize=fs) + plt.xlim([k_vec[0], k_vec[-1:]]) + plt.ylim([0, 2.0]) + # plt.legend(loc='lower left', fontsize=fs, prop={'size':fs-2}) + plt.gca().set_ylim([0.0, 2.0]) + plt.xticks([0, 1, 2, 3], fontsize=fs) + #plt.show() + filename = 'parareal-dispersion-svd.pdf' + plt.gcf().savefig(filename, bbox_inches='tight') + call(["pdfcrop", filename, filename]) diff --git a/src/parareal.py b/src/parareal.py index 88f3df9..f491f3b 100644 --- a/src/parareal.py +++ b/src/parareal.py @@ -96,6 +96,12 @@ def get_parareal_stab_function(self, k, ucoarse=None): Mat[:,i] = R.dot(M).flatten() return Mat + # Returns the largest singular value of the error propagation matrix + def get_max_svd(self, ucoarse=None): + Pmat, Bmat = self.get_parareal_matrix(ucoarse) + svds = linalg.svds(Pmat, k=1, return_singular_vectors=False) + return svds[0] + # Returns array containing all intermediate solutions def get_parareal_vector(self): b = np.zeros((self.u0.ndof*(self.timemesh.nslices+1),1))