Coverage for pySDC/implementations/convergence_controller_classes/estimate_extrapolation_error.py: 98%
163 statements
« prev ^ index » next coverage.py v7.17.0, created at 2026-11-02 12:30 +0000
1import numpy as np
2from scipy.special import factorial
4from pySDC.core.convergence_controller import ConvergenceController, Status
5from pySDC.core.errors import DataError
6from pySDC.implementations.datatype_classes.mesh import mesh, imex_mesh
7from pySDC.implementations.hooks.log_extrapolated_error_estimate import LogExtrapolationErrorEstimate
10class EstimateExtrapolationErrorBase(ConvergenceController):
11 """
12 Abstract base class for error estimates by Taylor extrapolation of solutions and right-hand sides from other times.
14 This error estimate extrapolates a solution based on Taylor expansions using solutions at other times. Child
15 classes decide where these come from: EstimateExtrapolationErrorNonMPI uses previous time steps, and
16 EstimateExtrapolationErrorWithinQ the collocation nodes of the current step. There is no MPI version.
17 """
19 def __init__(self, controller, params, description, **kwargs):
20 """
21 Initialization routine
23 Args:
24 controller (pySDC.Controller): The controller
25 params (dict): The params passed for this specific convergence controller
26 description (dict): The description object used to instantiate the controller
27 """
28 self.prev = Status(["t", "u", "f", "dt"]) # store solutions etc. of previous steps here
29 self.coeff = Status(["u", "f", "prefactor"]) # store coefficients for extrapolation here
30 super().__init__(controller, params, description)
31 controller.add_hook(LogExtrapolationErrorEstimate)
33 def setup(self, controller, params, description, **kwargs):
34 """
35 The extrapolation based method requires storage of previous values of u, f, t and dt and also requires solving
36 a linear system of equations to compute the Taylor expansion finite difference style. Here, all variables are
37 initialized which are needed for this process.
39 Args:
40 controller (pySDC.Controller): The controller
41 params (dict): The params passed for this specific convergence controller
42 description (dict): The description object used to instantiate the controller
44 Returns:
45 dict: Updated parameters with default values
46 """
48 from pySDC.implementations.convergence_controller_classes.hotrod import HotRod
49 from pySDC.implementations.convergence_controller_classes.adaptivity import (
50 Adaptivity,
51 )
53 default_params = {
54 "control_order": -75,
55 "use_adaptivity": True in [me == Adaptivity for me in description.get("convergence_controllers", {})],
56 "use_HotRod": True in [me == HotRod for me in description.get("convergence_controllers", {})],
57 "order_time_marching": description["step_params"]["maxiter"],
58 }
60 new_params = {**default_params, **super().setup(controller, params, description, **kwargs)}
62 # Do a sufficiently high order Taylor expansion
63 new_params["Taylor_order"] = new_params["order_time_marching"] + 2
65 # Estimate and store values from this iteration
66 new_params["estimate_iter"] = new_params["order_time_marching"] - (1 if new_params["use_HotRod"] else 0)
68 # Store n values. Since we store u and f, we need only half of each (the +1 is for rounding)
69 new_params["n"] = (new_params["Taylor_order"] + 1) // 2
70 new_params["n_per_proc"] = new_params["n"] * 1
72 return new_params
74 def setup_status_variables(self, controller, **kwargs):
75 """
76 Initialize coefficient variables and add variable to the levels for extrapolated error
78 Args:
79 controller (pySDC.controller): The controller
81 Returns:
82 None
83 """
84 self.coeff.u = [None] * self.params.n
85 self.coeff.f = [0.0] * self.params.n
87 self.add_status_variable_to_level('error_extrapolation_estimate')
89 def check_parameters(self, controller, params, description, **kwargs):
90 """
91 Check whether parameters are compatible with whatever assumptions went into the step size functions etc.
93 Args:
94 controller (pySDC.Controller): The controller
95 params (dict): The params passed for this specific convergence controller
96 description (dict): The description object used to instantiate the controller
98 Returns:
99 bool: Whether the parameters are compatible
100 str: Error message
101 """
102 if controller.params.mssdc_jac:
103 return (
104 False,
105 "Extrapolation error estimator needs the same order on all steps, please activate Gauss-Seid\
106el multistep mode!",
107 )
109 return True, ""
111 def store_values(self, S, **kwargs):
112 """
113 Store the required attributes of the step to do the extrapolation. We only care about the last collocation
114 node on the finest level at the moment.
116 Args:
117 S (pySDC.Step): The current step
119 Returns:
120 None
121 """
122 # figure out which values are to be replaced by the new ones
123 if None in self.prev.t:
124 oldest_val = len(self.prev.t) - len(self.prev.t[self.prev.t == [None]])
125 else:
126 oldest_val = np.argmin(self.prev.t)
128 # figure out how to store the right hand side
129 f = S.levels[0].f[-1]
130 if type(f) == imex_mesh:
131 self.prev.f[oldest_val] = f.impl + f.expl
132 elif type(f) == mesh:
133 self.prev.f[oldest_val] = f
134 else:
135 raise DataError(f"Unable to store f from datatype {type(f)}, extrapolation based error estimate only\
136 works with types imex_mesh and mesh")
138 # store the rest of the values
139 self.prev.u[oldest_val] = S.levels[0].u[-1]
140 self.prev.t[oldest_val] = S.time + S.dt
141 self.prev.dt[oldest_val] = S.dt
143 return None
145 def get_extrapolation_coefficients(self, t, dt, t_eval):
146 """
147 This function solves a linear system where in the matrix A, the row index reflects the order of the derivative
148 in the Taylor expansion and the column index reflects the particular step and whether its u or f from that
149 step. The vector b on the other hand, contains a 1 in the first entry and zeros elsewhere, since we want to
150 compute the value itself and all the derivatives should vanish after combining the Taylor expansions. This
151 works to the order of the number of rows and since we want a square matrix for solving, we need the same amount
152 of columns, which determines the memory overhead, since it is equal to the solutions / rhs that we need in
153 memory at the time of evaluation.
155 This is enough to get the extrapolated solution, but if we want to compute the local error, we have to compute
156 a prefactor. This is based on error accumulation between steps (first step's solution is exact plus 1 LTE,
157 second solution is exact plus 2 LTE and so on), which can be computed for adaptive step sizes as well. However,
158 this is only true for linear problems, which means we expect the error estimate to work less well for non-linear
159 problems.
161 Since only time differences are important for computing the coefficients, we need to compute this only once when
162 using constant step sizes. When we allow the step size to change, however, we need to recompute this in every
163 step, which is activated by the `use_adaptivity` parameter.
165 Solving for the coefficients requires solving a dense linear system of equations. The number of unknowns is
166 equal to the order of the Taylor expansion, so this step should be cheap compared to the solves in each SDC
167 iteration.
169 The function stores the computed coefficients in the `self.coeff` variables.
171 Args:
172 t (list): The list of times at which we have solutions available
173 dt (list): The step sizes used for computing these solutions (needed for the prefactor)
174 t_eval (float): The time we want to extrapolate to
176 Returns:
177 None
178 """
180 # prepare A matrix
181 A = np.zeros((self.params.Taylor_order, self.params.Taylor_order))
182 A[0, 0 : self.params.n] = 1.0
183 j = np.arange(self.params.Taylor_order)
184 inv_facs = 1.0 / factorial(j)
186 # get the steps backwards from the point of evaluation
187 idx = np.argsort(t)
188 steps_from_now = t[idx] - t_eval
190 # fill A matrix
191 for i in range(1, self.params.Taylor_order):
192 # Taylor expansions of the solutions
193 A[i, : self.params.n] = steps_from_now ** j[i] * inv_facs[i]
195 # Taylor expansions of the first derivatives a.k.a. right hand side evaluations
196 A[i, self.params.n : self.params.Taylor_order] = (
197 steps_from_now[2 * self.params.n - self.params.Taylor_order :] ** (j[i] - 1) * inv_facs[i - 1]
198 )
200 # prepare rhs
201 b = np.zeros(self.params.Taylor_order)
202 b[0] = 1.0
204 # solve linear system for the coefficients
205 coeff = np.linalg.solve(A, b)
206 self.coeff.u = coeff[: self.params.n]
207 self.coeff.f[self.params.n * 2 - self.params.Taylor_order :] = coeff[self.params.n : self.params.Taylor_order]
209 # determine prefactor
210 step_size_ratios = abs(dt[len(dt) - len(self.coeff.u) :] / dt[-1]) ** (self.params.Taylor_order - 1)
211 inv_prefactor = -sum(step_size_ratios[1:]) - 1.0
212 for i in range(len(self.coeff.u)):
213 inv_prefactor += sum(step_size_ratios[1 : i + 1]) * self.coeff.u[i]
214 self.coeff.prefactor = 1.0 / abs(inv_prefactor)
216 return None
219class EstimateExtrapolationErrorNonMPI(EstimateExtrapolationErrorBase):
220 """
221 Implementation of the extrapolation error estimate for the non-MPI controller.
222 """
224 def check_parameters(self, controller, params, description, **kwargs):
225 """
226 Check whether parameters are compatible with whatever assumptions went into the step size functions etc.
228 Args:
229 controller (pySDC.Controller): The controller
230 params (dict): The params passed for this specific convergence controller
231 description (dict): The description object used to instantiate the controller
233 Returns:
234 bool: Whether the parameters are compatible
235 str: Error message
236 """
237 if description["level_params"].get("restol", -1.0) >= 0:
238 return (
239 False,
240 "Extrapolation error needs constant order in time and hence restol in the level parameters \
241has to be smaller than 0!",
242 )
244 return super().check_parameters(controller, params, description, **kwargs)
246 def setup(self, controller, params, description, **kwargs):
247 """
248 Add a no parameter 'no_storage' which decides whether the standard or the no-memory-overhead version is run,
249 where only values are used for extrapolation which are in memory of other processes
251 Args:
252 controller (pySDC.Controller): The controller
253 params (dict): The params passed for this specific convergence controller
254 description (dict): The description object used to instantiate the controller
256 Returns:
257 dict: Updated parameters with default values
258 """
259 default_params = super().setup(controller, params, description)
261 non_mpi_defaults = {
262 "no_storage": False,
263 }
265 return {**non_mpi_defaults, **default_params}
267 def setup_status_variables(self, controller, **kwargs):
268 """
269 Initialize storage variables.
271 Args:
272 controller (pySDC.controller): The controller
274 Returns:
275 None
276 """
277 super().setup_status_variables(controller, **kwargs)
279 self.prev.t = np.array([None] * self.params.n)
280 self.prev.dt = np.array([None] * self.params.n)
281 self.prev.u = [None] * self.params.n
282 self.prev.f = [None] * self.params.n
284 return None
286 def post_iteration_processing(self, controller, S, **kwargs):
287 """
288 We perform three key operations here in the last iteration:
289 - Compute the error estimate
290 - Compute the coefficients if needed
291 - Store the values of the step if we pretend not to for the no-memory overhead version
293 Args:
294 controller (pySDC.Controller): The controller
295 S (pySDC.Step): The current step
297 Returns:
298 None
299 """
300 if S.status.iter == self.params.estimate_iter:
301 t_eval = S.time + S.dt
303 # compute the extrapolation coefficients if needed
304 if (
305 (
306 None in self.coeff.u
307 or self.params.use_adaptivity
308 or (not self.params.no_storage and S.status.time_size > 1)
309 )
310 and None not in self.prev.t
311 and t_eval > max(self.prev.t)
312 ):
313 self.get_extrapolation_coefficients(self.prev.t, self.prev.dt, t_eval)
315 # compute the error if we can
316 if None not in self.coeff.u and None not in self.prev.t:
317 self.get_extrapolated_error(S)
319 # store the solution and pretend we didn't because in the non MPI version we take a few shortcuts
320 if self.params.no_storage:
321 self.store_values(S)
323 return None
325 def prepare_next_block(self, controller, S, size, time, Tend, MS, **kwargs):
326 """
327 If the no-memory-overhead version is used, we need to delete stuff that shouldn't be available. Otherwise, we
328 need to store all the stuff that we can.
330 Args:
331 controller (pySDC.Controller): The controller
332 S (pySDC.step): The current step
333 size (int): Number of ranks
334 time (float): The current time
335 Tend (float): The final time
336 MS (list): Active steps
338 Returns:
339 None
340 """
342 # delete values that should not be available in the next step
343 if self.params.no_storage:
344 self.prev.t = np.array([None] * self.params.n)
345 self.prev.dt = np.array([None] * self.params.n)
346 self.prev.u = [None] * self.params.n
347 self.prev.f = [None] * self.params.n
349 else:
350 # decide where we need to restart to store everything up to that point
351 restarts = [S.status.restart for S in MS]
352 restart_at = np.where(restarts)[0][0] if True in restarts else len(MS)
354 # store values in the current block that don't need restarting
355 if restart_at > S.status.slot:
356 self.store_values(S)
358 return None
360 def get_extrapolated_solution(self, S, **kwargs):
361 """
362 Combine values from previous steps to extrapolate.
364 Args:
365 S (pySDC.Step): The current step
367 Returns:
368 dtype_u: The extrapolated solution
369 """
370 if len(S.levels) > 1:
371 raise NotImplementedError("Extrapolated estimate only works on the finest level for now")
373 # prepare variables
374 u_ex = S.levels[0].u[-1] * 0.0
375 idx = np.argsort(self.prev.t)
377 # see if we have a solution for the current step already stored
378 if (abs(S.time + S.dt - self.prev.t) < 10.0 * np.finfo(float).eps).any():
379 idx_step = idx[np.argmin(abs(self.prev.t - S.time - S.dt))]
380 else:
381 idx_step = max(idx) + 1
383 # make a mask of all the steps we want to include in the extrapolation
384 mask = np.logical_and(idx < idx_step, idx >= idx_step - self.params.n)
386 # do the extrapolation by summing everything up
387 for i in range(self.params.n):
388 u_ex += self.coeff.u[i] * self.prev.u[idx[mask][i]] + self.coeff.f[i] * self.prev.f[idx[mask][i]]
390 return u_ex
392 def get_extrapolated_error(self, S, **kwargs):
393 """
394 The extrapolation estimate combines values of u and f from multiple steps to extrapolate and compare to the
395 solution obtained by the time marching scheme.
397 Args:
398 S (pySDC.Step): The current step
400 Returns:
401 None
402 """
403 u_ex = self.get_extrapolated_solution(S)
404 if u_ex is not None:
405 difference = u_ex - S.levels[0].u[-1]
407 # algebraic components of a DAE have no time derivative to extrapolate
408 diff_mask = getattr(S.levels[0].prob, 'diff_mask', None)
409 if diff_mask is not None and not all(diff_mask):
410 difference = difference[diff_mask]
412 S.levels[0].status.error_extrapolation_estimate = abs(difference) * self.coeff.prefactor
413 else:
414 S.levels[0].status.error_extrapolation_estimate = None
417class EstimateExtrapolationErrorWithinQ(EstimateExtrapolationErrorBase):
418 """
419 Estimate the local error by comparing the converged SDC solution to an extrapolation from other quadrature nodes.
420 Collocation methods compute a high order solution from a linear combination
421 of solutions at intermediate time points. While the intermediate solutions (a.k.a. stages) don't share the order of
422 accuracy with the solution at the end of the interval, for SDC we know that the order is equal to the number of
423 nodes + 1 (locally). This is because the solution to the collocation problem is a polynomial approximation of order
424 of the number of nodes.
425 That means we can do a Taylor expansion around the end point of the interval to higher order and after cancelling
426 terms just like we are used to with the extrapolation based error estimate across multiple steps, we get an error
427 estimate that is of the order of accuracy of the stages.
428 This can be used for adaptivity, for instance, with the nice property that it doesn't matter how we arrived at the
429 converged collocation solution, as long as we did. We don't rely on knowing the order of accuracy after every sweep,
430 only after convergence of the collocation problem has been achieved, which we can check from the residual.
431 """
433 def setup(self, controller, params, description, **kwargs):
434 """
435 We need this convergence controller to become active after the check for convergence, because we need the step
436 to be converged.
438 Args:
439 controller (pySDC.Controller): The controller
440 params (dict): The params passed for this specific convergence controller
441 description (dict): The description object used to instantiate the controller
443 Returns:
444 dict: Updated parameters with default values
445 """
446 from pySDC.core.check_convergence import CheckConvergence
448 num_nodes = description['sweeper_params']['num_nodes']
449 self.comm = description['sweeper_params'].get('comm', None)
450 if self.comm:
451 from mpi4py import MPI
453 self.sum = MPI.SUM
455 self.check_convergence = CheckConvergence.check_convergence
457 default_params = {
458 'Taylor_order': 2 * num_nodes,
459 'n': num_nodes,
460 'recompute_coefficients': False,
461 **params,
462 }
464 return {**super().setup(controller, params, description, **kwargs), **default_params}
466 def post_iteration_processing(self, controller, S, **kwargs):
467 """
468 Compute the extrapolated error estimate here if the step is converged.
470 Args:
471 controller (pySDC.Controller): The controller
472 S (pySDC.Step): The current step
474 Returns:
475 None
476 """
478 if not self.check_convergence(S):
479 return None
481 lvl = S.levels[0]
483 nodes_ = lvl.sweep.coll.nodes * S.dt
484 nodes = S.time + np.append(0, nodes_[:-1])
485 t_eval = S.time + nodes_[-1]
487 dts = np.append(nodes_[0], nodes_[1:] - nodes_[:-1])
488 self.params.Taylor_order = len(nodes)
489 self.params.n = len(nodes)
491 # compute the extrapolation coefficients
492 if None in self.coeff.u or self.params.recompute_coefficients:
493 self.get_extrapolation_coefficients(nodes, dts, t_eval)
495 # compute the extrapolated solution
496 if lvl.f[0] is None:
497 lvl.f[0] = lvl.prob.eval_f(lvl.u[0], lvl.time)
499 if type(lvl.f[0]) == imex_mesh:
500 f = [lvl.f[i].impl + lvl.f[i].expl if self.coeff.f[i] and lvl.f[i] else 0.0 for i in range(len(lvl.f) - 1)]
501 elif type(lvl.f[0]) == mesh:
502 f = [lvl.f[i] if self.coeff.f[i] else 0.0 for i in range(len(lvl.f) - 1)]
503 else:
504 raise DataError(f"Unable to store f from datatype {type(lvl.f[0])}, extrapolation based error estimate only\
505 works with types imex_mesh and mesh")
507 # compute the error with the weighted sum
508 if self.comm:
509 idx = (self.comm.rank + 1) % self.comm.size
510 sendbuf = self.coeff.u[idx] * lvl.u[idx] + self.coeff.f[idx] * f[idx]
511 u_ex = lvl.prob.dtype_u(lvl.prob.init, val=0.0) if self.comm.rank == self.comm.size - 1 else None
512 self.comm.Reduce(sendbuf, u_ex, op=self.sum, root=self.comm.size - 1)
513 else:
514 u_ex = lvl.prob.dtype_u(lvl.prob.init, val=0.0)
515 for i in range(self.params.n):
516 u_ex += self.coeff.u[i] * lvl.u[i] + self.coeff.f[i] * f[i]
518 # store the error
519 if self.comm:
520 error = (
521 abs(u_ex - lvl.u[self.comm.rank + 1]) * self.coeff.prefactor
522 if self.comm.rank == self.comm.size - 1
523 else None
524 )
525 lvl.status.error_extrapolation_estimate = self.comm.bcast(error, root=self.comm.size - 1)
526 else:
527 lvl.status.error_extrapolation_estimate = abs(u_ex - lvl.u[-1]) * self.coeff.prefactor
528 return None