Coverage for pySDC/projects/AsympConv/PFASST_conv_tests.py: 100%

1import os 

2import pickle 

3 

4import matplotlib 

5import numpy as np 

6 

7matplotlib.use('Agg') 

8import matplotlib.pyplot as plt 

9 

10from pySDC.implementations.problem_classes.HeatEquation_ND_FD import heatNd_unforced 

11from pySDC.implementations.problem_classes.AdvectionEquation_ND_FD import advectionNd 

12from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit 

13from pySDC.implementations.transfer_classes.TransferMesh import mesh_to_mesh 

14from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI 

15 

16from pySDC.helpers.stats_helper import get_sorted 

17 

18 

19def main(): 

20 """ 

21 Main driver running diffusion and advection tests 

22 """ 

23 nsweeps = 3 

24 run_diffusion(nsweeps=nsweeps) 

25 run_advection(nsweeps=nsweeps) 

26 plot_results(nsweeps=nsweeps) 

27 

28 nsweeps = 2 

29 run_diffusion(nsweeps=nsweeps) 

30 run_advection(nsweeps=nsweeps) 

31 plot_results(nsweeps=nsweeps) 

32 

33 

34def run_diffusion(nsweeps): 

35 """ 

36 A simple test program to test PFASST convergence for the heat equation with random initial data 

37 

38 Args: 

39 nsweeps: number of fine sweeps to perform 

40 """ 

41 

42 # initialize level parameters 

43 level_params = dict() 

44 level_params['restol'] = 1e-08 

45 level_params['dt'] = 0.25 

46 level_params['nsweeps'] = [nsweeps, 1] 

47 

48 # initialize sweeper parameters 

49 sweeper_params = dict() 

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

51 sweeper_params['num_nodes'] = [3] 

52 sweeper_params['QI'] = ['LU'] 

53 sweeper_params['initial_guess'] = 'zero' 

54 

55 # initialize problem parameters 

56 problem_params = dict() 

57 problem_params['freq'] = 1 # frequency for the test value 

58 problem_params['nvars'] = [127, 63] # number of degrees of freedom for each level 

59 problem_params['bc'] = 'dirichlet-zero' # boundary conditions 

60 

61 # initialize step parameters 

62 step_params = dict() 

63 step_params['maxiter'] = 50 

64 

65 # initialize space transfer parameters 

66 space_transfer_params = dict() 

67 space_transfer_params['rorder'] = 2 

68 space_transfer_params['iorder'] = 2 

69 space_transfer_params['periodic'] = False 

70 

71 # initialize controller parameters 

72 controller_params = dict() 

73 controller_params['logger_level'] = 30 

74 

75 # fill description dictionary for easy step instantiation 

76 description = dict() 

77 description['problem_class'] = heatNd_unforced # pass problem class 

78 description['sweeper_class'] = generic_implicit # pass sweeper (see part B) 

79 description['sweeper_params'] = sweeper_params # pass sweeper parameters 

80 description['level_params'] = level_params # pass level parameters 

81 description['step_params'] = step_params # pass step parameters 

82 description['space_transfer_class'] = mesh_to_mesh # pass spatial transfer class 

83 description['space_transfer_params'] = space_transfer_params # pass paramters for spatial transfer 

84 

85 # set time parameters 

86 t0 = 0.0 

87 Tend = 4 * level_params['dt'] 

88 

89 # set up number of parallel time-steps to run PFASST with 

90 num_proc = 4 

91 

92 results = dict() 

93 

94 for i in range(-3, 10): 

95 ratio = level_params['dt'] / (1.0 / (problem_params['nvars'][0] + 1)) ** 2 

96 

97 problem_params['nu'] = 10.0**i / ratio # diffusion coefficient 

98 description['problem_params'] = problem_params # pass problem parameters 

99 

100 out = 'Working on c = %6.4e' % problem_params['nu'] 

101 print(out) 

102 cfl = ratio * problem_params['nu'] 

103 out = ' CFL number: %4.2e' % cfl 

104 print(out) 

105 

106 # instantiate controller 

107 controller = controller_nonMPI(num_procs=num_proc, controller_params=controller_params, description=description) 

108 

109 # get initial values on finest level 

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

111 uinit = P.u_exact(t0) 

112 

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

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

115 

116 # filter statistics by type (number of iterations) 

117 iter_counts = get_sorted(stats, type='niter', sortby='time') 

118 

119 niters = np.array([item[1] for item in iter_counts]) 

120 

121 out = ' Mean number of iterations: %4.2f' % np.mean(niters) 

122 print(out) 

123 

124 if nsweeps == 3 and (i == -3 or i == 9): 

125 assert np.mean(niters) <= 2, 'ERROR: too much iterations for diffusive asymptotics, got %s' % np.mean( 

126 niters 

127 ) 

128 

129 results[cfl] = np.mean(niters) 

130 

131 fname = 'data/results_conv_diffusion_NS' + str(nsweeps) + '.pkl' 

132 file = open(fname, 'wb') 

133 pickle.dump(results, file) 

134 file.close() 

135 

136 assert os.path.isfile(fname), 'ERROR: pickle did not create file' 

137 

138 

139def run_advection(nsweeps): 

140 """ 

141 A simple test program to test PFASST convergence for the periodic advection equation 

142 

143 Args: 

144 nsweeps: number of fine sweeps to perform 

145 """ 

146 

147 # initialize level parameters 

148 level_params = dict() 

149 level_params['restol'] = 1e-08 

150 level_params['dt'] = 0.25 

151 level_params['nsweeps'] = [nsweeps, 1] 

152 

153 # initialize sweeper parameters 

154 sweeper_params = dict() 

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

156 sweeper_params['num_nodes'] = [3] 

157 sweeper_params['QI'] = ['LU'] # For the IMEX sweeper, the LU-trick can be activated for the implicit part 

158 sweeper_params['initial_guess'] = 'zero' 

159 

160 # initialize problem parameters 

161 problem_params = dict() 

162 problem_params['freq'] = 64 # frequency for the test value 

163 problem_params['nvars'] = [128, 64] # number of degrees of freedom for each level 

164 problem_params['order'] = 2 

165 problem_params['stencil_type'] = 'center' 

166 problem_params['bc'] = 'periodic' # boundary conditions 

167 

168 # initialize step parameters 

169 step_params = dict() 

170 step_params['maxiter'] = 50 

171 

172 # initialize space transfer parameters 

173 space_transfer_params = dict() 

174 space_transfer_params['rorder'] = 2 

175 space_transfer_params['iorder'] = 2 

176 space_transfer_params['periodic'] = True 

177 

178 # initialize controller parameters 

179 controller_params = dict() 

180 controller_params['logger_level'] = 30 

181 

182 # fill description dictionary for easy step instantiation 

183 description = dict() 

184 description['problem_class'] = advectionNd # pass problem class 

185 description['sweeper_class'] = generic_implicit # pass sweeper (see part B) 

186 description['sweeper_params'] = sweeper_params # pass sweeper parameters 

187 description['level_params'] = level_params # pass level parameters 

188 description['step_params'] = step_params # pass step parameters 

189 description['space_transfer_class'] = mesh_to_mesh # pass spatial transfer class 

190 description['space_transfer_params'] = space_transfer_params # pass paramters for spatial transfer 

191 

192 # set time parameters 

193 t0 = 0.0 

194 Tend = 4 * level_params['dt'] 

195 

196 # set up number of parallel time-steps to run PFASST with 

197 num_proc = 4 

198 

199 results = dict() 

200 

201 for i in range(-3, 10): 

202 ratio = level_params['dt'] / (1.0 / (problem_params['nvars'][0] + 1)) 

203 

204 problem_params['c'] = 10.0**i / ratio # diffusion coefficient 

205 description['problem_params'] = problem_params # pass problem parameters 

206 

207 out = 'Working on nu = %6.4e' % problem_params['c'] 

208 print(out) 

209 cfl = ratio * problem_params['c'] 

210 out = ' CFL number: %4.2e' % cfl 

211 print(out) 

212 

213 # instantiate controller 

214 controller = controller_nonMPI(num_procs=num_proc, controller_params=controller_params, description=description) 

215 

216 # get initial values on finest level 

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

218 uinit = P.u_exact(t0) 

219 

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

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

222 

223 # filter statistics by type (number of iterations) 

224 iter_counts = get_sorted(stats, type='niter', sortby='time') 

225 

226 niters = np.array([item[1] for item in iter_counts]) 

227 

228 out = ' Mean number of iterations: %4.2f' % np.mean(niters) 

229 print(out) 

230 

231 if nsweeps == 3 and (i == -3 or i == 9): 

232 assert np.mean(niters) <= 2, 'ERROR: too much iterations for advective asymptotics, got %s' % np.mean( 

233 niters 

234 ) 

235 

236 results[cfl] = np.mean(niters) 

237 

238 fname = 'data/results_conv_advection_NS' + str(nsweeps) + '.pkl' 

239 file = open(fname, 'wb') 

240 pickle.dump(results, file) 

241 file.close() 

242 

243 assert os.path.isfile(fname), 'ERROR: pickle did not create file' 

244 

245 

246def plot_results(nsweeps): 

247 """ 

248 Plotting routine for iteration counts 

249 

250 Args: 

251 nsweeps: number of fine sweeps used 

252 """ 

253 

254 fname = 'data/results_conv_diffusion_NS' + str(nsweeps) + '.pkl' 

255 file = open(fname, 'rb') 

256 results_diff = pickle.load(file) 

257 file.close() 

258 

259 fname = 'data/results_conv_advection_NS' + str(nsweeps) + '.pkl' 

260 file = open(fname, 'rb') 

261 results_adv = pickle.load(file) 

262 file.close() 

263 

264 xvalues_diff = sorted(results_diff.keys()) 

265 niter_diff = [] 

266 for x in xvalues_diff: 

267 niter_diff.append(results_diff[x]) 

268 

269 xvalues_adv = sorted(results_adv.keys()) 

270 niter_adv = [] 

271 for x in xvalues_adv: 

272 niter_adv.append(results_adv[x]) 

273 

274 # set up plotting parameters 

275 params = { 

276 'legend.fontsize': 20, 

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

278 'axes.labelsize': 20, 

279 'axes.titlesize': 20, 

280 'xtick.labelsize': 16, 

281 'ytick.labelsize': 16, 

282 'lines.linewidth': 3, 

283 } 

284 plt.rcParams.update(params) 

285 

286 # set up figure 

287 plt.figure() 

288 plt.xlabel(r'$\mu$') 

289 plt.ylabel('no. of iterations') 

290 plt.xlim(min(xvalues_diff + xvalues_adv) / 10.0, max(xvalues_diff + xvalues_adv) * 10.0) 

291 plt.ylim(min(niter_diff + niter_adv) - 1, max(niter_diff + niter_adv) + 1) 

292 plt.grid() 

293 

294 # plot 

295 plt.semilogx(xvalues_diff, niter_diff, 'r-', marker='s', markersize=10, label='diffusion') 

296 plt.semilogx(xvalues_adv, niter_adv, 'b-', marker='o', markersize=10, label='advection') 

297 

298 plt.legend(loc=1, ncol=1, numpoints=1) 

299 

300 # plt.show() 

301 # save plot, beautify 

302 fname = 'data/conv_test_niter_NS' + str(nsweeps) + '.png' 

303 plt.savefig(fname, bbox_inches='tight') 

304 

305 assert os.path.isfile(fname), 'ERROR: plotting did not create file' 

306 

307 

308if __name__ == "__main__": 

309 main()