|
| 1 | +import sys |
| 2 | +sys.path.append('../src') |
| 3 | + |
| 4 | +from parareal import parareal |
| 5 | +from impeuler import impeuler |
| 6 | +from intexact import intexact |
| 7 | +from trapezoidal import trapezoidal |
| 8 | +from special_integrator import special_integrator |
| 9 | +from solution_linear import solution_linear |
| 10 | +import numpy as np |
| 11 | +import scipy.sparse as sparse |
| 12 | +import math |
| 13 | + |
| 14 | +from pylab import rcParams |
| 15 | +import matplotlib.pyplot as plt |
| 16 | +from matplotlib.patches import Polygon |
| 17 | +from subprocess import call |
| 18 | +import sympy |
| 19 | +from pylab import rcParams |
| 20 | + |
| 21 | +if __name__ == "__main__": |
| 22 | + |
| 23 | + Tend = 16.0 |
| 24 | + nslices = int(Tend) # Make sure each time slice has length 1 |
| 25 | + U_speed = 1.0 |
| 26 | + nu_v = [0.0, 1e-1, 0.5] |
| 27 | + ncoarse = 1 |
| 28 | + nfine = 10 |
| 29 | + dx = 1.0 |
| 30 | + Nsamples = 80 |
| 31 | + u0_val = np.array([[1.0]], dtype='complex') |
| 32 | + |
| 33 | + k_vec = np.linspace(0.0, np.pi, Nsamples+1, endpoint=False) |
| 34 | + |
| 35 | + svds = np.zeros((np.size(nu_v),np.size(k_vec))) |
| 36 | + |
| 37 | + for j in range(0,np.size(nu_v)): |
| 38 | + nu = nu_v[j] |
| 39 | + for i in range(0,np.size(k_vec)): |
| 40 | + waveno = k_vec[i] |
| 41 | + symb = -(1j*U_speed*waveno + nu*waveno**2) |
| 42 | + symb_coarse = symb |
| 43 | +# symb_coarse = -(1.0/dx)*(1.0 - np.exp(-1j*waveno*dx)) |
| 44 | + |
| 45 | + # Solution objects define the problem |
| 46 | + u0 = solution_linear(u0_val, np.array([[symb]],dtype='complex')) |
| 47 | + ucoarse = solution_linear(u0_val, np.array([[symb_coarse]],dtype='complex')) |
| 48 | + |
| 49 | + para = parareal(0.0, Tend, nslices, intexact, impeuler, nfine, ncoarse, 0.0, 0, u0) |
| 50 | + svds[j,i] = para.get_max_svd(ucoarse=ucoarse) |
| 51 | + |
| 52 | + rcParams['figure.figsize'] = 2.5, 2.5 |
| 53 | + fs = 8 |
| 54 | + fig = plt.figure() |
| 55 | + plt.plot(k_vec, svds[0,:], label=(r'$\nu$=%3.2f' % nu_v[0])) |
| 56 | + plt.plot(k_vec, svds[1,:], label=(r'$\nu$=%3.2f' % nu_v[1])) |
| 57 | + plt.plot(k_vec, svds[2,:], label=(r'$\nu$=%3.2f' % nu_v[2])) |
| 58 | + plt.gca().tick_params(axis='both', which='major', labelsize=fs-2) |
| 59 | + plt.gca().tick_params(axis='x', which='minor', bottom='off') |
| 60 | + plt.xlabel(r'Wave number $\kappa$', fontsize=fs, labelpad=1) |
| 61 | + plt.xlim([k_vec[0], k_vec[-1]]) |
| 62 | + plt.ylabel(r'Maximum singular value $\sigma$', fontsize=fs, labelpad=0) |
| 63 | + plt.legend(loc='upper left', fontsize=fs, prop={'size':fs-2}) |
| 64 | + filename='svd_vs_waveno.pdf' |
| 65 | + fig.savefig(filename, bbox_inches='tight') |
| 66 | + call(["pdfcrop", filename, filename]) |
| 67 | +# plt.show() |
0 commit comments