Coverage for pySDC / implementations / convergence_controller_classes / estimate_embedded_error.py: 77%
113 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-03-27 07:06 +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(f"Don't know how to estimate embedded error for sweeper type \
118\"{self.params.sweeper_type}\"")
120 def setup_status_variables(self, controller, **kwargs):
121 """
122 Add the embedded error to the level status
124 Args:
125 controller (pySDC.Controller): The controller
126 """
127 self.add_status_variable_to_level('error_embedded_estimate')
128 self.add_status_variable_to_level('increment')
130 def post_iteration_processing(self, controller, S, **kwargs):
131 """
132 Estimate the local error here.
134 If you are doing MSSDC, this is the global error within the block in Gauss-Seidel mode.
135 In Jacobi mode, I haven't thought about what this is.
137 Args:
138 controller (pySDC.Controller): The controller
139 S (pySDC.Step): The current step
141 Returns:
142 None
143 """
145 if S.status.iter > 0 or self.params.sweeper_type == "RK":
146 for L in S.levels:
147 L.status.error_embedded_estimate = max([self.estimate_embedded_error_serial(L), np.finfo(float).eps])
148 L.status.increment = L.status.error_embedded_estimate * 1
149 self.debug(f'L.status.error_embedded_estimate={L.status.error_embedded_estimate:.5e}', S)
151 return None
154class EstimateEmbeddedErrorLinearizedNonMPI(EstimateEmbeddedError):
155 def __init__(self, controller, params, description, **kwargs):
156 """
157 Initialisation routine. Add the buffers for communication.
159 Args:
160 controller (pySDC.Controller): The controller
161 params (dict): Parameters for the convergence controller
162 description (dict): The description object used to instantiate the controller
163 """
164 super().__init__(controller, params, description, **kwargs)
165 self.buffers = Pars({'e_em_last': 0.0})
167 def setup(self, controller, params, description, **kwargs):
168 """
169 Add option for averaging the local errors.
171 Args:
172 controller (pySDC.Controller): The controller
173 params (dict): Parameters for the convergence controller
174 description (dict): The description object used to instantiate the controller
176 Returns:
177 dict: Updated parameters
178 """
179 return {
180 'averaged': False,
181 **super().setup(controller, params, description, **kwargs),
182 }
184 def reset_buffers_nonMPI(self, controller, **kwargs):
185 """
186 Reset buffers for imitated communication.
188 Args:
189 controller (pySDC.controller): The controller
191 Returns:
192 None
193 """
194 self.buffers.e_em_last = 0.0
195 return None
197 def post_iteration_processing(self, controller, S, **kwargs):
198 """
199 Compute embedded error estimate on the last node of each level
200 In serial this is the local error, but in block Gauss-Seidel MSSDC this is a semi-global error in each block
202 Args:
203 controller (pySDC.Controller): The controller
204 S (pySDC.Step): The current step
206 Returns:
207 None
208 """
209 if len(S.levels) > 1 and len(controller.MS) > 1:
210 raise NotImplementedError("Embedded error estimate only works for serial multi-level or parallel single \
211level")
213 if S.status.iter > 0 or self.params.sweeper_type == "RK":
214 if self.params.averaged:
215 averaging = float(S.status.slot + 1)
216 else:
217 averaging = 1.0
219 for L in S.levels:
220 temp = self.estimate_embedded_error_serial(L)
221 L.status.error_embedded_estimate = max(
222 [abs(temp - self.buffers.e_em_last) / averaging, np.finfo(float).eps]
223 )
225 if not self.params.averaged:
226 self.buffers.e_em_last = temp * 1.0
228 return None
231class EstimateEmbeddedErrorLinearizedMPI(EstimateEmbeddedError):
232 def __init__(self, controller, params, description, **kwargs):
233 """
234 Initialisation routine. Add the buffers for communication.
236 Args:
237 controller (pySDC.Controller): The controller
238 params (dict): Parameters for the convergence controller
239 description (dict): The description object used to instantiate the controller
240 """
241 super().__init__(controller, params, description, **kwargs)
242 self.buffers = Pars({'e_em_last': 0.0})
244 def post_iteration_processing(self, controller, S, **kwargs):
245 """
246 Compute embedded error estimate on the last node of each level
247 In serial this is the local error, but in block Gauss-Seidel MSSDC this is a semi-global error in each block
249 Args:
250 controller (pySDC.Controller): The controller
251 S (pySDC.Step): The current step
253 Returns:
254 None
255 """
256 comm = kwargs['comm']
258 if S.status.iter > 0 or self.params.sweeper_type == "RK":
259 for L in S.levels:
260 # get accumulated local errors from previous steps
261 if not S.status.first:
262 if not S.status.prev_done:
263 self.buffers.e_em_last = self.recv(comm, S.status.slot - 1)
264 else:
265 self.buffers.e_em_last = 0.0
267 # estimate accumulated local error
268 temp = self.estimate_embedded_error_serial(L)
270 # estimate local error as difference of accumulated errors
271 L.status.error_embedded_estimate = max([abs(temp - self.buffers.e_em_last), np.finfo(float).eps])
273 # send the accumulated local errors forward
274 if not S.status.last:
275 self.send(comm, dest=S.status.slot + 1, data=temp, blocking=True)
277 return None
280class EstimateEmbeddedErrorCollocation(ConvergenceController):
281 """
282 Estimates an embedded error based on changing the underlying quadrature rule. The error estimate is stored as
283 `error_embedded_estimate_collocation` in the status of the level. Note that we only compute the estimate on the
284 finest level. The error is stored as a tuple with the first index denoting to which iteration it belongs. This
285 is useful since the error estimate is not available immediately after, but only when the next collocation problem
286 is converged to make sure the two solutions are of different accuracy.
288 Changing the collocation method between iterations happens using the `AdaptiveCollocation` convergence controller.
289 Please refer to that for documentation on how to use this. Just pass the parameters for that convergence controller
290 as `adaptive_coll_params` to the parameters for this one and they will be passed on when the `AdaptiveCollocation`
291 convergence controller is automatically added while loading dependencies.
292 """
294 def setup(self, controller, params, description, **kwargs):
295 """
296 Add a default value for control order to the parameters
298 Args:
299 controller (pySDC.Controller): The controller
300 params (dict): Parameters for the convergence controller
301 description (dict): The description object used to instantiate the controller
303 Returns:
304 dict: Updated parameters
305 """
306 defaults = {
307 "control_order": 210,
308 "adaptive_coll_params": {},
309 **super().setup(controller, params, description, **kwargs),
310 }
311 return defaults
313 def dependencies(self, controller, description, **kwargs):
314 """
315 Load the `AdaptiveCollocation` convergence controller to switch between collocation problems between iterations.
317 Args:
318 controller (pySDC.Controller): The controller
319 description (dict): The description object used to instantiate the controller
320 """
321 from pySDC.implementations.convergence_controller_classes.adaptive_collocation import AdaptiveCollocation
323 controller.add_convergence_controller(
324 AdaptiveCollocation, params=self.params.adaptive_coll_params, description=description
325 )
326 from pySDC.implementations.hooks.log_embedded_error_estimate import LogEmbeddedErrorEstimate
328 controller.add_hook(LogEmbeddedErrorEstimate)
330 def post_iteration_processing(self, controller, step, **kwargs):
331 """
332 Compute the embedded error as the difference between the interpolated and the current solution on the finest
333 level.
335 Args:
336 controller (pySDC.Controller.controller): The controller
337 step (pySDC.Step.step): The current step
338 """
339 if step.status.done:
340 lvl = step.levels[0]
341 lvl.sweep.compute_end_point()
342 self.status.u += [lvl.uend]
343 self.status.iter += [step.status.iter]
345 if len(self.status.u) > 1:
346 lvl.status.error_embedded_estimate_collocation = (
347 self.status.iter[-2],
348 max([np.finfo(float).eps, abs(self.status.u[-1] - self.status.u[-2])]),
349 )
351 def setup_status_variables(self, *args, **kwargs):
352 """
353 Add the embedded error variable to the levels and add a status variable for previous steps.
355 Args:
356 controller (pySDC.Controller): The controller
357 """
358 self.status = Status(['u', 'iter'])
359 self.status.u = [] # the solutions of converged collocation problems
360 self.status.iter = [] # the iteration in which the solution converged
362 self.add_status_variable_to_level('error_embedded_estimate_collocation')