Coverage for pySDC/implementations/sweeper_classes/Multistep.py: 93%

91 statements  

« prev     ^ index     » next       coverage.py v7.10.6, created at 2025-09-04 15:08 +0000

1from pySDC.core.sweeper import Sweeper, _Pars 

2from pySDC.core.level import Level 

3 

4 

5class Cache(object): 

6 """ 

7 Class for managing solutions and right hand side evaluations of previous steps for the MultiStep "sweeper". 

8 

9 Attributes: 

10 - u (list): Contains solution from previous steps 

11 - f (list): Contains right hand side evaluations from previous steps 

12 - t (list): Contains time of previous steps 

13 """ 

14 

15 def __init__(self, num_steps): 

16 """ 

17 Initialization routing 

18 

19 Args: 

20 num_steps (int): Number of entries for the cache 

21 """ 

22 self.num_steps = num_steps 

23 self.u = [None] * num_steps 

24 self.f = [None] * num_steps 

25 self.t = [None] * num_steps 

26 

27 def update(self, t, u, f): 

28 """ 

29 Add a new value to the cache and remove the oldest one. 

30 

31 Args: 

32 t (float): Time of new step 

33 u (dtype_u): Solution of new step 

34 f (dtype_f): Right hand side evaluation at new step 

35 """ 

36 self.u[:-1] = self.u[1:] 

37 self.f[:-1] = self.f[1:] 

38 self.t[:-1] = self.t[1:] 

39 

40 self.u[-1] = u 

41 self.f[-1] = f 

42 self.t[-1] = t 

43 

44 def __str__(self): 

45 """ 

46 Print the contents of the cache for debugging purposes. 

47 """ 

48 string = '' 

49 for i in range(self.num_steps): 

50 string = f'{string} t={self.t[i]}: u={self.u[i]}, f={self.f[i]}' 

51 

52 return string 

53 

54 

55class MultiStep(Sweeper): 

56 alpha = None 

57 beta = None 

58 

59 def __init__(self, params, level): 

60 """ 

61 Initialization routine for the base sweeper. 

62 

63 Multistep methods work by constructing Euleresque steps with the right hand side composed of a weighted sum of solutions and right hand side evaluations at previous steps. 

64 The coefficients in the weighted sum for the solutions are called alpha and the ones for the right hand sides are called beta in this implementation. 

65 So for an N step method, there need to be N alpha values and N + 1 beta values, where the last one is on the left hand side for implicit methods. 

66 The first element in the coefficients belongs to the value furthest in the past and vice versa. Values from previous time steps are stored in a `Cache` object. 

67 Be careful with the sign of the alpha values. You can look at the implementations of the Euler methods for guidance. 

68 

69 Class attributes: 

70 alpha (list): Coefficients for the solutions of previous steps 

71 beta (list): Coefficients for the right hand side evaluations 

72 

73 Args: 

74 params (dict): parameter object 

75 level (pySDC.Level.level): the level that uses this sweeper 

76 """ 

77 import logging 

78 from pySDC.core.collocation import CollBase 

79 

80 # set up logger 

81 self.logger = logging.getLogger('sweeper') 

82 

83 self.params = _Pars(params) 

84 

85 # check if some parameters are set which only apply to actual sweepers 

86 for key in ['initial_guess', 'collocation_class', 'num_nodes', 'quad_type']: 

87 if key in params: 

88 self.logger.warning(f'"{key}" will be ignored by multistep sweeper') 

89 

90 # we need a dummy collocation object to instantiate the levels. 

91 self.coll = CollBase(num_nodes=1, quad_type='RADAU-RIGHT') 

92 

93 self.__level = level 

94 

95 self.parallelizable = False 

96 

97 # proprietary variables for the multistep methods 

98 self.steps = len(self.alpha) 

99 self.cache = Cache(self.steps) 

100 

101 @property 

102 def level(self): 

103 """ 

104 Returns the current level 

105 

106 Returns: 

107 pySDC.Level.level: Current level 

108 """ 

109 return self.__level 

110 

111 @level.setter 

112 def level(self, lvl): 

113 """ 

114 Sets a reference to the current level (done in the initialization of the level) 

115 

116 Args: 

117 lvl (pySDC.Level.level): Current level 

118 """ 

119 assert isinstance(lvl, Level), f"You tried to set the sweeper's level with an instance of {type(lvl)}!" 

120 

121 self.__level = lvl 

122 

123 def predict(self): 

124 """ 

125 Add the initial conditions to the cache if needed. 

126 

127 Default prediction for the sweepers, only copies the values to all collocation nodes 

128 and evaluates the RHS of the ODE there 

129 """ 

130 lvl = self.level 

131 

132 if all(me is None for me in self.cache.t): 

133 prob = lvl.prob 

134 lvl.f[0] = prob.eval_f(lvl.u[0], lvl.time) 

135 self.cache.update(lvl.time, lvl.u[0], lvl.f[0]) 

136 

137 lvl.status.unlocked = True 

138 lvl.status.updated = True 

139 

140 def compute_residual(self, stage=None): 

141 """ 

142 Do nothing. 

143 

144 Args: 

145 stage (str): The current stage of the step the level belongs to 

146 """ 

147 lvl = self.level 

148 lvl.status.residual = 0.0 

149 lvl.status.updated = False 

150 

151 return None 

152 

153 def compute_end_point(self): 

154 """ 

155 The solution is stored in the single node that we have. 

156 """ 

157 self.level.uend = self.level.u[-1] 

158 

159 def update_nodes(self): 

160 """ 

161 Compute the solution to the current step. If the cache is not filled, we compute a provisional solution with a different method. 

162 """ 

163 

164 lvl = self.level 

165 prob = lvl.prob 

166 time = lvl.time + lvl.dt 

167 

168 if None in self.cache.t: 

169 self.generate_starting_values() 

170 

171 else: 

172 # build the right hand side from the previous solutions 

173 rhs = prob.dtype_u(prob.init) 

174 dts = [self.cache.t[i + 1] - self.cache.t[i] for i in range(self.steps - 1)] + [ 

175 time - self.cache.t[-1] 

176 ] # TODO: Does this work for adaptive step sizes? 

177 for i in range(len(self.alpha)): 

178 rhs -= self.alpha[i] * self.cache.u[i] 

179 rhs += dts[i] * self.beta[i] * self.cache.f[i] 

180 

181 # compute the solution to the current step "implicit Euler style" 

182 lvl.u[1] = prob.solve_system(rhs, lvl.dt * self.beta[-1], self.cache.u[-1], time) 

183 

184 lvl.f[1] = prob.eval_f(lvl.u[1], time) 

185 self.cache.update(time, lvl.u[1], lvl.f[1]) 

186 

187 def generate_starting_values(self): 

188 """ 

189 Compute solutions to the steps when not enough previous values are available for the multistep method. 

190 The initial conditions are added in `predict` since this is not bespoke behaviour to any method. 

191 """ 

192 raise NotImplementedError( 

193 "No implementation for generating solutions when not enough previous values are available!" 

194 ) 

195 

196 

197class AdamsBashforthExplicit1Step(MultiStep): 

198 """ 

199 This is just forward Euler. 

200 """ 

201 

202 alpha = [-1.0] 

203 beta = [1.0, 0.0] 

204 

205 

206class BackwardEuler(MultiStep): 

207 """ 

208 Almost as old, impressive and beloved as Koelner Dom. 

209 """ 

210 

211 alpha = [-1.0] 

212 beta = [0.0, 1.0] 

213 

214 

215class AdamsMoultonImplicit1Step(MultiStep): 

216 """ 

217 Trapezoidal method dressed up as a multistep method. 

218 """ 

219 

220 alpha = [-1.0] 

221 beta = [0.5, 0.5] 

222 

223 

224class AdamsMoultonImplicit2Step(MultiStep): 

225 """ 

226 Third order implicit scheme 

227 """ 

228 

229 alpha = [0.0, -1.0] 

230 beta = [-1.0 / 12.0, 8.0 / 12.0, 5.0 / 12.0] 

231 

232 def generate_starting_values(self): 

233 """ 

234 Generate starting value by trapezoidal rule. 

235 """ 

236 lvl = self.level 

237 prob = lvl.prob 

238 time = lvl.time + lvl.dt 

239 

240 # do trapezoidal rule 

241 rhs = lvl.u[0] + lvl.dt / 2 * lvl.f[0] 

242 

243 lvl.u[1] = prob.solve_system(rhs, lvl.dt / 2.0, lvl.u[0], time)