Coverage for pySDC/projects/FastWaveSlowWave/plot_stability.py: 95%

1import matplotlib 

2 

3matplotlib.use('Agg') 

4 

5import numpy as np 

6from matplotlib import pyplot as plt 

7from pylab import rcParams 

8from matplotlib.patches import Polygon 

9 

10from pySDC.implementations.problem_classes.FastWaveSlowWave_0D import swfw_scalar 

11from pySDC.implementations.sweeper_classes.imex_1st_order import imex_1st_order 

12 

13 

14from pySDC.core.Step import step 

15 

16 

17# noinspection PyShadowingNames 

18def compute_stability(): 

19 """ 

20 Routine to compute the stability domains of different configurations of fwsw-SDC 

21 

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 

31 

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) 

36 

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 

42 

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 

49 

50 # initialize level parameters 

51 level_params = dict() 

52 level_params['dt'] = 1.0 

53 

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 

62 

63 # SET NUMBER OF ITERATIONS - SET K=0 FOR COLLOCATION SOLUTION ### 

64 K = 3 

65 

66 # now the description contains more or less everything we need to create a step 

67 S = step(description=description) 

68 

69 L = S.levels[0] 

70 

71 Q = L.sweep.coll.Qmat[1:, 1:] 

72 nnodes = L.sweep.coll.num_nodes 

73 dt = L.params.dt 

74 

75 stab = np.zeros((N_f, N_s), dtype='complex') 

76 

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 

95 

96 return lambda_s, lambda_f, sweeper_params['num_nodes'], K, stab 

97 

98 

99# noinspection PyShadowingNames 

100def plot_stability(lambda_s, lambda_f, num_nodes, K, stab): 

101 """ 

102 Plotting routine of the stability domains 

103 

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 """ 

111 

112 lam_s_max = np.amax(lambda_s.imag) 

113 lam_f_max = np.amax(lambda_f.imag) 

114 

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') 

149 

150 

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)