|
| 1 | +import sys |
| 2 | +sys.path.append('../../src') |
| 3 | +import numpy as np |
| 4 | +from numpy import linalg as LA |
| 5 | +import scipy.sparse as sp |
| 6 | +import scipy.linalg as spla |
| 7 | +import scipy.sparse.linalg as linalg |
| 8 | + |
| 9 | +from parareal import parareal |
| 10 | +from impeuler import impeuler |
| 11 | +from intexact import intexact |
| 12 | +from integrator_dedalus import integrator_dedalus |
| 13 | +from trapezoidal import trapezoidal |
| 14 | +from solution_linear import solution_linear |
| 15 | +from solution_dedalus import solution_dedalus |
| 16 | +from get_matrix import get_upwind, get_centered, get_diffusion |
| 17 | +from parameter import parameter |
| 18 | + |
| 19 | +from pylab import rcParams |
| 20 | +import matplotlib.pyplot as plt |
| 21 | +from subprocess import call |
| 22 | + |
| 23 | +def ie(z): |
| 24 | + return 1.0/(1.0 - z) |
| 25 | + |
| 26 | +def trap(z): |
| 27 | + return (1.0 + 0.5*z)/(1.0 - 0.5*z) |
| 28 | + |
| 29 | +try: |
| 30 | + figure = int(sys.argv[1]) # 1 generates figure_1, 2 generates figure_2 |
| 31 | +except: |
| 32 | + print("No or wrong command line argument provided, creating figure 5. Use 5, 6, 7 or 8 as command line argument.") |
| 33 | + figure = 5 |
| 34 | +assert 5<= figure <= 8, "Figure should be 5, 6, 7 or 8" |
| 35 | + |
| 36 | +if figure==5 or figure==6: |
| 37 | + par = parameter(dedalus = False) |
| 38 | +elif figure==7 or figure==8: |
| 39 | + par = parameter(dedalus = True) |
| 40 | +else: |
| 41 | + sys.exit("This should have been caught above") |
| 42 | + |
| 43 | +Tend, nslices, maxiter, nfine, ncoarse, tol, epsilon, ndof_f = par.getpar() |
| 44 | + |
| 45 | +nsteps = [1, 2, 4, 8, 12, 16, 20] |
| 46 | +nsteps = [1, 10, 50, 100, 1000, 10000] |
| 47 | + |
| 48 | +ndof_c_v = [16, 24, 30] |
| 49 | +xaxis_f = np.linspace(0.0, 1.0, ndof_f+1)[0:ndof_f] |
| 50 | +dx_f = xaxis_f[1] - xaxis_f[0] |
| 51 | + |
| 52 | +if figure==5: |
| 53 | + A_f = get_upwind(ndof_f, dx_f) |
| 54 | + u0fine = solution_linear(np.zeros(ndof_f), A_f) |
| 55 | +elif figure==6 or figure==7: |
| 56 | + A_f = get_centered(ndof_f, dx_f) |
| 57 | + u0fine = solution_linear(np.zeros(ndof_f), A_f) |
| 58 | +elif figure==8: |
| 59 | + u0fine = solution_dedalus(np.zeros(ndof_f), ndof_f) |
| 60 | +else: |
| 61 | + sys.exit("This should have been caught above") |
| 62 | + |
| 63 | +norm_l2 = np.zeros((3,np.size(nsteps))) |
| 64 | +norm_inf = np.zeros((3,np.size(nsteps))) |
| 65 | +dt_f_v = np.zeros((1,np.size(nsteps))) |
| 66 | +dt_c_v = np.zeros((1,np.size(nsteps))) |
| 67 | +A_f_eig = np.zeros((3,np.size(nsteps)), dtype='complex') |
| 68 | + |
| 69 | +for nn in range(3): |
| 70 | + |
| 71 | + ndof_c = ndof_c_v[nn] |
| 72 | + xaxis_c = np.linspace(0.0, 1.0, ndof_c+1)[0:ndof_c] |
| 73 | + dx_c = xaxis_c[1] - xaxis_c[0] |
| 74 | + |
| 75 | + ### Eigenvalues of A_f |
| 76 | + eig_val, eig_vec = LA.eig(A_f.todense()) |
| 77 | + |
| 78 | + if figure==5: |
| 79 | + A_c = get_upwind(ndof_c, dx_c) |
| 80 | + u0coarse = solution_linear(np.zeros(ndof_c), A_c) |
| 81 | + filename = 'figure_5.pdf' |
| 82 | + elif figure==6: |
| 83 | + A_c = get_centered(ndof_c, dx_c) |
| 84 | + u0coarse = solution_linear(np.zeros(ndof_c), A_c) |
| 85 | + filename = 'figure_6.pdf' |
| 86 | + elif figure==7: |
| 87 | + A_c = get_centered(ndof_c, dx_c) |
| 88 | + u0coarse = solution_linear(np.zeros(ndof_c), A_c) |
| 89 | + filename = 'figure_7.pdf' |
| 90 | + elif figure==8: |
| 91 | + u0coarse = solution_dedalus(np.zeros(ndof_c), ndof_c) |
| 92 | + filename = 'figure_8.pdf' |
| 93 | + else: |
| 94 | + sys.exit("Value of figure should be") |
| 95 | + |
| 96 | + for mm in range(np.size(nsteps)): |
| 97 | + |
| 98 | + if figure==5 or figure==7: |
| 99 | + Rz = np.vectorize(ie) |
| 100 | + para = parareal(0.0, Tend, nslices, impeuler, impeuler, nsteps[mm], nsteps[mm], tol, maxiter, u0fine, u0coarse) |
| 101 | + elif figure==6: |
| 102 | + Rz = np.vectorize(trap) |
| 103 | + para = parareal(0.0, Tend, nslices, trapezoidal, trapezoidal, nsteps[mm], nsteps[mm], tol, maxiter, u0fine, u0coarse) |
| 104 | + elif figure==8: |
| 105 | + para = parareal(0.0, Tend, nslices, integrator_dedalus, integrator_dedalus, nsteps[mm], nsteps[mm], tol, maxiter, u0fine, u0coarse) |
| 106 | + else: |
| 107 | + sys.exit("Value of figure should be") |
| 108 | + Pmat, Bmat = para.get_parareal_matrix() |
| 109 | + dt_f_v[0,mm] = para.timemesh.slices[0].int_fine.dt |
| 110 | + dt_c_v[0,mm] = para.timemesh.slices[0].int_coarse.dt |
| 111 | + |
| 112 | + # Sort according to absolute values of R(z) with z = lambda*dt_f |
| 113 | + sort_index = np.argsort(np.abs(Rz(eig_val*dt_f_v[0,mm]))) |
| 114 | + |
| 115 | + # Store eigenvalue ndof_c+1. first "truncated" EV |
| 116 | + A_f_eig[nn,mm] = np.flip((eig_val[sort_index]))[ndof_c+1] |
| 117 | + |
| 118 | + ### Parareal iteration: y^k+1 = Pmat*y^k + Bmat*b |
| 119 | + norm_l2[nn,mm] = np.linalg.norm(Pmat.todense(), 2) |
| 120 | + |
| 121 | + |
| 122 | +rcParams['figure.figsize'] = 2.5, 2.5 |
| 123 | +fs = 8 |
| 124 | +ms = 4 |
| 125 | +fig = plt.figure(1) |
| 126 | +#plt.plot(dt_f_v[0,:], norm_l2[0,:], 'bo-', label='m='+str(ndof_c_v[0]), markersize=ms) |
| 127 | +#plt.plot(dt_f_v[0,:], norm_l2[1,:], 'rx-', label='m='+str(ndof_c_v[1]), markersize=ms) |
| 128 | +#plt.plot(dt_f_v[0,:], norm_l2[2,:], 'cd-', label='m='+str(ndof_c_v[2]), markersize=ms) |
| 129 | +#plt.loglog(dt_f_v[0,:], norm_l2[0,:]-np.abs(np.exp(np.multiply(A_f_eig[0,:],dt_f_v[0,:]))), 'bo-', label='m='+str(ndof_c_v[0]), markersize=ms) |
| 130 | +#plt.loglog(dt_f_v[0,:], norm_l2[1,:]-np.abs(np.exp(np.multiply(A_f_eig[1,:],dt_f_v[0,:]))), 'rx-', label='m='+str(ndof_c_v[1]), markersize=ms) |
| 131 | +print(norm_l2[1,:]) |
| 132 | +print(A_f_eig[1,:]) |
| 133 | +print("\n") |
| 134 | +print(norm_l2[2,:]) |
| 135 | +print(A_f_eig[2,:]) |
| 136 | + |
| 137 | +plt.loglog(dt_f_v[0,:], norm_l2[2,:]-0.0*np.abs(np.exp(np.multiply(A_f_eig[2,:],dt_f_v[0,:]))), 'cd-', label='m='+str(ndof_c_v[2]), markersize=ms) |
| 138 | +#plt.plot(dt_f_v[0,:], 1.0 + 0.0*dt_f_v[0,:], 'k:') |
| 139 | +plt.legend(loc='best', bbox_to_anchor=(0.5, 0.5), fontsize=fs, prop={'size':fs-2}, handlelength=3) |
| 140 | +plt.xlabel(r'$\delta t = \Delta t$', fontsize=fs) |
| 141 | +#plt.ylabel(r'$|| \mathbf{E} ||_2$', fontsize=fs) |
| 142 | +plt.ylabel(r'$|| \mathbf{E} ||_2 - | \exp(\lambda_{m+1} \delta t) |$', fontsize=fs) |
| 143 | + |
| 144 | +plt.xlim([0.0, dt_f_v[0,0]]) |
| 145 | +plt.ylim([1e-2, 1e1]) |
| 146 | +#plt.xlabel([0, maxiter]) |
| 147 | +plt.gcf().savefig(filename, bbox_inches='tight') |
| 148 | +call(["pdfcrop", filename, filename]) |
| 149 | +plt.show() |
0 commit comments