Coverage for pySDC/projects/FastWaveSlowWave/runitererror_acoustic.py: 100%

1import matplotlib 

2 

3matplotlib.use('Agg') 

4 

5import numpy as np 

6from matplotlib import pyplot as plt 

7from pylab import rcParams 

8 

9from pySDC.projects.FastWaveSlowWave.HookClass_acoustic import dump_energy 

10 

11from pySDC.implementations.problem_classes.AcousticAdvection_1D_FD_imex import acoustic_1d_imex 

12from pySDC.implementations.sweeper_classes.imex_1st_order import imex_1st_order 

13from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI 

14 

15from pySDC.helpers.stats_helper import filter_stats 

16 

17 

18def compute_and_plot_itererror(): 

19 """ 

20 Routine to compute and plot the error over the iterations for difference cs values 

21 """ 

22 

23 num_procs = 1 

24 

25 t0 = 0.0 

26 Tend = 0.025 

27 

28 # initialize level parameters 

29 level_params = dict() 

30 level_params['restol'] = 1e-10 

31 level_params['dt'] = Tend 

32 

33 # This comes as read-in for the step class 

34 step_params = dict() 

35 step_params['maxiter'] = 15 

36 

37 # This comes as read-in for the problem class 

38 problem_params = dict() 

39 problem_params['cadv'] = 0.1 

40 problem_params['nvars'] = [(2, 300)] 

41 problem_params['order_adv'] = 5 

42 problem_params['waveno'] = 5 

43 

44 # This comes as read-in for the sweeper class 

45 sweeper_params = dict() 

46 sweeper_params['quad_type'] = 'RADAU-RIGHT' 

47 sweeper_params['do_coll_update'] = True 

48 

49 # initialize controller parameters 

50 controller_params = dict() 

51 controller_params['logger_level'] = 30 

52 controller_params['hook_class'] = dump_energy 

53 

54 # Fill description dictionary for easy hierarchy creation 

55 description = dict() 

56 description['problem_class'] = acoustic_1d_imex 

57 description['sweeper_class'] = imex_1st_order 

58 description['step_params'] = step_params 

59 description['level_params'] = level_params 

60 

61 cs_v = [0.5, 1.0, 1.5, 5.0] 

62 nodes_v = [3] 

63 

64 residual = np.zeros((np.size(cs_v), np.size(nodes_v), step_params['maxiter'])) 

65 convrate = np.zeros((np.size(cs_v), np.size(nodes_v), step_params['maxiter'] - 1)) 

66 lastiter = np.zeros((np.size(cs_v), np.size(nodes_v))) + step_params['maxiter'] 

67 avg_convrate = np.zeros((np.size(cs_v), np.size(nodes_v))) 

68 

69 P = None 

70 for cs_ind in range(0, np.size(cs_v)): 

71 problem_params['cs'] = cs_v[cs_ind] 

72 description['problem_params'] = problem_params 

73 

74 for nodes_ind in np.arange(np.size(nodes_v)): 

75 sweeper_params['num_nodes'] = nodes_v[nodes_ind] 

76 description['sweeper_params'] = sweeper_params 

77 

78 # instantiate the controller 

79 controller = controller_nonMPI( 

80 num_procs=num_procs, controller_params=controller_params, description=description 

81 ) 

82 # get initial values on finest level 

83 P = controller.MS[0].levels[0].prob 

84 uinit = P.u_exact(t0) 

85 

86 print("Fast CFL number: %4.2f" % (problem_params['cs'] * level_params['dt'] / P.dx)) 

87 print("Slow CFL number: %4.2f" % (problem_params['cadv'] * level_params['dt'] / P.dx)) 

88 

89 # call main function to get things done... 

90 uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend) 

91 

92 extract_stats = filter_stats(stats, type='residual_post_iteration') 

93 

94 for k, v in extract_stats.items(): 

95 if k.iter != -1: 

96 residual[cs_ind, nodes_ind, k.iter - 2] = v 

97 

98 # Compute convergence rates 

99 for iter in range(0, step_params['maxiter'] - 1): 

100 if residual[cs_ind, nodes_ind, iter] < level_params['restol']: 

101 lastiter[cs_ind, nodes_ind] = iter 

102 else: 

103 convrate[cs_ind, nodes_ind, iter] = ( 

104 residual[cs_ind, nodes_ind, iter + 1] / residual[cs_ind, nodes_ind, iter] 

105 ) 

106 print(lastiter[cs_ind, nodes_ind]) 

107 avg_convrate[cs_ind, nodes_ind] = np.sum(convrate[cs_ind, nodes_ind, :]) / float( 

108 lastiter[cs_ind, nodes_ind] 

109 ) 

110 

111 # Plot the results 

112 fs = 8 

113 color = ['r', 'b', 'g', 'c'] 

114 shape = ['o', 'd', 's', 'v'] 

115 rcParams['figure.figsize'] = 2.5, 2.5 

116 rcParams['pgf.rcfonts'] = False 

117 fig = plt.figure() 

118 for ii in range(0, np.size(cs_v)): 

119 x = np.arange(1, lastiter[ii, 0] - 1) 

120 y = convrate[ii, 0, 0 : int(lastiter[ii, 0]) - 2] 

121 plt.plot( 

122 x, 

123 y, 

124 linestyle='-', 

125 marker=shape[ii], 

126 markersize=fs - 2, 

127 color=color[ii], 

128 label=r'$C_{fast}$=%4.2f' % (cs_v[ii] * level_params['dt'] / P.dx), 

129 ) 

130 

131 plt.legend(loc='upper right', fontsize=fs, prop={'size': fs - 2}) 

132 plt.xlabel('Iteration', fontsize=fs) 

133 plt.ylabel(r'$|| r^{k+1} ||_{\infty}/|| r^k ||_{\infty}$', fontsize=fs, labelpad=2) 

134 plt.xlim([0, step_params['maxiter']]) 

135 plt.ylim([0, 1.0]) 

136 plt.yticks(fontsize=fs) 

137 plt.xticks(fontsize=fs) 

138 filename = 'data/iteration.png' 

139 fig.savefig(filename, bbox_inches='tight') 

140 

141 

142if __name__ == "__main__": 

143 compute_and_plot_itererror()