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 run_GS, 

14 RECOVERY_THRESH_ABS, 

15) 

16from pySDC.projects.Resilience.strategies import ( 

17 BaseStrategy, 

18 AdaptivityStrategy, 

19 IterateStrategy, 

20 HotRodStrategy, 

21 DIRKStrategy, 

22 ERKStrategy, 

23 AdaptivityPolynomialError, 

24 cmap, 

25) 

26from pySDC.helpers.plot_helper import setup_mpl, figsize_by_journal 

27from pySDC.helpers.stats_helper import get_sorted 

28 

29 

30cm = 1 / 2.5 

31TEXTWIDTH = 11.9446244611 * cm 

32JOURNAL = 'Springer_Numerical_Algorithms' 

33BASE_PATH = 'data/paper' 

34 

35 

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

37 """ 

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

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

40 

41 Args: 

42 problem (function): A problem to run 

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

44 

45 Returns: 

46 FaultStats: Object to analyse resilience statistics from 

47 """ 

48 if strategy_type == 'SDC': 

49 strategies = [BaseStrategy(), AdaptivityStrategy(), IterateStrategy(), AdaptivityPolynomialError()] 

50 if JOURNAL not in ['JSC_beamer']: 

51 strategies += [HotRodStrategy()] 

52 elif strategy_type == 'RK': 

53 strategies = [DIRKStrategy()] 

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

55 strategies += [ERKStrategy()] 

56 

57 stats_analyser = FaultStats( 

58 prob=problem, 

59 strategies=strategies, 

60 faults=[False, True], 

61 reload=True, 

62 recovery_thresh=1.1, 

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

64 mode='default', 

65 stats_path=path, 

66 num_procs=num_procs, 

67 ) 

68 stats_analyser.get_recovered() 

69 return stats_analyser 

70 

71 

72def my_setup_mpl(**kwargs): 

73 setup_mpl(reset=True, font_size=8) 

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

75 

76 

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

78 """ 

79 Save a figure to some predefined location. 

80 

81 Args: 

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

83 name (str): The name of the plot 

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

85 Returns: 

86 None 

87 """ 

88 if tight_layout: 

89 fig.tight_layout() 

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

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

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

93 

94 

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

96 """ 

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

98 

99 Args: 

100 problem (function): A problem to run 

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

102 

103 Returns: 

104 None 

105 """ 

106 

107 stats_analyser = get_stats(problem, path) 

108 stats_analyser.get_recovered() 

109 

110 strategy = IterateStrategy() 

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

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

113 stats_analyser.print_faults(not_overflow) 

114 

115 compare_strategies(stats_analyser, **kwargs) 

116 plot_recovery_rate(stats_analyser, **kwargs) 

117 

118 

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

120 """ 

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

122 

123 Args: 

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

125 

126 Returns: 

127 None 

128 """ 

129 my_setup_mpl() 

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

131 stats_analyser.compare_strategies(ax=ax) 

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

133 

134 

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

136 """ 

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

138 

139 Args: 

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

141 

142 Returns: 

143 None 

144 """ 

145 my_setup_mpl() 

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

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

148 stats_analyser.plot_things_per_things( 

149 'recovered', 

150 'bit', 

151 False, 

152 op=stats_analyser.rec_rate, 

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

154 plotting_args={'markevery': 5}, 

155 ax=axs[0], 

156 ) 

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

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

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

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

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

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

163 

164 

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

166 """ 

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

168 

169 Args: 

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

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

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

173 

174 Returns: 

175 None 

176 """ 

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

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

179 

180 stats_analyser.plot_things_per_things( 

181 'recovered', 

182 'bit', 

183 False, 

184 op=stats_analyser.rec_rate, 

185 mask=fixable, 

186 args={**kwargs}, 

187 ax=ax, 

188 fig=fig, 

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

190 plotting_args={'markevery': 10 if stats_analyser.prob.__name__ in ['run_RBC', 'run_Schroedinger'] else 5}, 

191 ) 

192 

193 

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

195 """ 

196 Compare the recovery rate for various problems. 

197 Only faults that can be recovered are shown. 

198 

199 Returns: 

200 None 

201 """ 

202 if target == 'resilience': 

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

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

205 elif target in ['thesis', 'talk']: 

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

207 titles = ['Van der Pol', 'Lorenz', 'Gray-Scott', 'Rayleigh-Benard'] 

208 else: 

209 raise NotImplementedError() 

210 

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

212 

213 my_setup_mpl() 

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

215 [ 

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

217 for i in range(len(stats)) 

218 ] 

219 

220 for ax in axs.flatten(): 

221 ax.get_legend().remove() 

222 

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

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

225 else: 

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

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

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

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

230 

231 if target == 'talk': 

232 axs[0, 0].set_xlabel('') 

233 axs[0, 1].set_xlabel('') 

234 

235 name = '' 

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

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

238 

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

240 

241 

242def plot_recovery_rate_detailed_Lorenz(target='resilience'): # pragma: no cover 

243 stats = get_stats(run_Lorenz, num_procs=1, strategy_type='SDC') 

244 stats.get_recovered() 

245 mask = None 

246 

247 for x in ['node', 'iteration', 'bit']: 

248 if target == 'talk': 

249 fig, ax = plt.subplots(figsize=figsize_by_journal(JOURNAL, 0.6, 0.6)) 

250 else: 

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

252 

253 stats.plot_things_per_things( 

254 'recovered', 

255 x, 

256 False, 

257 op=stats.rec_rate, 

258 mask=mask, 

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

260 ax=ax, 

261 plotting_args={'markevery': 5 if x == 'bit' else 1}, 

262 ) 

263 ax.set_ylim((-0.05, 1.05)) 

264 

265 if x == 'node': 

266 ax.set_xticks([0, 1, 2, 3]) 

267 elif x == 'iteration': 

268 ax.set_xticks([1, 2, 3, 4, 5]) 

269 

270 savefig(fig, f'recovery_rate_Lorenz_{x}') 

271 

272 

273def plot_adaptivity_stuff(): # pragma: no cover 

274 """ 

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

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

277 

278 Returns: 

279 None 

280 """ 

281 from pySDC.implementations.hooks.log_errors import LogLocalErrorPostStep 

282 from pySDC.implementations.hooks.log_work import LogWork 

283 from pySDC.projects.Resilience.hook import LogData 

284 import pickle 

285 

286 my_setup_mpl() 

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

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

289 

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

291 """ 

292 Plot global error and cumulative sum of iterations 

293 

294 Args: 

295 stats (dict): Stats from pySDC run 

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

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

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

299 

300 Returns: 

301 None 

302 """ 

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

304 e = stats['e_local_post_step'] 

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

306 k = stats['work_newton'] 

307 iter_ax.plot( 

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

309 ) 

310 ax.set_yscale('log') 

311 ax.set_ylabel('local error') 

312 iter_ax.set_ylabel(r'Newton iterations') 

313 

314 run = False 

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

316 S = strategy(newton_inexactness=False) 

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

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

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

320 if strategy in [AdaptivityStrategy, BaseStrategy]: 

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

322 if strategy in [BaseStrategy, IterateStrategy]: 

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

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

325 if strategy in [IterateStrategy]: 

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

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

328 

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

330 if run: 

331 stats, _, _ = run_vdp( 

332 custom_description=desc, 

333 Tend=20, 

334 hook_class=[LogLocalErrorPostStep, LogWork, LogData], 

335 custom_controller_params={'logger_level': 15}, 

336 ) 

337 

338 data = { 

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

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

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

342 } 

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

344 pickle.dump(data, file) 

345 else: 

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

347 data = pickle.load(file) 

348 

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

350 

351 if strategy == BaseStrategy or True: 

352 u = data['u'] 

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

354 

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

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

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

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

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

360 savefig(fig, 'adaptivity') 

361 

362 

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

364 """ 

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

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

367 

368 Args: 

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

370 

371 Returns: 

372 None 

373 """ 

374 from pySDC.projects.Resilience.fault_stats import ( 

375 FaultStats, 

376 BaseStrategy, 

377 ) 

378 from pySDC.projects.Resilience.hook import LogData 

379 

380 stats_analyser = FaultStats( 

381 prob=run_vdp, 

382 strategies=[BaseStrategy()], 

383 faults=[False, True], 

384 reload=True, 

385 recovery_thresh=1.1, 

386 num_procs=1, 

387 mode='combination', 

388 ) 

389 

390 my_setup_mpl() 

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

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

393 ls = ['--', '-'] 

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

395 do_faults = [False, True] 

396 superscripts = ['*', ''] 

397 subscripts = ['', 't', ''] 

398 

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

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

401 

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

403 stats, controller, Tend = stats_analyser.single_run( 

404 strategy=BaseStrategy(), 

405 run=run, 

406 faults=do_faults[i], 

407 hook_class=[LogData], 

408 ) 

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

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

411 for j in [0, 1]: 

412 ax.plot( 

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

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

415 ls=ls[i], 

416 color=colors[j], 

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

418 marker=markers[j], 

419 markevery=60, 

420 ) 

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

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

423 print( 

424 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]}' 

425 ) 

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

427 

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

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

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

431 

432 

433def plot_fault_Lorenz(bit=0, target='resilience'): # pragma: no cover 

434 """ 

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

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

437 

438 Args: 

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

440 

441 Returns: 

442 None 

443 """ 

444 from pySDC.projects.Resilience.fault_stats import ( 

445 FaultStats, 

446 BaseStrategy, 

447 ) 

448 from pySDC.projects.Resilience.hook import LogData 

449 

450 stats_analyser = FaultStats( 

451 prob=run_Lorenz, 

452 strategies=[BaseStrategy()], 

453 faults=[False, True], 

454 reload=True, 

455 recovery_thresh=1.1, 

456 num_procs=1, 

457 mode='combination', 

458 ) 

459 

460 strategy = BaseStrategy() 

461 

462 my_setup_mpl() 

463 if target == 'resilience': 

464 fig, ax = plt.subplots(figsize=figsize_by_journal(JOURNAL, 0.4, 0.6)) 

465 else: 

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

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

468 ls = ['--', '-'] 

469 markers = [None, strategy.marker] 

470 do_faults = [False, True] 

471 superscripts = [r'\mathrm{no~faults}', ''] 

472 labels = ['x', 'x'] 

473 

474 run = 19 + 20 * bit 

475 

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

477 stats, controller, Tend = stats_analyser.single_run( 

478 strategy=BaseStrategy(), 

479 run=run, 

480 faults=do_faults[i], 

481 hook_class=[LogData], 

482 ) 

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

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

485 ax.plot( 

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

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

488 ls=ls[i], 

489 color=colors[i], 

490 label=rf'${ {labels[i]}} _{ {superscripts[i]}} $', 

491 marker=markers[i], 

492 markevery=500, 

493 ) 

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

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

496 print( 

497 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]}' 

498 ) 

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

500 

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

502 

503 h, l = ax.get_legend_handles_labels() 

504 fig.legend( 

505 h, 

506 l, 

507 loc='outside lower center', 

508 ncols=3, 

509 frameon=False, 

510 fancybox=True, 

511 borderaxespad=0.01, 

512 ) 

513 

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

515 

516 

517def plot_Lorenz_solution(): # pragma: no cover 

518 my_setup_mpl() 

519 

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

521 

522 strategy = BaseStrategy() 

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

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

525 

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

527 

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

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

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

531 

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

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

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

535 

536 for ax in axs: 

537 ax.set_box_aspect(1.0) 

538 

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

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

541 

542 

543def plot_quench_solution(): # pragma: no cover 

544 """ 

545 Plot the solution of Quench problem over time 

546 

547 Returns: 

548 None 

549 """ 

550 my_setup_mpl() 

551 if JOURNAL == 'JSC_beamer': 

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

553 else: 

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

555 

556 strategy = BaseStrategy() 

557 

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

559 

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

561 

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

563 

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

565 

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

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

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

569 

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

571 ax.legend(frameon=False) 

572 savefig(fig, 'quench_sol') 

573 

574 

575def plot_RBC_solution(setup='resilience'): # pragma: no cover 

576 """ 

577 Plot solution of Rayleigh-Benard convection 

578 """ 

579 my_setup_mpl() 

580 

581 from mpl_toolkits.axes_grid1 import make_axes_locatable 

582 

583 nplots = 3 if setup == 'thesis_intro' else 2 

584 aspect = 0.8 if nplots == 3 else 0.5 

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

586 fig, axs = plt.subplots(nplots, 1, sharex=True, sharey=True, figsize=figsize_by_journal(JOURNAL, 1.0, aspect)) 

587 caxs = [] 

588 for ax in axs: 

589 divider = make_axes_locatable(ax) 

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

591 

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

593 

594 prob = RayleighBenard(**PROBLEM_PARAMS) 

595 

596 def _plot(t, ax, cax): 

597 u_hat = prob.u_exact(t) 

598 u = prob.itransform(u_hat) 

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

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

601 

602 if setup == 'resilience': 

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

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

605 elif setup == 'work-precision': 

606 _plot(10, axs[0], caxs[0]) 

607 _plot(16, axs[1], caxs[1]) 

608 elif setup == 'resilience_thesis': 

609 _plot(20, axs[0], caxs[0]) 

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

611 elif setup == 'thesis_intro': 

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

613 _plot(18, axs[1], caxs[1]) 

614 _plot(30, axs[2], caxs[2]) 

615 

616 for ax in axs: 

617 ax.set_ylabel('$z$') 

618 ax.set_aspect(1) 

619 axs[-1].set_xlabel('$x$') 

620 

621 savefig(fig, f'RBC_sol_{setup}', tight_layout=False) 

622 

623 

624def plot_GS_solution(tend=500): # pragma: no cover 

625 my_setup_mpl() 

626 

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

628 

629 from mpl_toolkits.axes_grid1 import make_axes_locatable 

630 

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

632 cax = [] 

633 divider = make_axes_locatable(axs[0]) 

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

635 divider2 = make_axes_locatable(axs[1]) 

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

637 

638 from pySDC.projects.Resilience.GS import grayscott_imex_diffusion 

639 

640 problem_params = { 

641 'num_blobs': -48, 

642 'L': 2, 

643 'nvars': (128,) * 2, 

644 'A': 0.062, 

645 'B': 0.1229, 

646 'Du': 2e-5, 

647 'Dv': 1e-5, 

648 } 

649 P = grayscott_imex_diffusion(**problem_params) 

650 Tend = tend 

651 im = axs[0].pcolormesh(*P.X, P.u_exact(0)[1], rasterized=True, cmap='binary') 

652 im1 = axs[1].pcolormesh(*P.X, P.u_exact(Tend)[1], rasterized=True, cmap='binary') 

653 

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

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

656 axs[0].set_title(r'$v(t=0)$') 

657 axs[1].set_title(rf'$v(t={ {Tend}} )$') 

658 for ax in axs: 

659 ax.set_aspect(1) 

660 ax.set_xlabel('$x$') 

661 ax.set_ylabel('$y$') 

662 savefig(fig, f'GrayScott_sol{f"_{tend}" if tend != 500 else ""}') 

663 

664 

665def plot_Schroedinger_solution(): # pragma: no cover 

666 from pySDC.implementations.problem_classes.NonlinearSchroedinger_MPIFFT import nonlinearschroedinger_imex 

667 

668 my_setup_mpl() 

669 if JOURNAL == 'JSC_beamer': 

670 raise NotImplementedError 

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

672 else: 

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

674 

675 from mpl_toolkits.axes_grid1 import make_axes_locatable 

676 

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

678 cax = [] 

679 divider = make_axes_locatable(axs[0]) 

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

681 divider2 = make_axes_locatable(axs[1]) 

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

683 

684 problem_params = dict() 

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

686 problem_params['spectral'] = False 

687 problem_params['c'] = 1.0 

688 description = {'problem_params': problem_params} 

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

690 

691 P = nonlinearschroedinger_imex(**problem_params) 

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

693 

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

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

696 

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

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

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

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

701 for ax in axs: 

702 ax.set_aspect(1) 

703 ax.set_xlabel('$x$') 

704 ax.set_ylabel('$y$')