Coverage for pySDC/implementations/convergence_controller_classes/estimate_embedded_error.py: 77%

113 statements  

« prev     ^ index     » next       coverage.py v7.8.0, created at 2025-04-01 13:12 +0000

1import numpy as np 

2 

3from pySDC.core.convergence_controller import ConvergenceController, Pars, Status 

4from pySDC.implementations.convergence_controller_classes.store_uold import StoreUOld 

5 

6from pySDC.implementations.sweeper_classes.Runge_Kutta import RungeKutta 

7 

8 

9class EstimateEmbeddedError(ConvergenceController): 

10 """ 

11 The embedded error is obtained by computing two solutions of different accuracy and pretending the more accurate 

12 one is an exact solution from the point of view of the less accurate solution. In practice, we like to compute the 

13 solutions with different order methods, meaning that in SDC we can just subtract two consecutive sweeps, as long as 

14 you make sure your preconditioner is compatible, which you have to just try out... 

15 """ 

16 

17 @classmethod 

18 def get_implementation(cls, flavor='standard', useMPI=False): 

19 """ 

20 Retrieve the implementation for a specific flavor of this class. 

21 

22 Args: 

23 flavor (str): The implementation that you want 

24 

25 Returns: 

26 cls: The child class that implements the desired flavor 

27 """ 

28 if flavor == 'standard': 

29 return cls 

30 elif flavor == 'linearized': 

31 if useMPI: 

32 return EstimateEmbeddedErrorLinearizedMPI 

33 else: 

34 return EstimateEmbeddedErrorLinearizedNonMPI 

35 elif flavor == 'collocation': 

36 return EstimateEmbeddedErrorCollocation 

37 else: 

38 raise NotImplementedError(f'Flavor {flavor} of EstimateEmbeddedError is not implemented!') 

39 

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

41 """ 

42 Add a default value for control order to the parameters and check if we are using a Runge-Kutta sweeper 

43 

44 Args: 

45 controller (pySDC.Controller): The controller 

46 params (dict): Parameters for the convergence controller 

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

48 

49 Returns: 

50 dict: Updated parameters 

51 """ 

52 sweeper_type = 'SDC' 

53 if RungeKutta in description['sweeper_class'].__mro__: 

54 sweeper_type = 'RK' 

55 elif 'SweeperMPI' in [me.__name__ for me in description['sweeper_class'].__mro__]: 

56 sweeper_type = 'MPI' 

57 return { 

58 "control_order": -80, 

59 "sweeper_type": sweeper_type, 

60 "rel_error": False, 

61 **super().setup(controller, params, description, **kwargs), 

62 } 

63 

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

65 """ 

66 Load the convergence controller that stores the solution of the last sweep unless we are doing Runge-Kutta. 

67 Add the hook for recording the error. 

68 

69 Args: 

70 controller (pySDC.Controller): The controller 

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

72 

73 Returns: 

74 None 

75 """ 

76 if RungeKutta not in description["sweeper_class"].__bases__: 

77 controller.add_convergence_controller(StoreUOld, description=description) 

78 

79 from pySDC.implementations.hooks.log_embedded_error_estimate import LogEmbeddedErrorEstimate 

80 

81 controller.add_hook(LogEmbeddedErrorEstimate) 

82 return None 

83 

84 def estimate_embedded_error_serial(self, L): 

85 """ 

86 Estimate the serial embedded error, which may need to be modified for a parallel estimate. 

87 

88 Depending on the type of sweeper, the lower order solution is stored in a different place. 

89 

90 Args: 

91 L (pySDC.level): The level 

92 

93 Returns: 

94 dtype_u: The embedded error estimate 

95 """ 

96 if self.params.sweeper_type == "RK": 

97 L.sweep.compute_end_point() 

98 if self.params.rel_error: 

99 return abs(L.uend - L.sweep.u_secondary) / abs(L.uend) 

100 else: 

101 return abs(L.uend - L.sweep.u_secondary) 

102 elif self.params.sweeper_type == "SDC": 

103 # order rises by one between sweeps 

104 if self.params.rel_error: 

105 return abs(L.uold[-1] - L.u[-1]) / abs(L.u[-1]) 

106 else: 

107 return abs(L.uold[-1] - L.u[-1]) 

108 elif self.params.sweeper_type == 'MPI': 

109 comm = L.sweep.comm 

110 if self.params.rel_error: 

111 return comm.bcast( 

112 abs(L.uold[comm.rank + 1] - L.u[comm.rank + 1]) / abs(L.u[comm.rank + 1]), root=comm.size - 1 

113 ) 

114 else: 

115 return comm.bcast(abs(L.uold[comm.rank + 1] - L.u[comm.rank + 1]), root=comm.size - 1) 

116 else: 

117 raise NotImplementedError( 

118 f"Don't know how to estimate embedded error for sweeper type \ 

119\"{self.params.sweeper_type}\"" 

120 ) 

121 

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

123 """ 

124 Add the embedded error to the level status 

125 

126 Args: 

127 controller (pySDC.Controller): The controller 

128 """ 

129 self.add_status_variable_to_level('error_embedded_estimate') 

130 self.add_status_variable_to_level('increment') 

131 

132 def post_iteration_processing(self, controller, S, **kwargs): 

133 """ 

134 Estimate the local error here. 

135 

136 If you are doing MSSDC, this is the global error within the block in Gauss-Seidel mode. 

137 In Jacobi mode, I haven't thought about what this is. 

138 

139 Args: 

140 controller (pySDC.Controller): The controller 

141 S (pySDC.Step): The current step 

142 

143 Returns: 

144 None 

145 """ 

146 

147 if S.status.iter > 0 or self.params.sweeper_type == "RK": 

148 for L in S.levels: 

149 L.status.error_embedded_estimate = max([self.estimate_embedded_error_serial(L), np.finfo(float).eps]) 

150 L.status.increment = L.status.error_embedded_estimate * 1 

151 self.debug(f'L.status.error_embedded_estimate={L.status.error_embedded_estimate:.5e}', S) 

152 

153 return None 

154 

155 

156class EstimateEmbeddedErrorLinearizedNonMPI(EstimateEmbeddedError): 

157 def __init__(self, controller, params, description, **kwargs): 

158 """ 

159 Initialisation routine. Add the buffers for communication. 

160 

161 Args: 

162 controller (pySDC.Controller): The controller 

163 params (dict): Parameters for the convergence controller 

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

165 """ 

166 super().__init__(controller, params, description, **kwargs) 

167 self.buffers = Pars({'e_em_last': 0.0}) 

168 

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

170 """ 

171 Add option for averaging the local errors. 

172 

173 Args: 

174 controller (pySDC.Controller): The controller 

175 params (dict): Parameters for the convergence controller 

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

177 

178 Returns: 

179 dict: Updated parameters 

180 """ 

181 return { 

182 'averaged': False, 

183 **super().setup(controller, params, description, **kwargs), 

184 } 

185 

186 def reset_buffers_nonMPI(self, controller, **kwargs): 

187 """ 

188 Reset buffers for imitated communication. 

189 

190 Args: 

191 controller (pySDC.controller): The controller 

192 

193 Returns: 

194 None 

195 """ 

196 self.buffers.e_em_last = 0.0 

197 return None 

198 

199 def post_iteration_processing(self, controller, S, **kwargs): 

200 """ 

201 Compute embedded error estimate on the last node of each level 

202 In serial this is the local error, but in block Gauss-Seidel MSSDC this is a semi-global error in each block 

203 

204 Args: 

205 controller (pySDC.Controller): The controller 

206 S (pySDC.Step): The current step 

207 

208 Returns: 

209 None 

210 """ 

211 if len(S.levels) > 1 and len(controller.MS) > 1: 

212 raise NotImplementedError( 

213 "Embedded error estimate only works for serial multi-level or parallel single \ 

214level" 

215 ) 

216 

217 if S.status.iter > 0 or self.params.sweeper_type == "RK": 

218 if self.params.averaged: 

219 averaging = float(S.status.slot + 1) 

220 else: 

221 averaging = 1.0 

222 

223 for L in S.levels: 

224 temp = self.estimate_embedded_error_serial(L) 

225 L.status.error_embedded_estimate = max( 

226 [abs(temp - self.buffers.e_em_last) / averaging, np.finfo(float).eps] 

227 ) 

228 

229 if not self.params.averaged: 

230 self.buffers.e_em_last = temp * 1.0 

231 

232 return None 

233 

234 

235class EstimateEmbeddedErrorLinearizedMPI(EstimateEmbeddedError): 

236 def __init__(self, controller, params, description, **kwargs): 

237 """ 

238 Initialisation routine. Add the buffers for communication. 

239 

240 Args: 

241 controller (pySDC.Controller): The controller 

242 params (dict): Parameters for the convergence controller 

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

244 """ 

245 super().__init__(controller, params, description, **kwargs) 

246 self.buffers = Pars({'e_em_last': 0.0}) 

247 

248 def post_iteration_processing(self, controller, S, **kwargs): 

249 """ 

250 Compute embedded error estimate on the last node of each level 

251 In serial this is the local error, but in block Gauss-Seidel MSSDC this is a semi-global error in each block 

252 

253 Args: 

254 controller (pySDC.Controller): The controller 

255 S (pySDC.Step): The current step 

256 

257 Returns: 

258 None 

259 """ 

260 comm = kwargs['comm'] 

261 

262 if S.status.iter > 0 or self.params.sweeper_type == "RK": 

263 for L in S.levels: 

264 # get accumulated local errors from previous steps 

265 if not S.status.first: 

266 if not S.status.prev_done: 

267 self.buffers.e_em_last = self.recv(comm, S.status.slot - 1) 

268 else: 

269 self.buffers.e_em_last = 0.0 

270 

271 # estimate accumulated local error 

272 temp = self.estimate_embedded_error_serial(L) 

273 

274 # estimate local error as difference of accumulated errors 

275 L.status.error_embedded_estimate = max([abs(temp - self.buffers.e_em_last), np.finfo(float).eps]) 

276 

277 # send the accumulated local errors forward 

278 if not S.status.last: 

279 self.send(comm, dest=S.status.slot + 1, data=temp, blocking=True) 

280 

281 return None 

282 

283 

284class EstimateEmbeddedErrorCollocation(ConvergenceController): 

285 """ 

286 Estimates an embedded error based on changing the underlying quadrature rule. The error estimate is stored as 

287 `error_embedded_estimate_collocation` in the status of the level. Note that we only compute the estimate on the 

288 finest level. The error is stored as a tuple with the first index denoting to which iteration it belongs. This 

289 is useful since the error estimate is not available immediately after, but only when the next collocation problem 

290 is converged to make sure the two solutions are of different accuracy. 

291 

292 Changing the collocation method between iterations happens using the `AdaptiveCollocation` convergence controller. 

293 Please refer to that for documentation on how to use this. Just pass the parameters for that convergence controller 

294 as `adaptive_coll_params` to the parameters for this one and they will be passed on when the `AdaptiveCollocation` 

295 convergence controller is automatically added while loading dependencies. 

296 """ 

297 

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

299 """ 

300 Add a default value for control order to the parameters 

301 

302 Args: 

303 controller (pySDC.Controller): The controller 

304 params (dict): Parameters for the convergence controller 

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

306 

307 Returns: 

308 dict: Updated parameters 

309 """ 

310 defaults = { 

311 "control_order": 210, 

312 "adaptive_coll_params": {}, 

313 **super().setup(controller, params, description, **kwargs), 

314 } 

315 return defaults 

316 

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

318 """ 

319 Load the `AdaptiveCollocation` convergence controller to switch between collocation problems between iterations. 

320 

321 Args: 

322 controller (pySDC.Controller): The controller 

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

324 """ 

325 from pySDC.implementations.convergence_controller_classes.adaptive_collocation import AdaptiveCollocation 

326 

327 controller.add_convergence_controller( 

328 AdaptiveCollocation, params=self.params.adaptive_coll_params, description=description 

329 ) 

330 from pySDC.implementations.hooks.log_embedded_error_estimate import LogEmbeddedErrorEstimate 

331 

332 controller.add_hook(LogEmbeddedErrorEstimate) 

333 

334 def post_iteration_processing(self, controller, step, **kwargs): 

335 """ 

336 Compute the embedded error as the difference between the interpolated and the current solution on the finest 

337 level. 

338 

339 Args: 

340 controller (pySDC.Controller.controller): The controller 

341 step (pySDC.Step.step): The current step 

342 """ 

343 if step.status.done: 

344 lvl = step.levels[0] 

345 lvl.sweep.compute_end_point() 

346 self.status.u += [lvl.uend] 

347 self.status.iter += [step.status.iter] 

348 

349 if len(self.status.u) > 1: 

350 lvl.status.error_embedded_estimate_collocation = ( 

351 self.status.iter[-2], 

352 max([np.finfo(float).eps, abs(self.status.u[-1] - self.status.u[-2])]), 

353 ) 

354 

355 def setup_status_variables(self, *args, **kwargs): 

356 """ 

357 Add the embedded error variable to the levels and add a status variable for previous steps. 

358 

359 Args: 

360 controller (pySDC.Controller): The controller 

361 """ 

362 self.status = Status(['u', 'iter']) 

363 self.status.u = [] # the solutions of converged collocation problems 

364 self.status.iter = [] # the iteration in which the solution converged 

365 

366 self.add_status_variable_to_level('error_embedded_estimate_collocation')