Coverage for pySDC/projects/AsympConv/smoother_specrad_heatmap.py: 0%

1import matplotlib 

2import numpy as np 

3import scipy.linalg as LA 

4 

5matplotlib.use('Agg') 

6import matplotlib.pyplot as plt 

7from matplotlib.colors import LogNorm 

8 

9from pySDC.core.Collocation import CollBase 

10 

11 

12def compute_and_plot_specrad(): 

13 """ 

14 Compute and plot spectral radius of smoother iteration matrix for a whole range of eigenvalues 

15 """ 

16 # setup_list = [('LU', 'to0'), ('LU', 'toinf'), ('IE', 'to0'), ('IE', 'toinf')] 

17 # setup_list = [('LU', 'to0'), ('LU', 'toinf')] 

18 # setup_list = [('IE', 'to0'), ('IE', 'toinf')] 

19 # setup_list = [('LU', 'toinf'), ('IE', 'toinf')] 

20 setup_list = [('IE', 'full'), ('LU', 'full')] 

21 # setup_list = [('EX', 'to0'), ('PIC', 'to0')] 

22 

23 # set up plotting parameters 

24 params = { 

25 'legend.fontsize': 24, 

26 'figure.figsize': (12, 8), 

27 'axes.labelsize': 24, 

28 'axes.titlesize': 24, 

29 'xtick.labelsize': 24, 

30 'ytick.labelsize': 24, 

31 'lines.linewidth': 3, 

32 } 

33 plt.rcParams.update(params) 

34 

35 Nnodes = 3 

36 Nsteps = 4 

37 

38 coll = CollBase(Nnodes, 0, 1, quad_type='RADAU-RIGHT') 

39 Qmat = coll.Qmat[1:, 1:] 

40 

41 Nmat = np.zeros((Nnodes, Nnodes)) 

42 Nmat[:, -1] = 1 

43 

44 Emat = np.zeros((Nsteps, Nsteps)) 

45 np.fill_diagonal(Emat[1:, :], 1) 

46 

47 for qd_type, conv_type in setup_list: 

48 if qd_type == 'LU': 

49 QT = coll.Qmat[1:, 1:].T 

50 [_, _, U] = LA.lu(QT, overwrite_a=True) 

51 QDmat = U.T 

52 

53 elif qd_type == 'IE': 

54 QI = np.zeros(np.shape(coll.Qmat)) 

55 for m in range(coll.num_nodes + 1): 

56 QI[m, 1 : m + 1] = coll.delta_m[0:m] 

57 QDmat = QI[1:, 1:] 

58 

59 elif qd_type == 'EE': 

60 QE = np.zeros(np.shape(coll.Qmat)) 

61 for m in range(coll.num_nodes + 1): 

62 QE[m, 0:m] = coll.delta_m[0:m] 

63 QDmat = QE[1:, 1:] 

64 

65 elif qd_type == 'PIC': 

66 QDmat = np.zeros(np.shape(coll.Qmat[1:, 1:])) 

67 

68 elif qd_type == 'EX': 

69 QT = coll.Qmat[1:, 1:].T 

70 [_, _, U] = LA.lu(QT, overwrite_a=True) 

71 QDmat = np.tril(U.T, k=-1) 

72 print(QDmat) 

73 

74 else: 

75 raise NotImplementedError('qd_type %s is not implemented' % qd_type) 

76 

77 # lim_specrad = max(abs(np.linalg.eigvals(np.eye(Nnodes) - np.linalg.inv(QDmat).dot(Qmat)))) 

78 # print('qd_type: %s -- lim_specrad: %6.4e -- conv_type: %s' % (qd_type, lim_specrad, conv_type)) 

79 

80 if conv_type == 'to0': 

81 ilim_left = -4 

82 ilim_right = 2 

83 rlim_left = 2 

84 rlim_right = -4 

85 

86 elif conv_type == 'toinf': 

87 ilim_left = 0 

88 ilim_right = 11 

89 rlim_left = 6 

90 rlim_right = 0 

91 

92 elif conv_type == 'full': 

93 ilim_left = -10 

94 ilim_right = 11 

95 rlim_left = 10 

96 rlim_right = -11 

97 

98 else: 

99 raise NotImplementedError('conv_type %s is not implemented' % conv_type) 

100 

101 ilam_list = 1j * np.logspace(ilim_left, ilim_right, 201) 

102 rlam_list = -1 * np.logspace(rlim_left, rlim_right, 201) 

103 

104 assert (rlim_right - rlim_left + 1) % 5 == 0 

105 assert (ilim_right - ilim_left - 1) % 5 == 0 

106 assert (len(rlam_list) - 1) % 5 == 0 

107 assert (len(ilam_list) - 1) % 5 == 0 

108 

109 Prho = np.zeros((len(rlam_list), len(ilam_list))) 

110 

111 for idr, rlam in enumerate(rlam_list): 

112 for idi, ilam in enumerate(ilam_list): 

113 dxlam = rlam + ilam 

114 

115 mat = np.linalg.inv(np.eye(Nnodes * Nsteps) - dxlam * np.kron(np.eye(Nsteps), QDmat)).dot( 

116 dxlam * np.kron(np.eye(Nsteps), (Qmat - QDmat)) + np.kron(Emat, Nmat) 

117 ) 

118 mat = np.linalg.matrix_power(mat, Nnodes) 

119 

120 Prho[idr, idi] = max(abs(np.linalg.eigvals(mat))) 

121 

122 print(np.amax(Prho)) 

123 

124 fig, ax = plt.subplots(figsize=(15, 10)) 

125 

126 ax.set_xticks([i + 0.5 for i in range(0, len(rlam_list), int(len(rlam_list) / 5))]) 

127 ax.set_xticklabels( 

128 [r'-$10^{%d}$' % i for i in range(rlim_left, rlim_right, int((rlim_right - rlim_left + 1) / 5))] 

129 ) 

130 ax.set_yticks([i + 0.5 for i in range(0, len(ilam_list), int(len(ilam_list) / 5))]) 

131 ax.set_yticklabels( 

132 [r'$10^{%d}i$' % i for i in range(ilim_left, ilim_right, int((ilim_right - ilim_left - 1) / 5))] 

133 ) 

134 

135 cmap = plt.get_cmap('Reds') 

136 pcol = plt.pcolor(Prho.T, cmap=cmap, norm=LogNorm(vmin=1e-09, vmax=1e-00)) 

137 

138 plt.colorbar(pcol) 

139 

140 plt.xlabel(r'$Re(\Delta t\lambda)$') 

141 plt.ylabel(r'$Im(\Delta t\lambda)$') 

142 

143 fname = ( 

144 'data/heatmap_smoother_' + conv_type + '_Nsteps' + str(Nsteps) + '_M' + str(Nnodes) + '_' + qd_type + '.png' 

145 ) 

146 plt.savefig(fname, transparent=True, bbox_inches='tight') 

147 

148 

149if __name__ == "__main__": 

150 compute_and_plot_specrad()