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 

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(), AdaptivityPolynomialError()] 

49 if JOURNAL not in ['JSC_beamer']: 

50 strategies += [HotRodStrategy()] 

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 fig, axs = plt.subplots(1, 2, figsize=figsize_by_journal(JOURNAL, 1, 0.4), sharex=True) 

147 stats_analyser.plot_things_per_things( 

148 'recovered', 

149 'bit', 

150 False, 

151 op=stats_analyser.rec_rate, 

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

153 plotting_args={'markevery': 5}, 

154 ax=axs[0], 

155 ) 

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

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

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

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

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

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

162 

163 

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

165 """ 

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

167 

168 Args: 

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

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

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

172 

173 Returns: 

174 None 

175 """ 

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

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

178 

179 stats_analyser.plot_things_per_things( 

180 'recovered', 

181 'bit', 

182 False, 

183 op=stats_analyser.rec_rate, 

184 mask=fixable, 

185 args={**kwargs}, 

186 ax=ax, 

187 fig=fig, 

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

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

190 ) 

191 

192 

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

194 """ 

195 Compare the recovery rate for various problems. 

196 Only faults that can be recovered are shown. 

197 

198 Returns: 

199 None 

200 """ 

201 if target == 'resilience': 

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

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

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

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

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

207 else: 

208 raise NotImplementedError() 

209 

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

211 

212 my_setup_mpl() 

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

214 [ 

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

216 for i in range(len(stats)) 

217 ] 

218 

219 for ax in axs.flatten(): 

220 ax.get_legend().remove() 

221 

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

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

224 else: 

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

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

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

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

229 

230 if target == 'talk': 

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

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

233 

234 name = '' 

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

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

237 

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

239 

240 

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

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

243 stats.get_recovered() 

244 mask = None 

245 

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

247 if target == 'talk': 

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

249 else: 

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

251 

252 stats.plot_things_per_things( 

253 'recovered', 

254 x, 

255 False, 

256 op=stats.rec_rate, 

257 mask=mask, 

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

259 ax=ax, 

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

261 ) 

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

263 

264 if x == 'node': 

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

266 elif x == 'iteration': 

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

268 

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

270 

271 

272def plot_adaptivity_stuff(): # pragma: no cover 

273 """ 

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

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

276 

277 Returns: 

278 None 

279 """ 

280 from pySDC.implementations.hooks.log_errors import LogLocalErrorPostStep 

281 from pySDC.implementations.hooks.log_work import LogWork 

282 from pySDC.projects.Resilience.hook import LogData 

283 import pickle 

284 

285 my_setup_mpl() 

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

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

288 

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

290 """ 

291 Plot global error and cumulative sum of iterations 

292 

293 Args: 

294 stats (dict): Stats from pySDC run 

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

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

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

298 

299 Returns: 

300 None 

301 """ 

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

303 e = stats['e_local_post_step'] 

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

305 k = stats['work_newton'] 

306 iter_ax.plot( 

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

308 ) 

309 ax.set_yscale('log') 

310 ax.set_ylabel('local error') 

311 iter_ax.set_ylabel(r'Newton iterations') 

312 

313 run = False 

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

315 S = strategy(newton_inexactness=False) 

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

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

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

319 if strategy in [AdaptivityStrategy, BaseStrategy]: 

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

321 if strategy in [BaseStrategy, IterateStrategy]: 

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

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

324 if strategy in [IterateStrategy]: 

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

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

327 

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

329 if run: 

330 stats, _, _ = run_vdp( 

331 custom_description=desc, 

332 Tend=20, 

333 hook_class=[LogLocalErrorPostStep, LogWork, LogData], 

334 custom_controller_params={'logger_level': 15}, 

335 ) 

336 

337 data = { 

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

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

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

341 } 

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

343 pickle.dump(data, file) 

344 else: 

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

346 data = pickle.load(file) 

347 

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

349 

350 if strategy == BaseStrategy or True: 

351 u = data['u'] 

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

353 

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

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

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

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

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

359 savefig(fig, 'adaptivity') 

360 

361 

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

363 """ 

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

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

366 

367 Args: 

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

369 

370 Returns: 

371 None 

372 """ 

373 from pySDC.projects.Resilience.fault_stats import ( 

374 FaultStats, 

375 BaseStrategy, 

376 ) 

377 from pySDC.projects.Resilience.hook import LogData 

378 

379 stats_analyser = FaultStats( 

380 prob=run_vdp, 

381 strategies=[BaseStrategy()], 

382 faults=[False, True], 

383 reload=True, 

384 recovery_thresh=1.1, 

385 num_procs=1, 

386 mode='combination', 

387 ) 

388 

389 my_setup_mpl() 

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

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

392 ls = ['--', '-'] 

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

394 do_faults = [False, True] 

395 superscripts = ['*', ''] 

396 subscripts = ['', 't', ''] 

397 

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

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

400 

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

402 stats, controller, Tend = stats_analyser.single_run( 

403 strategy=BaseStrategy(), 

404 run=run, 

405 faults=do_faults[i], 

406 hook_class=[LogData], 

407 ) 

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

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

410 for j in [0, 1]: 

411 ax.plot( 

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

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

414 ls=ls[i], 

415 color=colors[j], 

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

417 marker=markers[j], 

418 markevery=60, 

419 ) 

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

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

422 print( 

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

424 ) 

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

426 

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

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

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

430 

431 

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

433 """ 

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

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

436 

437 Args: 

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

439 

440 Returns: 

441 None 

442 """ 

443 from pySDC.projects.Resilience.fault_stats import ( 

444 FaultStats, 

445 BaseStrategy, 

446 ) 

447 from pySDC.projects.Resilience.hook import LogData 

448 

449 stats_analyser = FaultStats( 

450 prob=run_Lorenz, 

451 strategies=[BaseStrategy()], 

452 faults=[False, True], 

453 reload=True, 

454 recovery_thresh=1.1, 

455 num_procs=1, 

456 mode='combination', 

457 ) 

458 

459 strategy = BaseStrategy() 

460 

461 my_setup_mpl() 

462 if target == 'resilience': 

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

464 else: 

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

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

467 ls = ['--', '-'] 

468 markers = [None, strategy.marker] 

469 do_faults = [False, True] 

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

471 labels = ['x', 'x'] 

472 

473 run = 19 + 20 * bit 

474 

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

476 stats, controller, Tend = stats_analyser.single_run( 

477 strategy=BaseStrategy(), 

478 run=run, 

479 faults=do_faults[i], 

480 hook_class=[LogData], 

481 ) 

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

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

484 ax.plot( 

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

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

487 ls=ls[i], 

488 color=colors[i], 

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

490 marker=markers[i], 

491 markevery=500, 

492 ) 

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

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

495 print( 

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

497 ) 

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

499 

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

501 

502 h, l = ax.get_legend_handles_labels() 

503 fig.legend( 

504 h, 

505 l, 

506 loc='outside lower center', 

507 ncols=3, 

508 frameon=False, 

509 fancybox=True, 

510 borderaxespad=0.01, 

511 ) 

512 

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

514 

515 

516def plot_Lorenz_solution(): # pragma: no cover 

517 my_setup_mpl() 

518 

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

520 

521 strategy = BaseStrategy() 

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

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

524 

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

526 

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

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

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

530 

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

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

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

534 

535 for ax in axs: 

536 ax.set_box_aspect(1.0) 

537 

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

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

540 

541 

542def plot_quench_solution(): # pragma: no cover 

543 """ 

544 Plot the solution of Quench problem over time 

545 

546 Returns: 

547 None 

548 """ 

549 my_setup_mpl() 

550 if JOURNAL == 'JSC_beamer': 

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

552 else: 

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

554 

555 strategy = BaseStrategy() 

556 

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

558 

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

560 

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

562 

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

564 

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

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

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

568 

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

570 ax.legend(frameon=False) 

571 savefig(fig, 'quench_sol') 

572 

573 

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

575 """ 

576 Plot solution of Rayleigh-Benard convection 

577 """ 

578 my_setup_mpl() 

579 

580 from mpl_toolkits.axes_grid1 import make_axes_locatable 

581 

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

583 aspect = 0.8 if nplots == 3 else 0.5 

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

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

586 caxs = [] 

587 for ax in axs: 

588 divider = make_axes_locatable(ax) 

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

590 

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

592 

593 prob = RayleighBenard(**PROBLEM_PARAMS) 

594 

595 def _plot(t, ax, cax): 

596 u_hat = prob.u_exact(t) 

597 u = prob.itransform(u_hat) 

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

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

600 

601 if setup == 'resilience': 

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

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

604 elif setup == 'work-precision': 

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

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

607 elif setup == 'resilience_thesis': 

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

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

610 elif setup == 'thesis_intro': 

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

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

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

614 

615 for ax in axs: 

616 ax.set_ylabel('$z$') 

617 ax.set_aspect(1) 

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

619 

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

621 

622 

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

624 my_setup_mpl() 

625 

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

627 

628 from mpl_toolkits.axes_grid1 import make_axes_locatable 

629 

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

631 cax = [] 

632 divider = make_axes_locatable(axs[0]) 

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

634 divider2 = make_axes_locatable(axs[1]) 

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

636 

637 from pySDC.projects.Resilience.GS import grayscott_imex_diffusion 

638 

639 problem_params = { 

640 'num_blobs': -48, 

641 'L': 2, 

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

643 'A': 0.062, 

644 'B': 0.1229, 

645 'Du': 2e-5, 

646 'Dv': 1e-5, 

647 } 

648 P = grayscott_imex_diffusion(**problem_params) 

649 Tend = tend 

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

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

652 

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

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

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

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

657 for ax in axs: 

658 ax.set_aspect(1) 

659 ax.set_xlabel('$x$') 

660 ax.set_ylabel('$y$') 

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

662 

663 

664def plot_Schroedinger_solution(): # pragma: no cover 

665 from pySDC.implementations.problem_classes.NonlinearSchroedinger_MPIFFT import nonlinearschroedinger_imex 

666 

667 my_setup_mpl() 

668 if JOURNAL == 'JSC_beamer': 

669 raise NotImplementedError 

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

671 else: 

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

673 

674 from mpl_toolkits.axes_grid1 import make_axes_locatable 

675 

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

677 cax = [] 

678 divider = make_axes_locatable(axs[0]) 

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

680 divider2 = make_axes_locatable(axs[1]) 

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

682 

683 problem_params = dict() 

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

685 problem_params['spectral'] = False 

686 problem_params['c'] = 1.0 

687 description = {'problem_params': problem_params} 

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

689 

690 P = nonlinearschroedinger_imex(**problem_params) 

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

692 

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

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

695 

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

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

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

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

700 for ax in axs: 

701 ax.set_aspect(1) 

702 ax.set_xlabel('$x$') 

703 ax.set_ylabel('$y$')