Coverage for pySDC/implementations/convergence_controller_classes/estimate_extrapolation_error.py: 97%
157 statements
« prev ^ index » next coverage.py v7.6.9, created at 2024-12-20 14:51 +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 extrapolated error estimates
13 ----------------------------------------------------
14 This error estimate extrapolates a solution based on Taylor expansions using solutions of previous time steps.
15 In particular, child classes need to implement how to make these solutions available, which works differently for
16 MPI and non-MPI versions.
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 description["step_params"].get("restol", -1.0) >= 0:
103 return (
104 False,
105 "Extrapolation error needs constant order in time and hence restol in the step parameters \
106has to be smaller than 0!",
107 )
109 if controller.params.mssdc_jac:
110 return (
111 False,
112 "Extrapolation error estimator needs the same order on all steps, please activate Gauss-Seid\
113el multistep mode!",
114 )
116 return True, ""
118 def store_values(self, S, **kwargs):
119 """
120 Store the required attributes of the step to do the extrapolation. We only care about the last collocation
121 node on the finest level at the moment.
123 Args:
124 S (pySDC.Step): The current step
126 Returns:
127 None
128 """
129 # figure out which values are to be replaced by the new ones
130 if None in self.prev.t:
131 oldest_val = len(self.prev.t) - len(self.prev.t[self.prev.t == [None]])
132 else:
133 oldest_val = np.argmin(self.prev.t)
135 # figure out how to store the right hand side
136 f = S.levels[0].f[-1]
137 if type(f) == imex_mesh:
138 self.prev.f[oldest_val] = f.impl + f.expl
139 elif type(f) == mesh:
140 self.prev.f[oldest_val] = f
141 else:
142 raise DataError(
143 f"Unable to store f from datatype {type(f)}, extrapolation based error estimate only\
144 works with types imex_mesh and mesh"
145 )
147 # store the rest of the values
148 self.prev.u[oldest_val] = S.levels[0].u[-1]
149 self.prev.t[oldest_val] = S.time + S.dt
150 self.prev.dt[oldest_val] = S.dt
152 return None
154 def get_extrapolation_coefficients(self, t, dt, t_eval):
155 """
156 This function solves a linear system where in the matrix A, the row index reflects the order of the derivative
157 in the Taylor expansion and the column index reflects the particular step and whether its u or f from that
158 step. The vector b on the other hand, contains a 1 in the first entry and zeros elsewhere, since we want to
159 compute the value itself and all the derivatives should vanish after combining the Taylor expansions. This
160 works to the order of the number of rows and since we want a square matrix for solving, we need the same amount
161 of columns, which determines the memory overhead, since it is equal to the solutions / rhs that we need in
162 memory at the time of evaluation.
164 This is enough to get the extrapolated solution, but if we want to compute the local error, we have to compute
165 a prefactor. This is based on error accumulation between steps (first step's solution is exact plus 1 LTE,
166 second solution is exact plus 2 LTE and so on), which can be computed for adaptive step sizes as well. However,
167 this is only true for linear problems, which means we expect the error estimate to work less well for non-linear
168 problems.
170 Since only time differences are important for computing the coefficients, we need to compute this only once when
171 using constant step sizes. When we allow the step size to change, however, we need to recompute this in every
172 step, which is activated by the `use_adaptivity` parameter.
174 Solving for the coefficients requires solving a dense linear system of equations. The number of unknowns is
175 equal to the order of the Taylor expansion, so this step should be cheap compared to the solves in each SDC
176 iteration.
178 The function stores the computed coefficients in the `self.coeff` variables.
180 Args:
181 t (list): The list of times at which we have solutions available
182 dt (list): The step sizes used for computing these solutions (needed for the prefactor)
183 t_eval (float): The time we want to extrapolate to
185 Returns:
186 None
187 """
189 # prepare A matrix
190 A = np.zeros((self.params.Taylor_order, self.params.Taylor_order))
191 A[0, 0 : self.params.n] = 1.0
192 j = np.arange(self.params.Taylor_order)
193 inv_facs = 1.0 / factorial(j)
195 # get the steps backwards from the point of evaluation
196 idx = np.argsort(t)
197 steps_from_now = t[idx] - t_eval
199 # fill A matrix
200 for i in range(1, self.params.Taylor_order):
201 # Taylor expansions of the solutions
202 A[i, : self.params.n] = steps_from_now ** j[i] * inv_facs[i]
204 # Taylor expansions of the first derivatives a.k.a. right hand side evaluations
205 A[i, self.params.n : self.params.Taylor_order] = (
206 steps_from_now[2 * self.params.n - self.params.Taylor_order :] ** (j[i] - 1) * inv_facs[i - 1]
207 )
209 # prepare rhs
210 b = np.zeros(self.params.Taylor_order)
211 b[0] = 1.0
213 # solve linear system for the coefficients
214 coeff = np.linalg.solve(A, b)
215 self.coeff.u = coeff[: self.params.n]
216 self.coeff.f[self.params.n * 2 - self.params.Taylor_order :] = coeff[self.params.n : self.params.Taylor_order]
218 # determine prefactor
219 step_size_ratios = abs(dt[len(dt) - len(self.coeff.u) :] / dt[-1]) ** (self.params.Taylor_order - 1)
220 inv_prefactor = -sum(step_size_ratios[1:]) - 1.0
221 for i in range(len(self.coeff.u)):
222 inv_prefactor += sum(step_size_ratios[1 : i + 1]) * self.coeff.u[i]
223 self.coeff.prefactor = 1.0 / abs(inv_prefactor)
225 return None
228class EstimateExtrapolationErrorNonMPI(EstimateExtrapolationErrorBase):
229 """
230 Implementation of the extrapolation error estimate for the non-MPI controller.
231 """
233 def setup(self, controller, params, description, **kwargs):
234 """
235 Add a no parameter 'no_storage' which decides whether the standard or the no-memory-overhead version is run,
236 where only values are used for extrapolation which are in memory of other processes
238 Args:
239 controller (pySDC.Controller): The controller
240 params (dict): The params passed for this specific convergence controller
241 description (dict): The description object used to instantiate the controller
243 Returns:
244 dict: Updated parameters with default values
245 """
246 default_params = super().setup(controller, params, description)
248 non_mpi_defaults = {
249 "no_storage": False,
250 }
252 return {**non_mpi_defaults, **default_params}
254 def setup_status_variables(self, controller, **kwargs):
255 """
256 Initialize storage variables.
258 Args:
259 controller (pySDC.controller): The controller
261 Returns:
262 None
263 """
264 super().setup_status_variables(controller, **kwargs)
266 self.prev.t = np.array([None] * self.params.n)
267 self.prev.dt = np.array([None] * self.params.n)
268 self.prev.u = [None] * self.params.n
269 self.prev.f = [None] * self.params.n
271 return None
273 def post_iteration_processing(self, controller, S, **kwargs):
274 """
275 We perform three key operations here in the last iteration:
276 - Compute the error estimate
277 - Compute the coefficients if needed
278 - Store the values of the step if we pretend not to for the no-memory overhead version
280 Args:
281 controller (pySDC.Controller): The controller
282 S (pySDC.Step): The current step
284 Returns:
285 None
286 """
287 if S.status.iter == self.params.estimate_iter:
288 t_eval = S.time + S.dt
290 # compute the extrapolation coefficients if needed
291 if (
292 (
293 None in self.coeff.u
294 or self.params.use_adaptivity
295 or (not self.params.no_storage and S.status.time_size > 1)
296 )
297 and None not in self.prev.t
298 and t_eval > max(self.prev.t)
299 ):
300 self.get_extrapolation_coefficients(self.prev.t, self.prev.dt, t_eval)
302 # compute the error if we can
303 if None not in self.coeff.u and None not in self.prev.t:
304 self.get_extrapolated_error(S)
306 # store the solution and pretend we didn't because in the non MPI version we take a few shortcuts
307 if self.params.no_storage:
308 self.store_values(S)
310 return None
312 def prepare_next_block(self, controller, S, size, time, Tend, MS, **kwargs):
313 """
314 If the no-memory-overhead version is used, we need to delete stuff that shouldn't be available. Otherwise, we
315 need to store all the stuff that we can.
317 Args:
318 controller (pySDC.Controller): The controller
319 S (pySDC.step): The current step
320 size (int): Number of ranks
321 time (float): The current time
322 Tend (float): The final time
323 MS (list): Active steps
325 Returns:
326 None
327 """
329 # delete values that should not be available in the next step
330 if self.params.no_storage:
331 self.prev.t = np.array([None] * self.params.n)
332 self.prev.dt = np.array([None] * self.params.n)
333 self.prev.u = [None] * self.params.n
334 self.prev.f = [None] * self.params.n
336 else:
337 # decide where we need to restart to store everything up to that point
338 restarts = [S.status.restart for S in MS]
339 restart_at = np.where(restarts)[0][0] if True in restarts else len(MS)
341 # store values in the current block that don't need restarting
342 if restart_at > S.status.slot:
343 self.store_values(S)
345 return None
347 def get_extrapolated_solution(self, S, **kwargs):
348 """
349 Combine values from previous steps to extrapolate.
351 Args:
352 S (pySDC.Step): The current step
354 Returns:
355 dtype_u: The extrapolated solution
356 """
357 if len(S.levels) > 1:
358 raise NotImplementedError("Extrapolated estimate only works on the finest level for now")
360 # prepare variables
361 u_ex = S.levels[0].u[-1] * 0.0
362 idx = np.argsort(self.prev.t)
364 # see if we have a solution for the current step already stored
365 if (abs(S.time + S.dt - self.prev.t) < 10.0 * np.finfo(float).eps).any():
366 idx_step = idx[np.argmin(abs(self.prev.t - S.time - S.dt))]
367 else:
368 idx_step = max(idx) + 1
370 # make a mask of all the steps we want to include in the extrapolation
371 mask = np.logical_and(idx < idx_step, idx >= idx_step - self.params.n)
373 # do the extrapolation by summing everything up
374 for i in range(self.params.n):
375 u_ex += self.coeff.u[i] * self.prev.u[idx[mask][i]] + self.coeff.f[i] * self.prev.f[idx[mask][i]]
377 return u_ex
379 def get_extrapolated_error(self, S, **kwargs):
380 """
381 The extrapolation estimate combines values of u and f from multiple steps to extrapolate and compare to the
382 solution obtained by the time marching scheme.
384 Args:
385 S (pySDC.Step): The current step
387 Returns:
388 None
389 """
390 u_ex = self.get_extrapolated_solution(S)
391 if u_ex is not None:
392 S.levels[0].status.error_extrapolation_estimate = abs(u_ex - S.levels[0].u[-1]) * self.coeff.prefactor
393 else:
394 S.levels[0].status.error_extrapolation_estimate = None
397class EstimateExtrapolationErrorWithinQ(EstimateExtrapolationErrorBase):
398 """
399 This convergence controller estimates the local error based on comparing the SDC solution to an extrapolated
400 solution within the quadrature matrix. Collocation methods compute a high order solution from a linear combination
401 of solutions at intermediate time points. While the intermediate solutions (a.k.a. stages) don't share the order of
402 accuracy with the solution at the end of the interval, for SDC we know that the order is equal to the number of
403 nodes + 1 (locally). This is because the solution to the collocation problem is a polynomial approximation of order
404 of the number of nodes.
405 That means we can do a Taylor expansion around the end point of the interval to higher order and after cancelling
406 terms just like we are used to with the extrapolation based error estimate across multiple steps, we get an error
407 estimate that is of the order of accuracy of the stages.
408 This can be used for adaptivity, for instance, with the nice property that it doesn't matter how we arrived at the
409 converged collocation solution, as long as we did. We don't rely on knowing the order of accuracy after every sweep,
410 only after convergence of the collocation problem has been achieved, which we can check from the residual.
411 """
413 def setup(self, controller, params, description, **kwargs):
414 """
415 We need this convergence controller to become active after the check for convergence, because we need the step
416 to be converged.
418 Args:
419 controller (pySDC.Controller): The controller
420 params (dict): The params passed for this specific convergence controller
421 description (dict): The description object used to instantiate the controller
423 Returns:
424 dict: Updated parameters with default values
425 """
426 from pySDC.implementations.convergence_controller_classes.check_convergence import CheckConvergence
428 num_nodes = description['sweeper_params']['num_nodes']
429 self.comm = description['sweeper_params'].get('comm', None)
430 if self.comm:
431 from mpi4py import MPI
433 self.sum = MPI.SUM
435 self.check_convergence = CheckConvergence.check_convergence
437 default_params = {
438 'Taylor_order': 2 * num_nodes,
439 'n': num_nodes,
440 'recompute_coefficients': False,
441 **params,
442 }
444 return {**super().setup(controller, params, description, **kwargs), **default_params}
446 def post_iteration_processing(self, controller, S, **kwargs):
447 """
448 Compute the extrapolated error estimate here if the step is converged.
450 Args:
451 controller (pySDC.Controller): The controller
452 S (pySDC.Step): The current step
454 Returns:
455 None
456 """
458 if not self.check_convergence(S):
459 return None
461 lvl = S.levels[0]
463 nodes_ = lvl.sweep.coll.nodes * S.dt
464 nodes = S.time + np.append(0, nodes_[:-1])
465 t_eval = S.time + nodes_[-1]
467 dts = np.append(nodes_[0], nodes_[1:] - nodes_[:-1])
468 self.params.Taylor_order = len(nodes)
469 self.params.n = len(nodes)
471 # compute the extrapolation coefficients
472 if None in self.coeff.u or self.params.recompute_coefficients:
473 self.get_extrapolation_coefficients(nodes, dts, t_eval)
475 # compute the extrapolated solution
476 if lvl.f[0] is None:
477 lvl.f[0] = lvl.prob.eval_f(lvl.u[0], lvl.time)
479 if type(lvl.f[0]) == imex_mesh:
480 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)]
481 elif type(lvl.f[0]) == mesh:
482 f = [lvl.f[i] if self.coeff.f[i] else 0.0 for i in range(len(lvl.f) - 1)]
483 else:
484 raise DataError(
485 f"Unable to store f from datatype {type(lvl.f[0])}, extrapolation based error estimate only\
486 works with types imex_mesh and mesh"
487 )
489 # compute the error with the weighted sum
490 if self.comm:
491 idx = (self.comm.rank + 1) % self.comm.size
492 sendbuf = self.coeff.u[idx] * lvl.u[idx] + self.coeff.f[idx] * f[idx]
493 u_ex = lvl.prob.dtype_u(lvl.prob.init, val=0.0) if self.comm.rank == self.comm.size - 1 else None
494 self.comm.Reduce(sendbuf, u_ex, op=self.sum, root=self.comm.size - 1)
495 else:
496 u_ex = lvl.prob.dtype_u(lvl.prob.init, val=0.0)
497 for i in range(self.params.n):
498 u_ex += self.coeff.u[i] * lvl.u[i] + self.coeff.f[i] * f[i]
500 # store the error
501 if self.comm:
502 error = (
503 abs(u_ex - lvl.u[self.comm.rank + 1]) * self.coeff.prefactor
504 if self.comm.rank == self.comm.size - 1
505 else None
506 )
507 lvl.status.error_extrapolation_estimate = self.comm.bcast(error, root=self.comm.size - 1)
508 else:
509 lvl.status.error_extrapolation_estimate = abs(u_ex - lvl.u[-1]) * self.coeff.prefactor
510 return None