Coverage for pySDC/projects/FastWaveSlowWave/plot_stifflimit_specrad.py: 95%
82 statements
« prev ^ index » next coverage.py v7.5.0, created at 2024-04-29 09:02 +0000
1import matplotlib
3matplotlib.use('Agg')
5import numpy as np
6from matplotlib import pyplot as plt
7from pylab import rcParams
9from pySDC.implementations.problem_classes.FastWaveSlowWave_0D import swfw_scalar
10from pySDC.implementations.sweeper_classes.imex_1st_order import imex_1st_order
13from pySDC.core.Step import step
16# noinspection PyShadowingNames
17def compute_specrad():
18 """
19 Routine to compute spectral radius and norm of the error propagation matrix E
21 Returns:
22 numpy.nparray: list of number of nodes
23 numpy.nparray: list of fast lambdas
24 numpy.nparray: list of spectral radii
25 numpy.nparray: list of norms
27 """
28 problem_params = dict()
29 # SET VALUE FOR lambda_slow AND VALUES FOR lambda_fast ###
30 problem_params['lambda_s'] = np.array([1.0 * 1j], dtype='complex')
31 problem_params['lambda_f'] = np.array([50.0 * 1j, 100.0 * 1j], dtype='complex')
32 problem_params['u0'] = 1.0
34 # initialize sweeper parameters
35 sweeper_params = dict()
36 # SET TYPE OF QUADRATURE NODES ###
37 sweeper_params['quad_type'] = 'RADAU-RIGHT'
39 # initialize level parameters
40 level_params = dict()
41 level_params['dt'] = 1.0
42 t0 = 0.0
44 # fill description dictionary for easy step instantiation
45 description = dict()
46 description['problem_class'] = swfw_scalar # pass problem class
47 description['problem_params'] = problem_params # pass problem parameters
48 description['sweeper_class'] = imex_1st_order # pass sweeper (see part B)
49 description['level_params'] = level_params # pass level parameters
50 description['step_params'] = dict() # pass step parameters
52 nodes_v = np.arange(2, 10)
53 specrad = np.zeros((3, np.size(nodes_v)))
54 norm = np.zeros((3, np.size(nodes_v)))
56 for i in range(0, np.size(nodes_v)):
57 sweeper_params['num_nodes'] = nodes_v[i]
58 description['sweeper_params'] = sweeper_params # pass sweeper parameters
60 # now the description contains more or less everything we need to create a step
61 S = step(description=description)
63 L = S.levels[0]
64 P = L.prob
66 u0 = S.levels[0].prob.u_exact(t0)
67 S.init_step(u0)
68 QE = L.sweep.QE[1:, 1:]
69 QI = L.sweep.QI[1:, 1:]
70 Q = L.sweep.coll.Qmat[1:, 1:]
71 nnodes = L.sweep.coll.num_nodes
72 dt = L.params.dt
74 assert nnodes == nodes_v[i], 'Something went wrong during instantiation, nnodes is not correct, got %s' % nnodes
76 for j in range(0, 2):
77 LHS = np.eye(nnodes) - dt * (P.lambda_f[j] * QI + P.lambda_s[0] * QE)
78 RHS = dt * ((P.lambda_f[j] + P.lambda_s[0]) * Q - (P.lambda_f[j] * QI + P.lambda_s[0] * QE))
79 evals, evecs = np.linalg.eig(np.linalg.inv(LHS).dot(RHS))
80 specrad[j + 1, i] = np.linalg.norm(evals, np.inf)
81 norm[j + 1, i] = np.linalg.norm(np.linalg.inv(LHS).dot(RHS), np.inf)
83 if L.sweep.coll.left_is_node:
84 # For Lobatto nodes, first column and row are all zeros, since q_1 = q_0; hence remove them
85 QI = QI[1:, 1:]
86 Q = Q[1:, 1:]
87 # Eigenvalue of error propagation matrix in stiff limit: E = I - inv(QI)*Q
88 evals, evecs = np.linalg.eig(np.eye(nnodes - 1) - np.linalg.inv(QI).dot(Q))
89 norm[0, i] = np.linalg.norm(np.eye(nnodes - 1) - np.linalg.inv(QI).dot(Q), np.inf)
90 else:
91 evals, evecs = np.linalg.eig(np.eye(nnodes) - np.linalg.inv(QI).dot(Q))
92 norm[0, i] = np.linalg.norm(np.eye(nnodes) - np.linalg.inv(QI).dot(Q), np.inf)
93 specrad[0, i] = np.linalg.norm(evals, np.inf)
95 print("Spectral radius of infinitely fast wave case > 1.0 for M=%2i" % nodes_v[np.argmax(specrad[0, :] > 1.0)])
96 print("Spectral radius of > 1.0 for M=%2i" % nodes_v[np.argmax(specrad[1, :] > 1.0)])
98 return nodes_v, problem_params['lambda_f'], specrad, norm
101# noinspection PyShadowingNames
102def plot_specrad(nodes_v, lambda_f, specrad, norm):
103 """
104 Plotting function for spectral radii and norms
106 Args:
107 nodes_v (numpy.nparray): list of number of nodes
108 lambda_f (numpy.nparray): list of fast lambdas
109 specrad (numpy.nparray): list of spectral radii
110 norm (numpy.nparray): list of norms
111 """
112 fs = 8
113 rcParams['figure.figsize'] = 2.5, 2.5
114 rcParams['pgf.rcfonts'] = False
115 fig = plt.figure()
116 plt.plot(nodes_v, specrad[0, :], 'rd-', markersize=fs - 2, label=r'$\lambda_{fast} = \infty$')
117 plt.plot(nodes_v, specrad[1, :], 'bo-', markersize=fs - 2, label=r'$\lambda_{fast} = %2.0f $' % lambda_f[0].imag)
118 plt.plot(nodes_v, specrad[2, :], 'gs-', markersize=fs - 2, label=r'$\lambda_{fast} = %2.0f $' % lambda_f[1].imag)
119 plt.xlabel(r'Number of nodes $M$', fontsize=fs)
120 plt.ylabel(r'Spectral radius $\sigma\left( \mathbf{E} \right)$', fontsize=fs, labelpad=2)
121 plt.legend(loc='lower right', fontsize=fs, prop={'size': fs})
122 plt.xlim([np.min(nodes_v), np.max(nodes_v)])
123 plt.ylim([0, 1.0])
124 plt.yticks(fontsize=fs)
125 plt.xticks(fontsize=fs)
126 filename = 'data/stifflimit-specrad.png'
127 fig.savefig(filename, bbox_inches='tight')
129 fig = plt.figure()
130 plt.plot(nodes_v, norm[0, :], 'rd-', markersize=fs - 2, label=r'$\lambda_{fast} = \infty$')
131 plt.plot(nodes_v, norm[1, :], 'bo-', markersize=fs - 2, label=r'$\lambda_{fast} = %2.0f $' % lambda_f[0].imag)
132 plt.plot(nodes_v, norm[2, :], 'gs-', markersize=fs - 2, label=r'$\lambda_{fast} = %2.0f $' % lambda_f[1].imag)
133 plt.xlabel(r'Number of nodes $M$', fontsize=fs)
134 plt.ylabel(r'Norm $\left|| \mathbf{E} \right||_{\infty}$', fontsize=fs, labelpad=2)
135 plt.legend(loc='lower right', fontsize=fs, prop={'size': fs})
136 plt.xlim([np.min(nodes_v), np.max(nodes_v)])
137 plt.ylim([0, 2.4])
138 plt.yticks(fontsize=fs)
139 plt.xticks(fontsize=fs)
140 filename = 'data/stifflimit-norm.png'
141 fig.savefig(filename, bbox_inches='tight')
144if __name__ == "__main__":
145 nodes_v, lambda_f, specrad, norm = compute_specrad()
146 plot_specrad(nodes_v, lambda_f, specrad, norm)