Coverage for pySDC/projects/Resilience/paper_plots.py: 0%

1# script to make pretty plots for papers or talks 

2import numpy as np 

3import matplotlib as mpl 

4import matplotlib.pyplot as plt 

5from pySDC.projects.Resilience.fault_stats import ( 

6 FaultStats, 

7 run_Lorenz, 

8 run_Schroedinger, 

9 run_vdp, 

10 run_quench, 

11 run_AC, 

12 run_RBC, 

13 RECOVERY_THRESH_ABS, 

14) 

15from pySDC.projects.Resilience.strategies import ( 

16 BaseStrategy, 

17 AdaptivityStrategy, 

18 IterateStrategy, 

19 HotRodStrategy, 

20 DIRKStrategy, 

21 ERKStrategy, 

22 AdaptivityPolynomialError, 

23 cmap, 

24) 

25from pySDC.helpers.plot_helper import setup_mpl, figsize_by_journal 

26from pySDC.helpers.stats_helper import get_sorted 

27 

28 

29cm = 1 / 2.5 

30TEXTWIDTH = 11.9446244611 * cm 

31JOURNAL = 'Springer_Numerical_Algorithms' 

32BASE_PATH = 'data/paper' 

33 

34 

35def get_stats(problem, path='data/stats-jusuf', num_procs=1, strategy_type='SDC'): 

36 """ 

37 Create a FaultStats object for a given problem to use for the plots. 

38 Note that the statistics need to be already generated somewhere else, this function will only load them. 

39 

40 Args: 

41 problem (function): A problem to run 

42 path (str): Path to the associated stats for the problem 

43 

44 Returns: 

45 FaultStats: Object to analyse resilience statistics from 

46 """ 

47 if strategy_type == 'SDC': 

48 strategies = [BaseStrategy(), AdaptivityStrategy(), IterateStrategy()] 

49 if JOURNAL not in ['JSC_beamer']: 

50 strategies += [HotRodStrategy(), AdaptivityPolynomialError()] 

51 elif strategy_type == 'RK': 

52 strategies = [DIRKStrategy()] 

53 if problem.__name__ in ['run_Lorenz', 'run_vdp']: 

54 strategies += [ERKStrategy()] 

55 

56 stats_analyser = FaultStats( 

57 prob=problem, 

58 strategies=strategies, 

59 faults=[False, True], 

60 reload=True, 

61 recovery_thresh=1.1, 

62 recovery_thresh_abs=RECOVERY_THRESH_ABS.get(problem, 0), 

63 mode='default', 

64 stats_path=path, 

65 num_procs=num_procs, 

66 ) 

67 stats_analyser.get_recovered() 

68 return stats_analyser 

69 

70 

71def my_setup_mpl(**kwargs): 

72 setup_mpl(reset=True, font_size=8) 

73 mpl.rcParams.update({'lines.markersize': 6}) 

74 

75 

76def savefig(fig, name, format='pdf', tight_layout=True): # pragma: no cover 

77 """ 

78 Save a figure to some predefined location. 

79 

80 Args: 

81 fig (Matplotlib.Figure): The figure of the plot 

82 name (str): The name of the plot 

83 tight_layout (bool): Apply tight layout or leave as is 

84 Returns: 

85 None 

86 """ 

87 if tight_layout: 

88 fig.tight_layout() 

89 path = f'{BASE_PATH}/{name}.{format}' 

90 fig.savefig(path, bbox_inches='tight', transparent=True, dpi=200) 

91 print(f'saved "{path}"') 

92 

93 

94def analyse_resilience(problem, path='data/stats', **kwargs): # pragma: no cover 

95 """ 

96 Generate some stats for resilience / load them if already available and make some plots. 

97 

98 Args: 

99 problem (function): A problem to run 

100 path (str): Path to the associated stats for the problem 

101 

102 Returns: 

103 None 

104 """ 

105 

106 stats_analyser = get_stats(problem, path) 

107 stats_analyser.get_recovered() 

108 

109 strategy = IterateStrategy() 

110 not_fixed = stats_analyser.get_mask(strategy=strategy, key='recovered', val=False) 

111 not_overflow = stats_analyser.get_mask(strategy=strategy, key='bit', val=1, op='uneq', old_mask=not_fixed) 

112 stats_analyser.print_faults(not_overflow) 

113 

114 compare_strategies(stats_analyser, **kwargs) 

115 plot_recovery_rate(stats_analyser, **kwargs) 

116 

117 

118def compare_strategies(stats_analyser, **kwargs): # pragma: no cover 

119 """ 

120 Make a plot showing local error and iteration number of time for all strategies 

121 

122 Args: 

123 stats_analyser (FaultStats): Fault stats object, which contains some stats 

124 

125 Returns: 

126 None 

127 """ 

128 my_setup_mpl() 

129 fig, ax = plt.subplots(figsize=(TEXTWIDTH, 5 * cm)) 

130 stats_analyser.compare_strategies(ax=ax) 

131 savefig(fig, 'compare_strategies', **kwargs) 

132 

133 

134def plot_recovery_rate(stats_analyser, **kwargs): # pragma: no cover 

135 """ 

136 Make a plot showing recovery rate for all faults and only for those that can be recovered. 

137 

138 Args: 

139 stats_analyser (FaultStats): Fault stats object, which contains some stats 

140 

141 Returns: 

142 None 

143 """ 

144 my_setup_mpl() 

145 fig, axs = plt.subplots(1, 2, figsize=(TEXTWIDTH, 5 * cm), sharex=True, sharey=True) 

146 stats_analyser.plot_things_per_things( 

147 'recovered', 

148 'bit', 

149 False, 

150 op=stats_analyser.rec_rate, 

151 args={'ylabel': 'recovery rate'}, 

152 plotting_args={'markevery': 5}, 

153 ax=axs[0], 

154 ) 

155 plot_recovery_rate_recoverable_only(stats_analyser, fig, axs[1], ylabel='') 

156 axs[0].get_legend().remove() 

157 axs[0].set_title('All faults') 

158 axs[1].set_title('Only recoverable faults') 

159 axs[0].set_ylim((-0.05, 1.05)) 

160 savefig(fig, 'recovery_rate_compared', **kwargs) 

161 

162 

163def plot_recovery_rate_recoverable_only(stats_analyser, fig, ax, **kwargs): # pragma: no cover 

164 """ 

165 Plot the recovery rate considering only faults that can be recovered theoretically. 

166 

167 Args: 

168 stats_analyser (FaultStats): Fault stats object, which contains some stats 

169 fig (matplotlib.pyplot.figure): Figure in which to plot 

170 ax (matplotlib.pyplot.axes): Somewhere to plot 

171 

172 Returns: 

173 None 

174 """ 

175 for i in range(len(stats_analyser.strategies)): 

176 fixable = stats_analyser.get_fixable_faults_only(strategy=stats_analyser.strategies[i]) 

177 

178 stats_analyser.plot_things_per_things( 

179 'recovered', 

180 'bit', 

181 False, 

182 op=stats_analyser.rec_rate, 

183 mask=fixable, 

184 args={**kwargs}, 

185 ax=ax, 

186 fig=fig, 

187 strategies=[stats_analyser.strategies[i]], 

188 plotting_args={'markevery': 5}, 

189 ) 

190 

191 

192def compare_recovery_rate_problems(target='resilience', **kwargs): # pragma: no cover 

193 """ 

194 Compare the recovery rate for various problems. 

195 Only faults that can be recovered are shown. 

196 

197 Returns: 

198 None 

199 """ 

200 if target == 'resilience': 

201 problems = [run_Lorenz, run_Schroedinger, run_AC, run_RBC] 

202 titles = ['Lorenz', r'Schr\"odinger', 'Allen-Cahn', 'Rayleigh-Benard'] 

203 elif target == 'thesis': 

204 problems = [run_vdp, run_Lorenz, run_AC, run_RBC] # TODO: swap in Gray-Scott 

205 titles = ['Van der Pol', 'Lorenz', 'Allen-Cahn', 'Rayleigh-Benard'] 

206 else: 

207 raise NotImplementedError() 

208 

209 stats = [get_stats(problem, **kwargs) for problem in problems] 

210 

211 my_setup_mpl() 

212 fig, axs = plt.subplots(2, 2, figsize=figsize_by_journal(JOURNAL, 1, 0.8), sharey=True) 

213 [ 

214 plot_recovery_rate_recoverable_only(stats[i], fig, axs.flatten()[i], ylabel='', title=titles[i]) 

215 for i in range(len(stats)) 

216 ] 

217 

218 for ax in axs.flatten(): 

219 ax.get_legend().remove() 

220 

221 if kwargs.get('strategy_type', 'SDC') == 'SDC': 

222 axs[1, 1].legend(frameon=False, loc="lower right") 

223 else: 

224 axs[0, 1].legend(frameon=False, loc="lower right") 

225 axs[0, 0].set_ylim((-0.05, 1.05)) 

226 axs[1, 0].set_ylabel('recovery rate') 

227 axs[0, 0].set_ylabel('recovery rate') 

228 

229 name = '' 

230 for key, val in kwargs.items(): 

231 name = f'{name}_{key}-{val}' 

232 

233 savefig(fig, f'compare_equations{name}.pdf') 

234 

235 

236def plot_adaptivity_stuff(): # pragma: no cover 

237 """ 

238 Plot the solution for a van der Pol problem as well as the local error and cost associated with the base scheme and 

239 adaptivity in k and dt in order to demonstrate that adaptivity is useful. 

240 

241 Returns: 

242 None 

243 """ 

244 from pySDC.implementations.hooks.log_errors import LogLocalErrorPostStep 

245 from pySDC.implementations.hooks.log_work import LogWork 

246 from pySDC.projects.Resilience.hook import LogData 

247 import pickle 

248 

249 my_setup_mpl() 

250 scale = 0.5 if JOURNAL == 'JSC_beamer' else 1.0 

251 fig, axs = plt.subplots(3, 1, figsize=figsize_by_journal(JOURNAL, scale, 1), sharex=True, sharey=False) 

252 

253 def plot_error(stats, ax, iter_ax, strategy, **kwargs): 

254 """ 

255 Plot global error and cumulative sum of iterations 

256 

257 Args: 

258 stats (dict): Stats from pySDC run 

259 ax (Matplotlib.pyplot.axes): Somewhere to plot the error 

260 iter_ax (Matplotlib.pyplot.axes): Somewhere to plot the iterations 

261 strategy (pySDC.projects.Resilience.fault_stats.Strategy): The resilience strategy 

262 

263 Returns: 

264 None 

265 """ 

266 markevery = 1 if type(strategy) in [AdaptivityStrategy, AdaptivityPolynomialError] else 10000 

267 e = stats['e_local_post_step'] 

268 ax.plot([me[0] for me in e], [me[1] for me in e], markevery=markevery, **strategy.style, **kwargs) 

269 k = stats['work_newton'] 

270 iter_ax.plot( 

271 [me[0] for me in k], np.cumsum([me[1] for me in k]), **strategy.style, markevery=markevery, **kwargs 

272 ) 

273 ax.set_yscale('log') 

274 ax.set_ylabel('local error') 

275 iter_ax.set_ylabel(r'Newton iterations') 

276 

277 run = False 

278 for strategy in [BaseStrategy, IterateStrategy, AdaptivityStrategy, AdaptivityPolynomialError]: 

279 S = strategy(newton_inexactness=False) 

280 desc = S.get_custom_description(problem=run_vdp, num_procs=1) 

281 desc['problem_params']['mu'] = 1000 

282 desc['problem_params']['u0'] = (1.1, 0) 

283 if strategy in [AdaptivityStrategy, BaseStrategy]: 

284 desc['step_params']['maxiter'] = 5 

285 if strategy in [BaseStrategy, IterateStrategy]: 

286 desc['level_params']['dt'] = 1e-4 

287 desc['sweeper_params']['QI'] = 'LU' 

288 if strategy in [IterateStrategy]: 

289 desc['step_params']['maxiter'] = 99 

290 desc['level_params']['restol'] = 1e-10 

291 

292 path = f'./data/adaptivity_paper_plot_data_{strategy.__name__}.pickle' 

293 if run: 

294 stats, _, _ = run_vdp( 

295 custom_description=desc, 

296 Tend=20, 

297 hook_class=[LogLocalErrorPostStep, LogWork, LogData], 

298 custom_controller_params={'logger_level': 15}, 

299 ) 

300 

301 data = { 

302 'u': get_sorted(stats, type='u', recomputed=False), 

303 'e_local_post_step': get_sorted(stats, type='e_local_post_step', recomputed=False), 

304 'work_newton': get_sorted(stats, type='work_newton', recomputed=None), 

305 } 

306 with open(path, 'wb') as file: 

307 pickle.dump(data, file) 

308 else: 

309 with open(path, 'rb') as file: 

310 data = pickle.load(file) 

311 

312 plot_error(data, axs[1], axs[2], strategy()) 

313 

314 if strategy == BaseStrategy or True: 

315 u = data['u'] 

316 axs[0].plot([me[0] for me in u], [me[1][0] for me in u], color='black', label=r'$u$') 

317 

318 axs[2].set_xlabel(r'$t$') 

319 axs[0].set_ylabel('solution') 

320 axs[2].legend(frameon=JOURNAL == 'JSC_beamer') 

321 axs[1].legend(frameon=True) 

322 axs[2].set_yscale('log') 

323 savefig(fig, 'adaptivity') 

324 

325 

326def plot_fault_vdp(bit=0): # pragma: no cover 

327 """ 

328 Make a plot showing the impact of a fault on van der Pol without any resilience. 

329 The faults are inserted in the last iteration in the last node in u_t such that you can best see the impact. 

330 

331 Args: 

332 bit (int): The bit that you want to flip 

333 

334 Returns: 

335 None 

336 """ 

337 from pySDC.projects.Resilience.fault_stats import ( 

338 FaultStats, 

339 BaseStrategy, 

340 ) 

341 from pySDC.projects.Resilience.hook import LogData 

342 

343 stats_analyser = FaultStats( 

344 prob=run_vdp, 

345 strategies=[BaseStrategy()], 

346 faults=[False, True], 

347 reload=True, 

348 recovery_thresh=1.1, 

349 num_procs=1, 

350 mode='combination', 

351 ) 

352 

353 my_setup_mpl() 

354 fig, ax = plt.subplots(figsize=figsize_by_journal(JOURNAL, 0.8, 0.5)) 

355 colors = ['blue', 'red', 'magenta'] 

356 ls = ['--', '-'] 

357 markers = ['*', '^'] 

358 do_faults = [False, True] 

359 superscripts = ['*', ''] 

360 subscripts = ['', 't', ''] 

361 

362 run = 779 + 12 * bit # for faults in u_t 

363 # run = 11 + 12 * bit # for faults in u 

364 

365 for i in range(len(do_faults)): 

366 stats, controller, Tend = stats_analyser.single_run( 

367 strategy=BaseStrategy(), 

368 run=run, 

369 faults=do_faults[i], 

370 hook_class=[LogData], 

371 ) 

372 u = get_sorted(stats, type='u') 

373 faults = get_sorted(stats, type='bitflip') 

374 for j in [0, 1]: 

375 ax.plot( 

376 [me[0] for me in u], 

377 [me[1][j] for me in u], 

378 ls=ls[i], 

379 color=colors[j], 

380 label=rf'$u^{ {superscripts[i]}} _{ {subscripts[j]}} $', 

381 marker=markers[j], 

382 markevery=60, 

383 ) 

384 for idx in range(len(faults)): 

385 ax.axvline(faults[idx][0], color='black', label='Fault', ls=':') 

386 print( 

387 f'Fault at t={faults[idx][0]:.2e}, iter={faults[idx][1][1]}, node={faults[idx][1][2]}, space={faults[idx][1][3]}, bit={faults[idx][1][4]}' 

388 ) 

389 ax.set_title(f'Fault in bit {faults[idx][1][4]}') 

390 

391 ax.legend(frameon=True, loc='lower left') 

392 ax.set_xlabel(r'$t$') 

393 savefig(fig, f'fault_bit_{bit}') 

394 

395 

396def plot_fault_Lorenz(bit=0): # pragma: no cover 

397 """ 

398 Make a plot showing the impact of a fault on the Lorenz attractor without any resilience. 

399 The faults are inserted in the last iteration in the last node in x such that you can best see the impact. 

400 

401 Args: 

402 bit (int): The bit that you want to flip 

403 

404 Returns: 

405 None 

406 """ 

407 from pySDC.projects.Resilience.fault_stats import ( 

408 FaultStats, 

409 BaseStrategy, 

410 ) 

411 from pySDC.projects.Resilience.hook import LogData 

412 

413 stats_analyser = FaultStats( 

414 prob=run_Lorenz, 

415 strategies=[BaseStrategy()], 

416 faults=[False, True], 

417 reload=True, 

418 recovery_thresh=1.1, 

419 num_procs=1, 

420 mode='combination', 

421 ) 

422 

423 strategy = BaseStrategy() 

424 

425 my_setup_mpl() 

426 fig, ax = plt.subplots(figsize=figsize_by_journal(JOURNAL, 0.8, 0.5)) 

427 colors = ['grey', strategy.color, 'magenta'] 

428 ls = ['--', '-'] 

429 markers = [None, strategy.marker] 

430 do_faults = [False, True] 

431 superscripts = ['*', ''] 

432 labels = ['x', 'x'] 

433 

434 run = 19 + 20 * bit 

435 

436 for i in range(len(do_faults)): 

437 stats, controller, Tend = stats_analyser.single_run( 

438 strategy=BaseStrategy(), 

439 run=run, 

440 faults=do_faults[i], 

441 hook_class=[LogData], 

442 ) 

443 u = get_sorted(stats, type='u') 

444 faults = get_sorted(stats, type='bitflip') 

445 ax.plot( 

446 [me[0] for me in u], 

447 [me[1][0] for me in u], 

448 ls=ls[i], 

449 color=colors[i], 

450 label=rf'${ {labels[i]}} ^{ {superscripts[i]}} $', 

451 marker=markers[i], 

452 markevery=500, 

453 ) 

454 for idx in range(len(faults)): 

455 ax.axvline(faults[idx][0], color='black', label='Fault', ls=':') 

456 print( 

457 f'Fault at t={faults[idx][0]:.2e}, iter={faults[idx][1][1]}, node={faults[idx][1][2]}, space={faults[idx][1][3]}, bit={faults[idx][1][4]}' 

458 ) 

459 ax.set_title(f'Fault in bit {faults[idx][1][4]}') 

460 

461 ax.legend(frameon=True, loc='lower left') 

462 ax.set_xlabel(r'$t$') 

463 savefig(fig, f'fault_bit_{bit}') 

464 

465 

466def plot_Lorenz_solution(): # pragma: no cover 

467 my_setup_mpl() 

468 

469 fig, axs = plt.subplots(1, 2, figsize=figsize_by_journal(JOURNAL, 1, 0.4), sharex=True) 

470 

471 strategy = BaseStrategy() 

472 desc = strategy.get_custom_description(run_Lorenz, num_procs=1) 

473 stats, controller, _ = run_Lorenz(custom_description=desc, Tend=strategy.get_Tend(run_Lorenz)) 

474 

475 u = get_sorted(stats, recomputed=False, type='u') 

476 

477 axs[0].plot([me[1][0] for me in u], [me[1][2] for me in u]) 

478 axs[0].set_ylabel('$z$') 

479 axs[0].set_xlabel('$x$') 

480 

481 axs[1].plot([me[1][0] for me in u], [me[1][1] for me in u]) 

482 axs[1].set_ylabel('$y$') 

483 axs[1].set_xlabel('$x$') 

484 

485 for ax in axs: 

486 ax.set_box_aspect(1.0) 

487 

488 path = 'data/paper/Lorenz_sol.pdf' 

489 fig.savefig(path, bbox_inches='tight', transparent=True, dpi=200) 

490 

491 

492def plot_quench_solution(): # pragma: no cover 

493 """ 

494 Plot the solution of Quench problem over time 

495 

496 Returns: 

497 None 

498 """ 

499 my_setup_mpl() 

500 if JOURNAL == 'JSC_beamer': 

501 fig, ax = plt.subplots(figsize=figsize_by_journal(JOURNAL, 0.5, 0.9)) 

502 else: 

503 fig, ax = plt.subplots(figsize=figsize_by_journal(JOURNAL, 1.0, 0.45)) 

504 

505 strategy = BaseStrategy() 

506 

507 custom_description = strategy.get_custom_description(run_quench, num_procs=1) 

508 

509 stats, controller, _ = run_quench(custom_description=custom_description, Tend=strategy.get_Tend(run_quench)) 

510 

511 prob = controller.MS[0].levels[0].prob 

512 

513 u = get_sorted(stats, type='u', recomputed=False) 

514 

515 ax.plot([me[0] for me in u], [max(me[1]) for me in u], color='black', label='$T$') 

516 ax.axhline(prob.u_thresh, label=r'$T_\mathrm{thresh}$', ls='--', color='grey', zorder=-1) 

517 ax.axhline(prob.u_max, label=r'$T_\mathrm{max}$', ls=':', color='grey', zorder=-1) 

518 

519 ax.set_xlabel(r'$t$') 

520 ax.legend(frameon=False) 

521 savefig(fig, 'quench_sol') 

522 

523 

524def plot_RBC_solution(): # pragma: no cover 

525 """ 

526 Plot solution of Rayleigh-Benard convection 

527 """ 

528 my_setup_mpl() 

529 

530 from mpl_toolkits.axes_grid1 import make_axes_locatable 

531 

532 plt.rcParams['figure.constrained_layout.use'] = True 

533 fig, axs = plt.subplots(2, 1, sharex=True, sharey=True, figsize=figsize_by_journal(JOURNAL, 1.0, 0.45)) 

534 caxs = [] 

535 divider = make_axes_locatable(axs[0]) 

536 caxs += [divider.append_axes('right', size='3%', pad=0.03)] 

537 divider2 = make_axes_locatable(axs[1]) 

538 caxs += [divider2.append_axes('right', size='3%', pad=0.03)] 

539 

540 from pySDC.projects.Resilience.RBC import RayleighBenard, PROBLEM_PARAMS 

541 

542 prob = RayleighBenard(**PROBLEM_PARAMS) 

543 

544 def _plot(t, ax, cax): 

545 u_hat = prob.u_exact(t) 

546 u = prob.itransform(u_hat) 

547 im = ax.pcolormesh(prob.X, prob.Z, u[prob.index('T')], rasterized=True, cmap='plasma') 

548 fig.colorbar(im, cax, label=f'$T(t={ {t}} )$') 

549 

550 _plot(0, axs[0], caxs[0]) 

551 _plot(21, axs[1], caxs[1]) 

552 

553 axs[1].set_xlabel('$x$') 

554 axs[0].set_ylabel('$z$') 

555 axs[1].set_ylabel('$z$') 

556 

557 savefig(fig, 'RBC_sol', tight_layout=False) 

558 

559 

560def plot_Schroedinger_solution(): # pragma: no cover 

561 from pySDC.implementations.problem_classes.NonlinearSchroedinger_MPIFFT import nonlinearschroedinger_imex 

562 

563 my_setup_mpl() 

564 if JOURNAL == 'JSC_beamer': 

565 raise NotImplementedError 

566 fig, ax = plt.subplots(figsize=figsize_by_journal(JOURNAL, 0.5, 0.9)) 

567 else: 

568 fig, axs = plt.subplots(1, 2, figsize=figsize_by_journal(JOURNAL, 1.0, 0.45), sharex=True, sharey=True) 

569 

570 from mpl_toolkits.axes_grid1 import make_axes_locatable 

571 

572 plt.rcParams['figure.constrained_layout.use'] = True 

573 cax = [] 

574 divider = make_axes_locatable(axs[0]) 

575 cax += [divider.append_axes('right', size='5%', pad=0.05)] 

576 divider2 = make_axes_locatable(axs[1]) 

577 cax += [divider2.append_axes('right', size='5%', pad=0.05)] 

578 

579 problem_params = dict() 

580 problem_params['nvars'] = (256, 256) 

581 problem_params['spectral'] = False 

582 problem_params['c'] = 1.0 

583 description = {'problem_params': problem_params} 

584 stats, _, _ = run_Schroedinger(Tend=1.0e0, custom_description=description) 

585 

586 P = nonlinearschroedinger_imex(**problem_params) 

587 u = get_sorted(stats, type='u') 

588 

589 im = axs[0].pcolormesh(*P.X, np.abs(u[0][1]), rasterized=True) 

590 im1 = axs[1].pcolormesh(*P.X, np.abs(u[-1][1]), rasterized=True) 

591 

592 fig.colorbar(im, cax=cax[0]) 

593 fig.colorbar(im1, cax=cax[1]) 

594 axs[0].set_title(r'$\|u(t=0)\|$') 

595 axs[1].set_title(r'$\|u(t=1)\|$') 

596 for ax in axs: 

597 ax.set_aspect(1) 

598 ax.set_xlabel('$x$') 

599 ax.set_ylabel('$y$') 

600 savefig(fig, 'Schroedinger_sol') 

601 

602 

603def plot_AC_solution(): # pragma: no cover 

604 from pySDC.projects.Resilience.AC import monitor 

605 

606 my_setup_mpl() 

607 if JOURNAL == 'JSC_beamer': 

608 raise NotImplementedError 

609 fig, ax = plt.subplots(figsize=figsize_by_journal(JOURNAL, 0.5, 0.9)) 

610 else: 

611 fig, axs = plt.subplots(1, 2, figsize=figsize_by_journal(JOURNAL, 1.0, 0.45)) 

612 

613 description = {'problem_params': {'nvars': (256, 256)}} 

614 stats, _, _ = run_AC(Tend=0.032, hook_class=monitor, custom_description=description) 

615 

616 u = get_sorted(stats, type='u') 

617 

618 computed_radius = get_sorted(stats, type='computed_radius') 

619 axs[1].plot([me[0] for me in computed_radius], [me[1] for me in computed_radius], ls='-') 

620 axs[1].axvline(0.025, ls=':', label=r'$t=0.025$', color='grey') 

621 axs[1].set_title('Radius over time') 

622 axs[1].set_xlabel('$t$') 

623 axs[1].legend(frameon=False) 

624 

625 im = axs[0].imshow(u[0][1], extent=(-0.5, 0.5, -0.5, 0.5)) 

626 fig.colorbar(im) 

627 axs[0].set_title(r'$u_0$') 

628 axs[0].set_xlabel('$x$') 

629 axs[0].set_ylabel('$y$') 

630 savefig(fig, 'AC_sol') 

631 

632 

633def plot_vdp_solution(): # pragma: no cover 

634 """ 

635 Plot the solution of van der Pol problem over time to illustrate the varying time scales. 

636 

637 Returns: 

638 None 

639 """ 

640 from pySDC.implementations.convergence_controller_classes.adaptivity import Adaptivity 

641 

642 my_setup_mpl() 

643 if JOURNAL == 'JSC_beamer': 

644 fig, ax = plt.subplots(figsize=figsize_by_journal(JOURNAL, 0.5, 0.9)) 

645 else: 

646 fig, ax = plt.subplots(figsize=figsize_by_journal(JOURNAL, 1.0, 0.33)) 

647 

648 custom_description = { 

649 'convergence_controllers': {Adaptivity: {'e_tol': 1e-7, 'dt_max': 1e0}}, 

650 'problem_params': {'mu': 1000, 'crash_at_maxiter': False}, 

651 'level_params': {'dt': 1e-3}, 

652 } 

653 

654 stats, _, _ = run_vdp(custom_description=custom_description, Tend=2000) 

655 

656 u = get_sorted(stats, type='u', recomputed=False) 

657 _u = np.array([me[1][0] for me in u]) 

658 _x = np.array([me[0] for me in u]) 

659 

660 x1 = _x[abs(_u - 1.1) < 1e-2][0] 

661 ax.plot(_x, _u, color='black') 

662 ax.axvspan(x1, x1 + 20, alpha=0.4) 

663 ax.set_ylabel(r'$u$') 

664 ax.set_xlabel(r'$t$') 

665 savefig(fig, 'vdp_sol') 

666 

667 

668def work_precision(): # pragma: no cover 

669 from pySDC.projects.Resilience.work_precision import ( 

670 all_problems, 

671 ) 

672 

673 all_params = { 

674 'record': False, 

675 'work_key': 't', 

676 'precision_key': 'e_global_rel', 

677 'plotting': True, 

678 'base_path': 'data/paper', 

679 } 

680 

681 for mode in ['compare_strategies', 'parallel_efficiency', 'RK_comp']: 

682 all_problems(**all_params, mode=mode) 

683 all_problems(**{**all_params, 'work_key': 'param'}, mode='compare_strategies') 

684 

685 

686def plot_recovery_rate_per_acceptance_threshold(problem): # pragma no cover 

687 stats_analyser = get_stats(problem) 

688 

689 stats_analyser.plot_recovery_thresholds(thresh_range=np.linspace(0.5, 1.5, 1000), recoverable_only=True) 

690 

691 

692def make_plots_for_TIME_X_website(): # pragma: no cover 

693 global JOURNAL, BASE_PATH 

694 JOURNAL = 'JSC_beamer' 

695 BASE_PATH = 'data/paper/time-x_website' 

696 

697 fig, ax = plt.subplots(figsize=figsize_by_journal(JOURNAL, 0.5, 2.0 / 3.0)) 

698 plot_recovery_rate_recoverable_only(get_stats(run_vdp), fig, ax) 

699 savefig(fig, 'recovery_rate', format='png') 

700 

701 from pySDC.projects.Resilience.work_precision import vdp_stiffness_plot 

702 

703 vdp_stiffness_plot(base_path=BASE_PATH, format='png') 

704 

705 

706def make_plots_for_SIAM_CSE23(): # pragma: no cover 

707 """ 

708 Make plots for the SIAM talk 

709 """