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

1import matplotlib 

2import numpy as np 

3import scipy.linalg as LA 

4 

5matplotlib.use('Agg') 

6import matplotlib.pylab as plt 

7from matplotlib import rc 

8 

9from pySDC.core.collocation import CollBase 

10 

11 

12def compute_and_plot_specrad(Nnodes, lam): 

13 """ 

14 Compute and plot the spectral radius of the smoother for different step-sizes 

15 

16 Args: 

17 Nnodes: number of collocation nodes 

18 lam: test parameter representing the spatial problem 

19 """ 

20 

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

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

23 

24 # do LU decomposition of QT (St. Martin's trick) 

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

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

27 QDmat = U.T 

28 

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

30 Nmat[:, -1] = 1 

31 

32 Nsteps_list = [64, 256] 

33 color_list = ['red', 'blue'] 

34 marker_list = ['s', 'o'] 

35 

36 setup_list = zip(Nsteps_list, color_list, marker_list) 

37 

38 xlist = [0.1**i for i in range(11)] 

39 

40 rc('font', **{"sans-serif": ["Arial"], "size": 24}) 

41 plt.subplots(figsize=(15, 10)) 

42 

43 for Nsteps, color, marker in setup_list: 

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

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

46 

47 Prho_list = [] 

48 predict_list = [] 

49 for x in xlist: 

50 mat = np.linalg.inv(np.eye(Nnodes * Nsteps) - x * lam * np.kron(np.eye(Nsteps), QDmat)).dot( 

51 x * lam * np.kron(np.eye(Nsteps), (Qmat - QDmat)) + np.kron(Emat, Nmat) 

52 ) 

53 

54 Prho_list.append(max(abs(np.linalg.eigvals(mat)))) 

55 # predict_list.append((1 + x) ** (1.0 - 1.0 / (Nnodes * Nsteps)) * x ** (1.0 / (Nnodes * Nsteps))) 

56 predict_list.append(x ** (1.0 / (Nsteps))) 

57 

58 if len(predict_list) > 1: 

59 print( 

60 x, 

61 predict_list[-1], 

62 Prho_list[-1], 

63 Prho_list[-2] / Prho_list[-1], 

64 predict_list[-2] / predict_list[-1], 

65 ) 

66 

67 plt.loglog( 

68 xlist, 

69 Prho_list, 

70 linestyle='-', 

71 linewidth=3, 

72 color=color, 

73 marker=marker, 

74 markersize=10, 

75 label='spectral radius, L=' + str(Nsteps), 

76 ) 

77 plt.loglog( 

78 xlist, 

79 predict_list, 

80 linestyle='--', 

81 linewidth=2, 

82 color=color, 

83 marker=marker, 

84 markersize=10, 

85 label='estimate, L=' + str(Nsteps), 

86 ) 

87 

88 ax = plt.gca() 

89 ax.invert_xaxis() 

90 

91 plt.xlabel('time-step size') 

92 plt.ylabel('spectral radius') 

93 plt.legend(loc=3, numpoints=1) 

94 plt.grid() 

95 plt.ylim([1e-02, 1e01]) 

96 

97 if type(lam) is complex: 

98 fname = 'data/smoother_specrad_to0_L64+256_M' + str(Nnodes) + 'LU_imag.png' 

99 else: 

100 fname = 'data/smoother_specrad_to0_L64+256_M' + str(Nnodes) + 'LU_real.png' 

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

102 

103 

104if __name__ == "__main__": 

105 compute_and_plot_specrad(Nnodes=3, lam=-1) 

106 compute_and_plot_specrad(Nnodes=3, lam=1j)