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-18 08:18 +0000
1import numpy as np
3from pySDC.core.convergence_controller import ConvergenceController, Pars, Status
4from pySDC.implementations.convergence_controller_classes.store_uold import StoreUOld
6from pySDC.implementations.sweeper_classes.Runge_Kutta import RungeKutta
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 """
17 @classmethod
18 def get_implementation(cls, flavor='standard', useMPI=False):
19 """
20 Retrieve the implementation for a specific flavor of this class.
22 Args:
23 flavor (str): The implementation that you want
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!')
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
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
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 }
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.
69 Args:
70 controller (pySDC.Controller): The controller
71 description (dict): The description object used to instantiate the controller
73 Returns:
74 None
75 """
76 if RungeKutta not in description["sweeper_class"].__bases__:
77 controller.add_convergence_controller(StoreUOld, description=description)
79 from pySDC.implementations.hooks.log_embedded_error_estimate import LogEmbeddedErrorEstimate
81 controller.add_hook(LogEmbeddedErrorEstimate)
82 return None
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.
88 Depending on the type of sweeper, the lower order solution is stored in a different place.
90 Args:
91 L (pySDC.level): The level
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 )
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 def __init__(self, controller, params, description, **kwargs):
158 """
159 Initialisation routine. Add the buffers for communication.
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})
169 def setup(self, controller, params, description, **kwargs):
170 """
171 Add option for averaging the local errors.
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
178 Returns:
179 dict: Updated parameters
180 """
181 return {
182 'averaged': False,
183 **super().setup(controller, params, description, **kwargs),
184 }
186 def reset_buffers_nonMPI(self, controller, **kwargs):
187 """
188 Reset buffers for imitated communication.
190 Args:
191 controller (pySDC.controller): The controller
193 Returns:
194 None
195 """
196 self.buffers.e_em_last = 0.0
197 return None
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
204 Args:
205 controller (pySDC.Controller): The controller
206 S (pySDC.Step): The current step
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 )
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
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 )
229 if not self.params.averaged:
230 self.buffers.e_em_last = temp * 1.0
232 return None
235class EstimateEmbeddedErrorLinearizedMPI(EstimateEmbeddedError):
236 def __init__(self, controller, params, description, **kwargs):
237 """
238 Initialisation routine. Add the buffers for communication.
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})
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
253 Args:
254 controller (pySDC.Controller): The controller
255 S (pySDC.Step): The current step
257 Returns:
258 None
259 """
260 comm = kwargs['comm']
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
271 # estimate accumulated local error
272 temp = self.estimate_embedded_error_serial(L)
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])
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)
281 return None
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.
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 """
298 def setup(self, controller, params, description, **kwargs):
299 """
300 Add a default value for control order to the parameters
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
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
317 def dependencies(self, controller, description, **kwargs):
318 """
319 Load the `AdaptiveCollocation` convergence controller to switch between collocation problems between iterations.
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
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
332 controller.add_hook(LogEmbeddedErrorEstimate)
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.
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]
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 )
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.
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
366 self.add_status_variable_to_level('error_embedded_estimate_collocation')