Coverage for pySDC/projects/AllenCahn_Bayreuth/run_simple_forcing_verification.py: 97%

1from argparse import ArgumentParser 

2import json 

3import glob 

4import numpy as np 

5from mpi4py import MPI 

6 

7import pySDC.helpers.plot_helper as plt_helper 

8import matplotlib.ticker as ticker 

9 

10from pySDC.helpers.stats_helper import get_sorted 

11 

12from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI 

13from pySDC.implementations.sweeper_classes.imex_1st_order import imex_1st_order 

14from pySDC.implementations.problem_classes.AllenCahn_MPIFFT import allencahn_imex, allencahn_imex_timeforcing 

15from pySDC.implementations.transfer_classes.TransferMesh_MPIFFT import fft_to_fft 

16 

17from pySDC.projects.AllenCahn_Bayreuth.AllenCahn_monitor import monitor 

18 

19 

20def run_simulation(name='', spectral=None, nprocs_space=None): 

21 """ 

22 A test program to do PFASST runs for the AC equation with different forcing 

23 

24 Args: 

25 name (str): name of the run, will be used to distinguish different setups 

26 spectral (bool): run in real or spectral space 

27 nprocs_space (int): number of processors in space (None if serial) 

28 """ 

29 

30 # set MPI communicator 

31 comm = MPI.COMM_WORLD 

32 

33 world_rank = comm.Get_rank() 

34 world_size = comm.Get_size() 

35 

36 # split world communicator to create space-communicators 

37 if nprocs_space is not None: 

38 color = int(world_rank / nprocs_space) 

39 else: 

40 color = int(world_rank / 1) 

41 space_comm = comm.Split(color=color) 

42 space_rank = space_comm.Get_rank() 

43 space_size = space_comm.Get_size() 

44 

45 assert world_size == space_size, 'This script cannot run parallel-in-time with MPI, only spatial parallelism' 

46 

47 # initialize level parameters 

48 level_params = dict() 

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

50 level_params['dt'] = 1e-03 

51 level_params['nsweeps'] = [1] 

52 

53 # initialize sweeper parameters 

54 sweeper_params = dict() 

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

56 sweeper_params['num_nodes'] = [3] 

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

58 sweeper_params['initial_guess'] = 'spread' 

59 

60 # initialize problem parameters 

61 problem_params = dict() 

62 problem_params['L'] = 1.0 

63 problem_params['nvars'] = [(128, 128), (32, 32)] 

64 problem_params['eps'] = [0.04] 

65 problem_params['radius'] = 0.25 

66 problem_params['comm'] = space_comm 

67 problem_params['init_type'] = 'circle' 

68 problem_params['spectral'] = spectral 

69 

70 if name == 'AC-test-constforce': 

71 problem_params['dw'] = [-23.59] 

72 

73 # initialize step parameters 

74 step_params = dict() 

75 step_params['maxiter'] = 50 

76 

77 # initialize controller parameters 

78 controller_params = dict() 

79 controller_params['logger_level'] = 30 if space_rank == 0 else 99 # set level depending on rank 

80 controller_params['hook_class'] = monitor 

81 controller_params['predict_type'] = 'pfasst_burnin' 

82 

83 # fill description dictionary for easy step instantiation 

84 description = dict() 

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

86 description['sweeper_class'] = imex_1st_order 

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

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

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

90 description['space_transfer_class'] = fft_to_fft 

91 

92 if name == 'AC-test-noforce' or name == 'AC-test-constforce': 

93 description['problem_class'] = allencahn_imex 

94 elif name == 'AC-test-timeforce': 

95 description['problem_class'] = allencahn_imex_timeforcing 

96 else: 

97 raise NotImplementedError(f'{name} is not implemented') 

98 

99 # set time parameters 

100 t0 = 0.0 

101 Tend = 32 * 0.001 

102 

103 if space_rank == 0: 

104 out = f'---------> Running {name} with spectral={spectral} and {space_size} process(es) in space...' 

105 print(out) 

106 

107 # instantiate controller 

108 controller = controller_nonMPI(num_procs=8, controller_params=controller_params, description=description) 

109 

110 # get initial values on finest level 

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

112 uinit = P.u_exact(t0) 

113 

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

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

116 

117 if space_rank == 0: 

118 # convert filtered statistics to list of computed radii, sorted by time 

119 computed_radii = get_sorted(stats, type='computed_radius', sortby='time') 

120 exact_radii = get_sorted(stats, type='exact_radius', sortby='time') 

121 computed_vol = get_sorted(stats, type='computed_volume', sortby='time') 

122 exact_vol = get_sorted(stats, type='exact_volume', sortby='time') 

123 

124 # print and store radii and error over time 

125 err_test = 0.0 

126 results = dict() 

127 for cr, er, cv, ev in zip(computed_radii, exact_radii, computed_vol, exact_vol): 

128 if name == 'AC-test-noforce': 

129 exrad = er[1] 

130 exvol = ev[1] 

131 else: 

132 exrad = computed_radii[0][1] 

133 exvol = computed_vol[0][1] 

134 if exrad > 0: 

135 errr = abs(cr[1] - exrad) / exrad 

136 errv = abs(cv[1] - exvol) / exvol 

137 else: 

138 errr = 1.0 

139 errv = 1.0 

140 if cr[0] == 0.025: 

141 err_test = errr 

142 out = f'Computed/exact/error radius for time {cr[0]:6.4f}: ' f'{cr[1]:8.6f} / {exrad:8.6f} / {errr:6.4e}' 

143 print(out) 

144 results[cr[0]] = (cr[1], exrad, errr, cv[1], exvol, errv) 

145 fname = f'./data/{name}_results.json' 

146 with open(fname, 'w') as fp: 

147 json.dump(results, fp, sort_keys=True, indent=4) 

148 

149 print() 

150 

151 # convert filtered statistics of iterations count, sorted by time 

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

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

154 out = f'Mean number of iterations: {niters:.4f}' 

155 print(out) 

156 

157 # get setup time 

158 timing = get_sorted(stats, type='timing_setup', sortby='time') 

159 out = f'Setup time: {timing[0][1]:.4f} sec.' 

160 print(out) 

161 

162 # get running time 

163 timing = get_sorted(stats, type='timing_run', sortby='time') 

164 out = f'Time to solution: {timing[0][1]:.4f} sec.' 

165 print(out) 

166 

167 out = '...Done <---------\n' 

168 print(out) 

169 

170 # Testing the output 

171 if name == 'AC-test-noforce': 

172 if spectral: 

173 exp_iters = 6.59375 

174 exp_err = 7.821e-02 

175 else: 

176 exp_iters = 7.8125 

177 exp_err = 7.85e-02 

178 elif name == 'AC-test-constforce': 

179 if spectral: 

180 exp_iters = 2.875 

181 exp_err = 4.678e-04 

182 else: 

183 exp_iters = 4.3125 

184 exp_err = 6.2384e-04 

185 elif name == 'AC-test-timeforce': 

186 if spectral: 

187 exp_iters = 1.65625 

188 exp_err = 6.2345e-04 

189 else: 

190 exp_iters = 2.40625 

191 exp_err = 6.2345e-04 

192 else: 

193 raise NotImplementedError(f'{name} is not implemented') 

194 

195 assert niters == exp_iters, f'Got deviating iteration counts of {niters} instead of {exp_iters}' 

196 assert err_test < exp_err, f'Got deviating errors of {err_test} instead of {exp_err}' 

197 

198 space_comm.Free() 

199 

200 

201def visualize_radii(): 

202 """ 

203 Routine to plot the radii of the runs vs. the exact radii 

204 """ 

205 

206 plt_helper.setup_mpl() 

207 

208 filelist = glob.glob('./data/*_results.json') 

209 

210 for file in filelist: 

211 # read in file with data 

212 with open(file, 'r') as fp: 

213 results = json.load(fp) 

214 

215 print(f'Working on {file}...') 

216 

217 # get times and radii 

218 xcoords = list(results) 

219 computed_radii = [v[0] for k, v in results.items()] 

220 exact_radii = [v[1] for k, v in results.items()] 

221 computed_vol = [v[3] for k, v in results.items()] 

222 exact_vol = [v[4] for k, v in results.items()] 

223 

224 # compute bound for y-axis 

225 max_rad = max(max(computed_radii), max(exact_radii)) 

226 max_vol = max(max(computed_vol), max(exact_vol)) 

227 

228 # set up plot for radii 

229 fig, ax = plt_helper.newfig(textwidth=238.96, scale=1.0) 

230 

231 # and plot 

232 ax.plot(xcoords, computed_radii, label='Computed radius') 

233 ax.plot(xcoords, exact_radii, color='k', linestyle='--', linewidth=1, label='Exact radius') 

234 

235 # beautify and save plot 

236 ax.set_ylim([-0.01, max_rad * 1.1]) 

237 ax.yaxis.set_major_formatter(ticker.FormatStrFormatter('%1.2f')) 

238 # ax.xaxis.set_major_formatter(ticker.FormatStrFormatter('%1.2f')) 

239 ax.xaxis.set_major_locator(ticker.MultipleLocator(4)) 

240 ax.set_ylabel('radius') 

241 ax.set_xlabel('time') 

242 ax.grid() 

243 ax.legend(loc=3) 

244 # ax.set_title(file.split('/')[-1].replace('_results.json', '')) 

245 f = file.replace('_results.json', '_radii') 

246 plt_helper.savefig(f) 

247 

248 # test if all went well 

249 assert glob.glob(f'{f}.pdf'), 'ERROR: plotting did not create PDF file' 

250 # assert glob.glob(f'{f}.pgf'), 'ERROR: plotting did not create PGF file' 

251 assert glob.glob(f'{f}.png'), 'ERROR: plotting did not create PNG file' 

252 

253 # set up plot for volumes 

254 fig, ax = plt_helper.newfig(textwidth=238.96, scale=1.0) 

255 

256 # and plot 

257 ax.plot(xcoords, computed_vol, label='Computed volume') 

258 ax.plot(xcoords, exact_vol, color='k', linestyle='--', linewidth=1, label='Exact volume') 

259 

260 # beautify and save plot 

261 ax.set_ylim([-0.01, max_vol * 1.1]) 

262 ax.yaxis.set_major_formatter(ticker.FormatStrFormatter('%1.2f')) 

263 # ax.xaxis.set_major_formatter(ticker.FormatStrFormatter('%1.2f')) 

264 ax.xaxis.set_major_locator(ticker.MultipleLocator(4)) 

265 ax.set_ylabel('radius') 

266 ax.set_xlabel('time') 

267 ax.grid() 

268 ax.legend(loc=3) 

269 # ax.set_title(file.split('/')[-1].replace('_results.json', '')) 

270 f = file.replace('_results.json', '_volume') 

271 plt_helper.savefig(f) 

272 

273 # test if all went well 

274 assert glob.glob(f'{f}.pdf'), 'ERROR: plotting did not create PDF file' 

275 # assert glob.glob(f'{f}.pgf'), 'ERROR: plotting did not create PGF file' 

276 assert glob.glob(f'{f}.png'), 'ERROR: plotting did not create PNG file' 

277 

278 

279def main(nprocs_space=None): 

280 """ 

281 Little helper routine to run the whole thing 

282 

283 Args: 

284 nprocs_space (int): number of processors in space (None if serial) 

285 

286 """ 

287 name_list = ['AC-test-noforce', 'AC-test-constforce', 'AC-test-timeforce'] 

288 

289 for name in name_list: 

290 run_simulation(name=name, spectral=False, nprocs_space=nprocs_space) 

291 run_simulation(name=name, spectral=True, nprocs_space=nprocs_space) 

292 

293 

294if __name__ == "__main__": 

295 # Add parser to get number of processors in space (have to do this here to enable automatic testing) 

296 parser = ArgumentParser() 

297 parser.add_argument("-n", "--nprocs_space", help='Specifies the number of processors in space', type=int) 

298 args = parser.parse_args() 

299 

300 main(nprocs_space=args.nprocs_space) 

301 visualize_radii()