Coverage for pySDC/implementations/convergence_controller_classes/adaptive_alpha.py: 100%

1import numpy as np 

2 

3from pySDC.core.convergence_controller import ConvergenceController 

4 

5 

6class AdaptiveAlpha(ConvergenceController): 

7 r""" 

8 Choose the ParaDiag :math:`\alpha` adaptively from the residual. 

9 

10 ParaDiag replaces the time-stepping matrix by an :math:`\alpha`-circulant approximation, and 

11 :math:`\alpha` trades two error sources against each other: a small value approximates the original 

12 problem better, but conditions the diagonalization worse, so round-off and inexact inner solves 

13 contaminate the result. A single fixed :math:`\alpha` therefore has to be a compromise for the whole 

14 run, even though the balance shifts as the residual falls. 

15 

16 This convergence controller updates :math:`\alpha` after every iteration instead, following the 

17 strategy in `Čaklović et al. <https://doi.org/10.2140/camcos.2023.18.55>`_: 

18 

19 .. math:: 

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

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

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

23 

24 with :math:`L` the number of steps in the block, :math:`\epsilon` machine precision, :math:`\tau` 

25 the inner solver tolerance, :math:`r_k` the residual and :math:`e_k` a running bound on the error. 

26 :math:`\gamma` is the accuracy floor: there is no point pushing :math:`\alpha` below the level at 

27 which round-off and the inner solver dominate anyway. 

28 

29 The residual is reduced over the whole block, so every rank computes the same :math:`\alpha` and the 

30 controllers stay in step. 

31 """ 

32 

33 def setup(self, controller, params, description, **kwargs): 

34 """ 

35 Define default parameters here. 

36 

37 Args: 

38 controller (pySDC.Controller): The controller 

39 params (dict): The params passed for this specific convergence controller 

40 description (dict): The description object used to instantiate the controller 

41 

42 Returns: 

43 (dict): The updated params dictionary 

44 """ 

45 defaults = { 

46 'control_order': +300, 

47 # accuracy floor: round-off plus whatever the inner solver leaves behind 

48 'inner_tol': 0.0, 

49 # initial bound on the error, before we have seen a residual 

50 'e0': 1.0, 

51 # keep alpha in a sane range no matter what the residual does 

52 'alpha_min': 1e-12, 

53 'alpha_max': 1.0, 

54 } 

55 return {**defaults, **super().setup(controller, params, description, **kwargs)} 

56 

57 def setup_status_variables(self, controller, **kwargs): 

58 """ 

59 Start the alpha history, which spans the whole run rather than a single block. 

60 

61 Args: 

62 controller (pySDC.Controller): The controller 

63 

64 Returns: 

65 None 

66 """ 

67 self.alphas = [] 

68 return None 

69 

70 def reset_status_variables(self, controller, **kwargs): 

71 """ 

72 Reset the error bound at the start of every block. 

73 

74 Args: 

75 controller (pySDC.Controller): The controller 

76 

77 Returns: 

78 None 

79 """ 

80 self.e = self.params.e0 

81 return None 

82 

83 def get_gamma(self, controller): 

84 r""" 

85 The accuracy floor :math:`\gamma = L (3 \epsilon + \tau)`. 

86 

87 Args: 

88 controller (pySDC.Controller): The controller 

89 

90 Returns: 

91 float: gamma 

92 """ 

93 eps = np.finfo(complex).eps 

94 return controller.n_steps * (3 * eps + self.params.inner_tol) 

95 

96 def post_iteration_processing_block(self, controller, **kwargs): 

97 r""" 

98 Compute the next :math:`\alpha` from the residual of the whole block. 

99 

100 This is a block hook rather than a per-step one because :math:`\alpha` belongs to the block: 

101 it parametrises the transform across the steps, so it has to advance once per iteration no 

102 matter how the block was decomposed. 

103 

104 Args: 

105 controller (pySDC.Controller): The controller 

106 

107 Returns: 

108 None 

109 """ 

110 # the residual of the composite problem is the largest one across the block 

111 residual = max(step.levels[0].status.residual for step in controller.steps) 

112 

113 comm = kwargs.get('comm', None) 

114 if comm is not None: 

115 residual = comm.allreduce(residual, op=self.MPI_MAX) 

116 

117 if residual <= 0: 

118 return None 

119 

120 gamma = self.get_gamma(controller) 

121 alpha = np.sqrt(gamma * residual / self.e) 

122 self.e = 2 * np.sqrt(gamma * self.e * residual) 

123 

124 controller.params.alpha = min(max(alpha, self.params.alpha_min), self.params.alpha_max) 

125 self.debug(f'Set alpha to {controller.params.alpha:.3e} from residual {residual:.3e}', controller.steps[0]) 

126 

127 # keep the history around; it is what you want to look at when tuning 

128 self.alphas.append(controller.params.alpha) 

129 

130 return None 

131 

132 def dependencies(self, controller, description, **kwargs): 

133 """ 

134 Prepare the MPI reduction we need for the block residual. 

135 

136 Args: 

137 controller (pySDC.Controller): The controller 

138 description (dict): The description object used to instantiate the controller 

139 

140 Returns: 

141 None 

142 """ 

143 if self.params.useMPI: 

144 self.prepare_MPI_datatypes() 

145 from mpi4py import MPI 

146 

147 self.MPI_MAX = MPI.MAX 

148 

149 super().dependencies(controller, description, **kwargs) 

150 return None