Coverage for pySDC/tutorial/step_9/D_adaptive_alpha.py: 100%

1# --- 

2# jupyter: 

3# jupytext: 

4# formats: py:percent 

5# kernelspec: 

6# display_name: Python 3 

7# name: python3 

8# --- 

9 

10# %% [markdown] 

11# # Part D: Adaptive alpha 

12# 

13# Parts A to C all picked $\alpha$ by hand and kept it fixed, which means committing to one compromise for the whole 

14# run. A small $\alpha$ approximates the original problem better and converges in fewer iterations, but conditions 

15# the diagonalization worse, so round-off and inexact inner solves get amplified. The right balance shifts as the 

16# residual falls, so a fixed value is wrong at one end of the run or the other. 

17# 

18# The `AdaptiveAlpha` convergence controller updates $\alpha$ after every iteration instead, following 

19# [Čaklović et al.](https://doi.org/10.2140/camcos.2023.18.55): 

20# 

21# $$ 

22# \gamma = L (3 \epsilon + \tau), \quad 

23# \alpha_{k} = \sqrt{\frac{\gamma r_k}{e_k}}, \quad 

24# e_{k+1} = 2 \sqrt{\gamma e_k r_k}, 

25# $$ 

26# 

27# with $L$ the number of steps in the block, $\epsilon$ machine precision, $\tau$ the inner solver tolerance, $r_k$ 

28# the residual and $e_k$ a running bound on the error. 

29# 

30# ## The setup 

31# 

32# We compare a few fixed values of $\alpha$ with the adaptive one on the advection problem of 

33# [Part C](C_paradiag_in_pySDC), but with a direct solver for the complex shifted systems ParaDiag produces instead 

34# of GMRES, which struggles with them. The setup lives in `paradiag_setup.py` next to this tutorial, because 

35# [Part E](E_paradiag_MPI) uses exactly the same. Switching to the adaptive strategy is one entry in the description: 

36# 

37# :::{literalinclude} paradiag_setup.py 

38# :pyobject: get_controller_params 

39# ::: 

40# 

41# And `run` sets up the controller, the virtually parallel one when it gets no communicator: 

42# 

43# :::{literalinclude} paradiag_setup.py 

44# :pyobject: run 

45# ::: 

46 

47# %% 

48from pathlib import Path 

49 

50import matplotlib.pyplot as plt 

51import numpy as np 

52 

53from pySDC.tutorial.step_9.paradiag_setup import alpha_settings, format_result, num_steps_total, run 

54 

55# %% [markdown] 

56# ## The comparison 

57# 

58# One block holding every step, with the virtually parallel controller: 

59 

60# %% 

61# one block holding every step; Part E runs the same settings across MPI ranks 

62block_size = num_steps_total 

63 

64results = {} 

65lines = [] 

66for alpha in alpha_settings: 

67 uend, niter, error, final_alpha = run(alpha, block_size) 

68 results[alpha] = (uend, niter, final_alpha) 

69 lines.append(format_result('virtual', alpha, niter, error, final_alpha)) 

70 

71# Part E compares with these results, so they go into a file, too 

72Path("data").mkdir(parents=True, exist_ok=True) 

73with open('data/step_9_D_out.txt', 'w') as f: 

74 for line in lines: 

75 f.write(line + '\n') 

76 print(line) 

77 

78# %% tags=["hide-input"] 

79fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 3)) 

80labels = [str(alpha) for alpha in alpha_settings] 

81ax1.bar(labels, [results[alpha][1] for alpha in alpha_settings], color='#4f6bed') 

82ax1.set_title('iterations', fontsize=10) 

83ax2.bar(labels, [results[alpha][2] for alpha in alpha_settings], color='#e8743b') 

84ax2.set_yscale('log') 

85ax2.set_title(r'$\alpha$ at the end of the run', fontsize=10) 

86for ax in (ax1, ax2): 

87 ax.set_xlabel(r'$\alpha$ setting') 

88fig.tight_layout() 

89 

90# %% [markdown] 

91# :::{admonition} Important things to note 

92# - $\gamma$ is an accuracy floor. There is no point pushing $\alpha$ below the level at which round-off and the 

93# inner solver dominate anyway, which is why `inner_tol` enters: a looser inner solve should get a larger $\alpha$. 

94# - The interesting result is not that the adaptive strategy wins on iteration count. It ties with the best fixed 

95# value we tried, but it gets there without being told, and it keeps $\alpha$ orders of magnitude larger while doing 

96# so, which is exactly the margin that protects you once the inner solves are inexact. 

97# - The residual is taken over the whole block, so every rank computes the same $\alpha$ and the controllers stay 

98# in step. 

99# - $\alpha$ is a property of the method, not of the parallelization, so everything here runs with the virtually 

100# parallel controller. [Part E](E_paradiag_MPI) takes exactly these settings across MPI ranks and checks they come 

101# out the same. 

102# ::: 

103# 

104# The checks the tests run: 

105 

106# %% 

107# the adaptive strategy should need no more iterations than the best fixed alpha we tried 

108best_fixed = min(results[a][1] for a in alpha_settings if a != 'adaptive') 

109assert ( 

110 results['adaptive'][1] <= best_fixed 

111), f"ERROR: adaptive alpha needed {results['adaptive'][1]} iterations, the best fixed alpha only {best_fixed}" 

112 

113# alpha changes the iteration, not the problem, so all settings solve the same thing 

114reference = results[alpha_settings[0]][0] 

115for alpha in alpha_settings[1:]: 

116 assert np.allclose(results[alpha][0], reference, atol=1e-5), f'ERROR: alpha {alpha} gives a different solution'