Coverage for pySDC/implementations/convergence_controller_classes/estimate_embedded_error.py: 78%
115 statements
« prev ^ index » next coverage.py v7.17.0, created at 2026-11-02 12:30 +0000
1import numpy as np
3from pySDC.core.convergence_controller import ConvergenceController, Pars
4from pySDC.implementations.convergence_controller_classes.store_uold import StoreUOld
6from pySDC.implementations.sweeper_classes.Runge_Kutta import RungeKutta
9class EstimateEmbeddedError(ConvergenceController):
10 """
11 Estimate the local error as the difference between two solutions of different order, such as two consecutive sweeps.
12 The embedded error is obtained by computing two solutions of different accuracy and pretending the more accurate
13 one is an exact solution from the point of view of the less accurate solution. In practice, we like to compute the
14 solutions with different order methods, meaning that in SDC we can just subtract two consecutive sweeps, as long as
15 you make sure your preconditioner is compatible, which you have to just try out...
16 """
18 @classmethod
19 def get_implementation(cls, flavor='standard', useMPI=False):
20 """
21 Retrieve the implementation for a specific flavor of this class.
23 Args:
24 flavor (str): The implementation that you want
25 useMPI (bool): Whether to return the MPI version, only relevant for the `'linearized'` flavor
27 Returns:
28 cls: The child class that implements the desired flavor
29 """
30 if flavor == 'standard':
31 return cls
32 elif flavor == 'linearized':
33 if useMPI:
34 return EstimateEmbeddedErrorLinearizedMPI
35 else:
36 return EstimateEmbeddedErrorLinearizedNonMPI
37 elif flavor == 'collocation':
38 return EstimateEmbeddedErrorCollocation
39 else:
40 raise NotImplementedError(f'Flavor {flavor} of EstimateEmbeddedError is not implemented!')
42 def setup(self, controller, params, description, **kwargs):
43 """
44 Add a default value for control order to the parameters and check if we are using a Runge-Kutta sweeper
46 Args:
47 controller (pySDC.Controller): The controller
48 params (dict): Parameters for the convergence controller
49 description (dict): The description object used to instantiate the controller
51 Returns:
52 dict: Updated parameters
53 """
54 sweeper_type = 'SDC'
55 if RungeKutta in description['sweeper_class'].__mro__:
56 sweeper_type = 'RK'
57 elif 'SweeperMPI' in [me.__name__ for me in description['sweeper_class'].__mro__]:
58 sweeper_type = 'MPI'
59 return {
60 "control_order": -80,
61 "sweeper_type": sweeper_type,
62 "rel_error": False,
63 **super().setup(controller, params, description, **kwargs),
64 }
66 def dependencies(self, controller, description, **kwargs):
67 """
68 Load the convergence controller that stores the solution of the last sweep unless we are doing Runge-Kutta.
69 Add the hook for recording the error.
71 Args:
72 controller (pySDC.Controller): The controller
73 description (dict): The description object used to instantiate the controller
75 Returns:
76 None
77 """
78 if RungeKutta not in description["sweeper_class"].__bases__:
79 controller.add_convergence_controller(StoreUOld, description=description)
81 from pySDC.implementations.hooks.log_embedded_error_estimate import LogEmbeddedErrorEstimate
83 controller.add_hook(LogEmbeddedErrorEstimate)
84 return None
86 def estimate_embedded_error_serial(self, L):
87 """
88 Estimate the serial embedded error, which may need to be modified for a parallel estimate.
90 Depending on the type of sweeper, the lower order solution is stored in a different place.
92 Args:
93 L (pySDC.level): The level
95 Returns:
96 dtype_u: The embedded error estimate
97 """
98 if self.params.sweeper_type == "RK":
99 L.sweep.compute_end_point()
100 if self.params.rel_error:
101 return abs(L.uend - L.sweep.u_secondary) / abs(L.uend)
102 else:
103 return abs(L.uend - L.sweep.u_secondary)
104 elif self.params.sweeper_type == "SDC":
105 # order rises by one between sweeps
106 if self.params.rel_error:
107 return abs(L.uold[-1] - L.u[-1]) / abs(L.u[-1])
108 else:
109 return abs(L.uold[-1] - L.u[-1])
110 elif self.params.sweeper_type == 'MPI':
111 comm = L.sweep.comm
112 if self.params.rel_error:
113 return comm.bcast(
114 abs(L.uold[comm.rank + 1] - L.u[comm.rank + 1]) / abs(L.u[comm.rank + 1]), root=comm.size - 1
115 )
116 else:
117 return comm.bcast(abs(L.uold[comm.rank + 1] - L.u[comm.rank + 1]), root=comm.size - 1)
118 else:
119 raise NotImplementedError(f"Don't know how to estimate embedded error for sweeper type \
120\"{self.params.sweeper_type}\"")
122 def setup_status_variables(self, controller, **kwargs):
123 """
124 Add the embedded error to the level status
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')
132 def post_iteration_processing(self, controller, S, **kwargs):
133 """
134 Estimate the local error here.
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.
139 Args:
140 controller (pySDC.Controller): The controller
141 S (pySDC.Step): The current step
143 Returns:
144 None
145 """
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)
153 return None
156class EstimateEmbeddedErrorLinearizedNonMPI(EstimateEmbeddedError):
157 """
158 Local error in a block of steps as the embedded estimate minus that of the step before, for the non-MPI controller.
159 """
161 def __init__(self, controller, params, description, **kwargs):
162 """
163 Initialisation routine. Add the buffers for communication.
165 Args:
166 controller (pySDC.Controller): The controller
167 params (dict): Parameters for the convergence controller
168 description (dict): The description object used to instantiate the controller
169 """
170 super().__init__(controller, params, description, **kwargs)
171 self.buffers = Pars({'e_em_last': 0.0})
173 def setup(self, controller, params, description, **kwargs):
174 """
175 Add option for averaging the local errors.
177 Args:
178 controller (pySDC.Controller): The controller
179 params (dict): Parameters for the convergence controller
180 description (dict): The description object used to instantiate the controller
182 Returns:
183 dict: Updated parameters
184 """
185 return {
186 'averaged': False,
187 **super().setup(controller, params, description, **kwargs),
188 }
190 def reset_buffers_nonMPI(self, controller, **kwargs):
191 """
192 Reset buffers for imitated communication.
194 Args:
195 controller (pySDC.controller): The controller
197 Returns:
198 None
199 """
200 self.buffers.e_em_last = 0.0
201 return None
203 def post_iteration_processing(self, controller, S, **kwargs):
204 """
205 Compute embedded error estimate on the last node of each level
206 In serial this is the local error, but in block Gauss-Seidel MSSDC this is a semi-global error in each block
208 Args:
209 controller (pySDC.Controller): The controller
210 S (pySDC.Step): The current step
212 Returns:
213 None
214 """
215 if len(S.levels) > 1 and len(controller.MS) > 1:
216 raise NotImplementedError("Embedded error estimate only works for serial multi-level or parallel single \
217level")
219 if S.status.iter > 0 or self.params.sweeper_type == "RK":
220 if self.params.averaged:
221 averaging = float(S.status.slot + 1)
222 else:
223 averaging = 1.0
225 for L in S.levels:
226 temp = self.estimate_embedded_error_serial(L)
227 L.status.error_embedded_estimate = max(
228 [abs(temp - self.buffers.e_em_last) / averaging, np.finfo(float).eps]
229 )
231 if not self.params.averaged:
232 self.buffers.e_em_last = temp * 1.0
234 return None
237class EstimateEmbeddedErrorLinearizedMPI(EstimateEmbeddedError):
238 """
239 Local error in a block of steps as the embedded estimate minus that of the step before, for the MPI controller.
240 """
242 def __init__(self, controller, params, description, **kwargs):
243 """
244 Initialisation routine. Add the buffers for communication.
246 Args:
247 controller (pySDC.Controller): The controller
248 params (dict): Parameters for the convergence controller
249 description (dict): The description object used to instantiate the controller
250 """
251 super().__init__(controller, params, description, **kwargs)
252 self.buffers = Pars({'e_em_last': 0.0})
254 def post_iteration_processing(self, controller, S, **kwargs):
255 """
256 Compute embedded error estimate on the last node of each level
257 In serial this is the local error, but in block Gauss-Seidel MSSDC this is a semi-global error in each block
259 Args:
260 controller (pySDC.Controller): The controller
261 S (pySDC.Step): The current step
263 Returns:
264 None
265 """
266 comm = kwargs['comm']
268 if S.status.iter > 0 or self.params.sweeper_type == "RK":
269 for L in S.levels:
270 # get accumulated local errors from previous steps
271 if not S.status.first:
272 if not S.status.prev_done:
273 self.buffers.e_em_last = self.recv(comm, S.status.slot - 1)
274 else:
275 self.buffers.e_em_last = 0.0
277 # estimate accumulated local error
278 temp = self.estimate_embedded_error_serial(L)
280 # estimate local error as difference of accumulated errors
281 L.status.error_embedded_estimate = max([abs(temp - self.buffers.e_em_last), np.finfo(float).eps])
283 # send the accumulated local errors forward
284 if not S.status.last:
285 self.send(comm, dest=S.status.slot + 1, data=temp, blocking=True)
287 return None
290class EstimateEmbeddedErrorCollocation(ConvergenceController):
291 """
292 Estimates an embedded error based on changing the underlying quadrature rule. The error estimate is stored as
293 `error_embedded_estimate_collocation` in the status of the level. Note that we only compute the estimate on the
294 finest level. The error is stored as a tuple with the first index denoting to which iteration it belongs. This
295 is useful since the error estimate is not available immediately after, but only when the next collocation problem
296 is converged to make sure the two solutions are of different accuracy.
298 Changing the collocation method between iterations happens using the `AdaptiveCollocation` convergence controller.
299 Please refer to that for documentation on how to use this. Just pass the parameters for that convergence controller
300 as `adaptive_coll_params` to the parameters for this one and they will be passed on when the `AdaptiveCollocation`
301 convergence controller is automatically added while loading dependencies.
302 """
304 def setup(self, controller, params, description, **kwargs):
305 """
306 Add a default value for control order to the parameters
308 Args:
309 controller (pySDC.Controller): The controller
310 params (dict): Parameters for the convergence controller
311 description (dict): The description object used to instantiate the controller
313 Returns:
314 dict: Updated parameters
315 """
316 defaults = {
317 "control_order": 210,
318 "adaptive_coll_params": {},
319 **super().setup(controller, params, description, **kwargs),
320 }
321 return defaults
323 def dependencies(self, controller, description, **kwargs):
324 """
325 Load the `AdaptiveCollocation` convergence controller to switch between collocation problems between iterations.
327 Args:
328 controller (pySDC.Controller): The controller
329 description (dict): The description object used to instantiate the controller
330 """
331 from pySDC.implementations.convergence_controller_classes.adaptive_collocation import AdaptiveCollocation
333 controller.add_convergence_controller(
334 AdaptiveCollocation, params=self.params.adaptive_coll_params, description=description
335 )
336 from pySDC.implementations.hooks.log_embedded_error_estimate import LogEmbeddedErrorEstimate
338 controller.add_hook(LogEmbeddedErrorEstimate)
340 def post_iteration_processing(self, controller, step, **kwargs):
341 """
342 Compute the embedded error as the difference between the interpolated and the current solution on the finest
343 level.
345 Args:
346 controller (pySDC.Controller.controller): The controller
347 step (pySDC.Step.step): The current step
348 """
349 if step.status.done:
350 lvl = step.levels[0]
351 lvl.sweep.compute_end_point()
352 step.status.coll_u += (lvl.uend,)
353 step.status.coll_iter += (step.status.iter,)
355 if len(step.status.coll_u) > 1:
356 lvl.status.error_embedded_estimate_collocation = (
357 step.status.coll_iter[-2],
358 max([np.finfo(float).eps, abs(step.status.coll_u[-1] - step.status.coll_u[-2])]),
359 )
361 def setup_status_variables(self, *args, **kwargs):
362 """
363 Add the embedded error variable to the levels and add status variables for the collocation problems solved in
364 the current step.
366 These live on the step rather than on this object, since there is one instance of this class per controller,
367 which is a whole block of steps when running with more than one process. They are tuples, so that the same empty
368 value can be set on all steps and appending to it on one step does not change it on the others.
370 Args:
371 controller (pySDC.Controller): The controller
372 """
373 self.add_status_variable_to_step('coll_u', ()) # the solutions of converged collocation problems
374 self.add_status_variable_to_step('coll_iter', ()) # the iteration in which the solution converged
376 self.add_status_variable_to_level('error_embedded_estimate_collocation')
378 def reset_status_variables(self, *args, **kwargs):
379 """
380 Forget the collocation problems solved in the previous step, or the previous attempt of this step.
382 Args:
383 controller (pySDC.Controller): The controller
384 """
385 self.set_step_status_variable('coll_u', ())
386 self.set_step_status_variable('coll_iter', ())