Coverage for pySDC/projects/FastWaveSlowWave/runconvergence_acoustic.py: 50%

1import matplotlib 

2 

3matplotlib.use('Agg') 

4 

5import numpy as np 

6from matplotlib import pyplot as plt 

7from pylab import rcParams 

8from matplotlib.ticker import ScalarFormatter 

9 

10from pySDC.projects.FastWaveSlowWave.HookClass_acoustic import dump_energy 

11 

12from pySDC.implementations.problem_classes.AcousticAdvection_1D_FD_imex import acoustic_1d_imex 

13from pySDC.implementations.sweeper_classes.imex_1st_order import imex_1st_order 

14from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI 

15 

16 

17def compute_convergence_data(cwd=''): 

18 """ 

19 Routine to run the 1d acoustic-advection example with different orders 

20 

21 Args: 

22 cwd (string): current working directory 

23 """ 

24 

25 num_procs = 1 

26 

27 # initialize level parameters 

28 level_params = dict() 

29 level_params['restol'] = 1e-14 

30 

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

32 step_params = dict() 

33 

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

35 problem_params = dict() 

36 problem_params['cadv'] = 0.1 

37 problem_params['cs'] = 1.00 

38 problem_params['order_adv'] = 5 

39 problem_params['waveno'] = 5 

40 

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

42 sweeper_params = dict() 

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

44 sweeper_params['num_nodes'] = 3 

45 sweeper_params['do_coll_update'] = True 

46 

47 # initialize controller parameters 

48 controller_params = dict() 

49 controller_params['logger_level'] = 30 

50 controller_params['hook_class'] = dump_energy 

51 

52 # Fill description dictionary for easy hierarchy creation 

53 description = dict() 

54 description['problem_class'] = acoustic_1d_imex 

55 description['sweeper_class'] = imex_1st_order 

56 description['sweeper_params'] = sweeper_params 

57 

58 nsteps = np.zeros((3, 9)) 

59 nsteps[0, :] = [20, 30, 40, 50, 60, 70, 80, 90, 100] 

60 nsteps[1, :] = nsteps[0, :] 

61 nsteps[2, :] = nsteps[0, :] 

62 

63 for order in [3, 4, 5]: 

64 error = np.zeros(np.shape(nsteps)[1]) 

65 

66 # setup parameters "in time" 

67 t0 = 0 

68 Tend = 1.0 

69 

70 if order == 3: 

71 file = open(cwd + 'data/conv-data.txt', 'w') 

72 else: 

73 file = open(cwd + 'data/conv-data.txt', 'a') 

74 

75 step_params['maxiter'] = order 

76 description['step_params'] = step_params 

77 

78 for ii in range(0, np.shape(nsteps)[1]): 

79 ns = nsteps[order - 3, ii] 

80 if (order == 3) or (order == 4): 

81 problem_params['nvars'] = [(2, int(2 * ns))] 

82 elif order == 5: 

83 problem_params['nvars'] = [(2, int(2 * ns))] 

84 

85 description['problem_params'] = problem_params 

86 

87 dt = Tend / float(ns) 

88 

89 level_params['dt'] = dt 

90 description['level_params'] = level_params 

91 

92 # instantiate the controller 

93 controller = controller_nonMPI( 

94 num_procs=num_procs, controller_params=controller_params, description=description 

95 ) 

96 # get initial values on finest level 

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

98 uinit = P.u_exact(t0) 

99 if ii == 0: 

100 print("Time step: %4.2f" % dt) 

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

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

103 

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

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

106 

107 # compute exact solution and compare 

108 uex = P.u_exact(Tend) 

109 

110 error[ii] = np.linalg.norm(uex - uend, np.inf) / np.linalg.norm(uex, np.inf) 

111 file.write(str(order) + " " + str(ns) + " " + str(error[ii]) + "\n") 

112 

113 file.close() 

114 

115 for ii in range(0, np.shape(nsteps)[1]): 

116 print('error for nsteps= %s: %s' % (nsteps[order - 3, ii], error[ii])) 

117 

118 

119def plot_convergence(cwd=''): 

120 """ 

121 Plotting routine for the convergence data 

122 

123 Args: 

124 cwd (string): current workign directory 

125 """ 

126 

127 fs = 8 

128 order = np.array([]) 

129 nsteps = np.array([]) 

130 error = np.array([]) 

131 

132 file = open(cwd + 'data/conv-data.txt', 'r') 

133 while True: 

134 line = file.readline() 

135 if not line: 

136 break 

137 items = str.split(line, " ", 3) 

138 order = np.append(order, int(items[0])) 

139 nsteps = np.append(nsteps, int(float(items[1]))) 

140 error = np.append(error, float(items[2])) 

141 

142 assert np.size(order) == np.size(nsteps), 'Found different number of entries in order and nsteps' 

143 assert np.size(nsteps) == np.size(error), 'Found different number of entries in nsteps and error' 

144 

145 assert np.size(nsteps) % 3 == 0, 'Number of entries not a multiple of three, got %s' % np.size(nsteps) 

146 

147 N = int(np.size(nsteps) / 3) 

148 

149 error_plot = np.zeros((3, N)) 

150 nsteps_plot = np.zeros((3, N)) 

151 convline = np.zeros((3, N)) 

152 order_plot = np.zeros(3) 

153 

154 for ii in range(0, 3): 

155 order_plot[ii] = order[N * ii] 

156 for jj in range(0, N): 

157 error_plot[ii, jj] = error[N * ii + jj] 

158 nsteps_plot[ii, jj] = nsteps[N * ii + jj] 

159 convline[ii, jj] = ( 

160 error_plot[ii, 0] * (float(nsteps_plot[ii, 0]) / float(nsteps_plot[ii, jj])) ** order_plot[ii] 

161 ) 

162 

163 color = ['r', 'b', 'g'] 

164 shape = ['o', 'd', 's'] 

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

166 rcParams['pgf.rcfonts'] = False 

167 fig = plt.figure() 

168 for ii in range(0, 3): 

169 plt.loglog(nsteps_plot[ii, :], convline[ii, :], '-', color=color[ii]) 

170 plt.loglog( 

171 nsteps_plot[ii, :], 

172 error_plot[ii, :], 

173 shape[ii], 

174 markersize=fs, 

175 color=color[ii], 

176 label='p=' + str(int(order_plot[ii])), 

177 ) 

178 

179 plt.legend(loc='lower left', fontsize=fs, prop={'size': fs}) 

180 plt.xlabel('Number of time steps', fontsize=fs) 

181 plt.ylabel('Relative error', fontsize=fs, labelpad=2) 

182 plt.xlim([0.9 * np.min(nsteps_plot), 1.1 * np.max(nsteps_plot)]) 

183 plt.ylim([1e-5, 1e0]) 

184 plt.yticks([1e-5, 1e-4, 1e-3, 1e-2, 1e-1, 1e0], fontsize=fs) 

185 plt.xticks([20, 30, 40, 60, 80, 100], fontsize=fs) 

186 plt.gca().get_xaxis().get_major_formatter().labelOnlyBase = False 

187 plt.gca().get_xaxis().set_major_formatter(ScalarFormatter()) 

188 filename = 'data/convergence.png' 

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

190 

191 

192if __name__ == "__main__": 

193 compute_convergence_data() 

194 plot_convergence()