Coverage for pySDC/projects/FastWaveSlowWave/plot_stability.py: 95%
75 statements
« prev ^ index » next coverage.py v7.6.9, created at 2024-12-20 14:51 +0000
1import matplotlib
3matplotlib.use('Agg')
5import numpy as np
6from matplotlib import pyplot as plt
7from pylab import rcParams
8from matplotlib.patches import Polygon
10from pySDC.implementations.problem_classes.FastWaveSlowWave_0D import swfw_scalar
11from pySDC.implementations.sweeper_classes.imex_1st_order import imex_1st_order
14from pySDC.core.step import Step
17# noinspection PyShadowingNames
18def compute_stability():
19 """
20 Routine to compute the stability domains of different configurations of fwsw-SDC
22 Returns:
23 numpy.ndarray: lambda_slow
24 numpy.ndarray: lambda_fast
25 int: number of collocation nodes
26 int: number of iterations
27 numpy.ndarray: stability numbers
28 """
29 N_s = 100
30 N_f = 400
32 lam_s_max = 5.0
33 lam_f_max = 12.0
34 lambda_s = 1j * np.linspace(0.0, lam_s_max, N_s)
35 lambda_f = 1j * np.linspace(0.0, lam_f_max, N_f)
37 problem_params = dict()
38 # SET VALUE FOR lambda_slow AND VALUES FOR lambda_fast ###
39 problem_params['lambda_s'] = np.array([0.0])
40 problem_params['lambda_f'] = np.array([0.0])
41 problem_params['u0'] = 1.0
43 # initialize sweeper parameters
44 sweeper_params = dict()
45 # SET TYPE AND NUMBER OF QUADRATURE NODES ###
46 sweeper_params['quad_type'] = 'RADAU-RIGHT'
47 sweeper_params['num_nodes'] = 3
48 sweeper_params['do_coll_update'] = True
50 # initialize level parameters
51 level_params = dict()
52 level_params['dt'] = 1.0
54 # fill description dictionary for easy step instantiation
55 description = dict()
56 description['problem_class'] = swfw_scalar # pass problem class
57 description['problem_params'] = problem_params # pass problem parameters
58 description['sweeper_class'] = imex_1st_order # pass sweeper (see part B)
59 description['sweeper_params'] = sweeper_params # pass sweeper parameters
60 description['level_params'] = level_params # pass level parameters
61 description['step_params'] = dict() # pass step parameters
63 # SET NUMBER OF ITERATIONS - SET K=0 FOR COLLOCATION SOLUTION ###
64 K = 3
66 # now the description contains more or less everything we need to create a step
67 S = Step(description=description)
69 L = S.levels[0]
71 Q = L.sweep.coll.Qmat[1:, 1:]
72 nnodes = L.sweep.coll.num_nodes
73 dt = L.params.dt
75 stab = np.zeros((N_f, N_s), dtype='complex')
77 for i in range(0, N_s):
78 for j in range(0, N_f):
79 lambda_fast = lambda_f[j]
80 lambda_slow = lambda_s[i]
81 if K != 0:
82 lambdas = [lambda_fast, lambda_slow]
83 # LHS, RHS = L.sweep.get_scalar_problems_sweeper_mats(lambdas=lambdas)
84 Mat_sweep = L.sweep.get_scalar_problems_manysweep_mat(nsweeps=K, lambdas=lambdas)
85 else:
86 # Compute stability function of collocation solution
87 Mat_sweep = np.linalg.inv(np.eye(nnodes) - dt * (lambda_fast + lambda_slow) * Q)
88 if L.sweep.params.do_coll_update:
89 stab_fh = 1.0 + (lambda_fast + lambda_slow) * L.sweep.coll.weights.dot(Mat_sweep.dot(np.ones(nnodes)))
90 else:
91 q = np.zeros(nnodes)
92 q[nnodes - 1] = 1.0
93 stab_fh = q.dot(Mat_sweep.dot(np.ones(nnodes)))
94 stab[j, i] = stab_fh
96 return lambda_s, lambda_f, sweeper_params['num_nodes'], K, stab
99# noinspection PyShadowingNames
100def plot_stability(lambda_s, lambda_f, num_nodes, K, stab):
101 """
102 Plotting routine of the stability domains
104 Args:
105 lambda_s (numpy.ndarray): lambda_slow
106 lambda_f (numpy.ndarray): lambda_fast
107 num_nodes (int): number of collocation nodes
108 K (int): number of iterations
109 stab (numpy.ndarray): stability numbers
110 """
112 lam_s_max = np.amax(lambda_s.imag)
113 lam_f_max = np.amax(lambda_f.imag)
115 rcParams['figure.figsize'] = 1.5, 1.5
116 fs = 8
117 fig = plt.figure()
118 levels = np.array([0.25, 0.5, 0.75, 0.9, 1.1])
119 CS1 = plt.contour(lambda_s.imag, lambda_f.imag, np.absolute(stab), levels, colors='k', linestyles='dashed')
120 CS2 = plt.contour(lambda_s.imag, lambda_f.imag, np.absolute(stab), [1.0], colors='k')
121 # Set markers at points used in plot_stab_vs_k
122 plt.plot(4, 10, 'x', color='k', markersize=fs - 4)
123 plt.plot(1, 10, 'x', color='k', markersize=fs - 4)
124 plt.clabel(CS1, inline=True, fmt='%3.2f', fontsize=fs - 2)
125 manual_locations = [(1.5, 2.5)]
126 if K > 0: # for K=0 and no 1.0 isoline, this crashes Matplotlib for somer reason
127 plt.clabel(CS2, inline=True, fmt='%3.2f', fontsize=fs - 2, manual=manual_locations)
128 plt.gca().add_patch(
129 Polygon(
130 [[0, 0], [lam_s_max, 0], [lam_s_max, lam_s_max]],
131 visible=True,
132 fill=True,
133 facecolor='.75',
134 edgecolor='k',
135 linewidth=1.0,
136 zorder=11,
137 )
138 )
139 plt.gca().set_xticks(np.arange(0, int(lam_s_max) + 1))
140 plt.gca().set_yticks(np.arange(0, int(lam_f_max) + 2, 2))
141 plt.gca().tick_params(axis='both', which='both', labelsize=fs)
142 plt.xlim([0.0, lam_s_max])
143 plt.ylim([0.0, lam_f_max])
144 plt.xlabel(r'$\Delta t \lambda_{slow}$', fontsize=fs, labelpad=0.0)
145 plt.ylabel(r'$\Delta t \lambda_{fast}$', fontsize=fs, labelpad=0.0)
146 plt.title(r'$M=%1i$, $K=%1i$' % (num_nodes, K), fontsize=fs)
147 filename = 'data/stability-K' + str(K) + '-M' + str(num_nodes) + '.png'
148 fig.savefig(filename, bbox_inches='tight')
151if __name__ == "__main__":
152 lambda_s, lambda_f, num_nodes, K, stab = compute_stability()
153 plot_stability(lambda_s, lambda_f, num_nodes, K, stab)