Coverage for pySDC/implementations/sweeper_classes/delta_form.py: 98%
189 statements
« prev ^ index » next coverage.py v7.16.1, created at 2026-09-25 20:28 +0000
1r"""
2Deferred-correction ("delta-form") SDC sweepers.
4A standard SDC sweep
6.. math::
7 u^{k+1}_m = u_0 + \tau_m + \Delta t (Q f^k)_m
8 + \Delta t \sum_j Q^\Delta_{mj}\,(f^{k+1}_j - f^k_j)
10is algebraically identical to, with :math:`\delta_m = u^{k+1}_m - u^k_m` and the collocation
11residual :math:`\varepsilon_m = u_0 + \tau_m + \Delta t (Q f^k)_m - u^k_m`,
13.. math::
14 \delta_m = \varepsilon_m + \Delta t \sum_j Q^\Delta_{mj}\,\Delta f_j,
15 \qquad \Delta f_j = f(u^k_j + \delta_j) - f(u^k_j).
17Written this way, every sweep is iterative refinement: a high-precision residual, a correction
18solve, and a high-precision update ``u <- u + delta``. No Jacobian appears, so an IMEX splitting
19survives unchanged.
21The point of the reformulation is that the quantity handed to the node-local solver is a
22*correction*. Its magnitude tends to zero as the sweeps converge, so a reduced-precision solve
23introduces an error proportional to :math:`|\delta|` rather than to :math:`|u|` and therefore does
24not cap the attainable accuracy.
26Three node-local strategies are supported, selected automatically:
28``solve_system_delta``
29 Used when the problem provides it. Solves
30 :math:`\delta - \alpha[f(w+\delta) - f(w)] = r` for the correction directly. This is the only
31 option for a nonlinear implicit operator, and the only one that hands a reduced-precision
32 solver a small unknown.
34``linear_implicit=True``
35 For a linear or affine implicit operator, :math:`f(w+\delta) - f(w) = A\delta`, so the stock
36 ``solve_system`` already solves the correction equation once the affine part
37 :math:`\alpha f(0, t)` is removed from the right-hand side. No problem class needs changing.
39fallback
40 Otherwise the substitution :math:`y = u^k_m + \delta_m` reduces the correction equation to the
41 ordinary implicit solve. Always correct and identical to :class:`generic_implicit`, but the
42 solver sees an :math:`\mathcal{O}(1)` unknown, so there is no precision benefit.
44``correction_precision`` additionally stores the small quantities
45(:math:`\varepsilon`, :math:`\delta`, :math:`\Delta f`) in a reduced-precision datatype built from
46the problem's own ``init`` tuple, scaled by the residual's magnitude so the format's mantissa is
47used however small the correction gets.
49The increment :math:`\Delta f` is formed by subtracting two stored right-hand sides unless the
50problem provides ``eval_f_increment(base, delta, t)``, which expands it analytically. That
51subtraction is a cancellation carrying the operator norm, so it binds as soon as a level runs below
52backend precision.
54The same sweeper serves any number of levels, and which it is doing is decided by the transfer, not
55by the choice of sweeper. On the finest level it computes its own residual, which is where the
56accuracy of the whole iteration is set. Given a transfer that hands one down -- by setting
57``eps_in`` on the level's sweeper -- it uses that instead and advances it in place, so nothing is
58ever rebuilt out of :math:`\mathcal{O}(1)` coarse state. :class:`BaseTransfer` hands nothing down,
59so on a stock hierarchy this sweeper reproduces MLSDC and PFASST exactly.
60"""
62import numpy as np
64from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit
65from pySDC.implementations.sweeper_classes.imex_1st_order import imex_1st_order
68class DeltaFormMixin:
69 """Shared machinery for the delta-form sweepers."""
71 def _delta_setup(self):
72 """Read the optional sweeper parameters, and arm the per-sweep recorders."""
73 token = getattr(self.params, 'correction_precision', None)
74 self._work_dtype = None if token is None else np.dtype(token)
75 self._linear_implicit = bool(getattr(self.params, 'linear_implicit', False))
76 self._deltas, self._dfs = [], []
78 eps_in = None
79 """Residual handed down by the transfer, or ``None`` on the finest level."""
81 delta_acc = None
82 """Corrections this level has accumulated since the last restriction."""
84 def sync_initial_value(self):
85 r"""
86 Follow a change of :math:`u_0` made after the residual was handed down.
88 PFASST receives the initial value from the predecessor *after* the restriction, directly
89 into ``u[0]``. The residual depends on it additively, so following it is one addition:
90 :math:`\varepsilon_m \leftarrow \varepsilon_m + (u_0 - u_0^{\mathrm{ref}})`. Exactly zero
91 when nothing arrived, which is every serial run.
93 This is the one place the hierarchy still differences two :math:`\mathcal{O}(1)` values, and
94 it is the reason a reduced-precision coarse level buys less under PFASST than under MLSDC:
95 the difference is small -- it is the coarse-versus-fine discrepancy at the step interface,
96 and it converges to zero -- but it is *formed* by cancelling two values of size
97 :math:`|u|`. Putting the step-to-step exchange itself in delta form is what would remove it.
99 Returns
100 -------
101 None
102 """
103 lvl = self.level
104 if self.eps_in is None or lvl.u0_reference is None:
105 return
106 shift = lvl.u[0] - lvl.u0_reference
107 self.eps_in = [eps + shift for eps in self.eps_in]
108 lvl.u0_reference = lvl.prob.dtype_u(lvl.u[0])
109 return None
111 def compute_residual(self, stage=''):
112 r"""
113 Report the residual this level is already tracking, instead of rebuilding one.
115 A level with an inherited residual carries it forward through every sweep and prolongation,
116 so recomputing it from :math:`\mathcal{O}(1)` state would be both redundant and less
117 accurate. It is also what makes the :math:`\tau` term unnecessary: the only other reader is
118 :meth:`compute_end_point`, and only when the end point comes from the quadrature update.
120 Parameters
121 ----------
122 stage : str
123 The stage of the step this level belongs to.
125 Returns
126 -------
127 None
128 """
129 if self.eps_in is None:
130 return super().compute_residual(stage=stage)
132 lvl = self.level
133 if stage in self.params.skip_residual_computation:
134 lvl.status.residual = 0.0 if lvl.status.residual is None else lvl.status.residual
135 return None
137 lvl.residual = [lvl.prob.dtype_u(eps) for eps in self.eps_in]
138 norms = [abs(eps) for eps in self.eps_in]
139 kind = lvl.params.residual_type
140 if kind not in ('full_abs', 'last_abs', 'full_rel', 'last_rel'):
141 raise NotImplementedError(f'residual type "{kind}" not implemented!')
142 value = norms[-1] if kind.startswith('last') else max(norms)
143 lvl.status.residual = value / abs(lvl.u[0]) if kind.endswith('rel') else value
144 lvl.status.updated = False
145 return None
147 def advance_residual(self, eps, deltas, dfs):
148 r"""
149 Advance a residual by an update: :math:`\varepsilon \leftarrow \varepsilon - \delta
150 + \Delta t (Q \Delta f)`.
152 Every term is small, so this never cancels. It is exact for any update to the nodal values,
153 which is why it serves both a sweep and a prolongation.
155 Parameters
156 ----------
157 eps : list
158 The residual, one entry per node.
159 deltas : list
160 The update applied to the nodal values, one entry per node.
161 dfs : list
162 The resulting right-hand side increments, one entry per node.
164 Returns
165 -------
166 list
167 The advanced residual.
168 """
169 dt, Q = self.level.dt, self.coll.Qmat
170 out = []
171 for m in range(len(eps)):
172 acc = eps[m] - deltas[m]
173 for j in range(len(dfs)):
174 if Q[m + 1, j + 1] != 0.0:
175 acc += self._coeff(dt * Q[m + 1, j + 1]) * dfs[j]
176 out.append(acc)
177 return out
179 def accumulate(self, deltas):
180 """Add one round of corrections to what this level owes upwards."""
181 self.delta_acc = (
182 deltas if self.delta_acc is None else [a + d for a, d in zip(self.delta_acc, deltas, strict=True)]
183 )
185 def update_nodes(self):
186 r"""
187 Sweep, and keep the level's bookkeeping straight if it is part of a hierarchy.
189 The sweep itself is :meth:`_sweep_nodes`, which each concrete sweeper supplies. Around it:
190 follow any change of :math:`u_0` that arrived after the residual was handed down, advance
191 that residual by the update just applied, and bank the corrections for a transfer to prolong.
193 A level that computes its own residual has no bookkeeping to do, so on a single-level run
194 every line below the sweep is a no-op. That is why there is one sweeper rather than two.
196 Returns
197 -------
198 None
199 """
200 self.sync_initial_value()
201 self._sweep_nodes()
202 if self.eps_in is None:
203 return None
205 deltas = [self._to_work(self.level.prob, d) for d in self._deltas]
206 self.accumulate(deltas)
207 self.eps_in = self.advance_residual(self.eps_in, deltas, self._dfs)
208 return None
210 _work_scale = 1.0
211 r"""Shared divisor applied to the correction quantities before they are stored."""
213 def _work_init(self, prob):
214 """Build the problem's ``init`` tuple with the correction dtype substituted."""
215 return (prob.init[0], prob.init[1], self._work_dtype)
217 def _scales_corrections(self):
218 """
219 Whether this sweeper may choose a scale for its correction quantities.
221 A level that computes its own residual may: everything it stores is derived from that
222 residual within the same sweep, so one divisor keeps them all commensurate and ordinary
223 arithmetic on them stays correct. A level that *inherits* a residual may not, because the
224 inherited value was scaled by whoever produced it and is carried across sweeps.
225 """
226 return self.eps_in is None
228 def _set_work_scale(self, values):
229 r"""
230 Choose the divisor for this sweep, from the residual the corrections will be built out of.
232 This is what lets a correction be stored below ``float16``'s smallest normal, 6.1e-5. The
233 delta form drives :math:`\varepsilon` and :math:`\delta` towards zero on purpose, and half
234 precision has almost no mantissa left down there -- 1.3e-2 relative at 1e-6, 1.9e-1 at 1e-7 --
235 so an unscaled correction turns to noise exactly when it starts to matter. Dividing by the
236 residual's own magnitude keeps the stored values at :math:`\mathcal{O}(1)`, which is block
237 floating point, and is what half-precision hardware does anyway.
239 Parameters
240 ----------
241 values : list
242 The residual at the collocation nodes, in backend units.
243 """
244 if self._work_dtype is None or not self._scales_corrections():
245 self._work_scale = 1.0
246 return
247 biggest = max((abs(value) for value in values), default=0.0)
248 self._work_scale = float(biggest) if biggest > 0.0 else 1.0
250 def _to_work(self, prob, value):
251 """
252 Store a small correction quantity in a reduced-precision datatype.
254 Returns the value unchanged when no reduced precision was requested, so the default path
255 makes no assumption about the datatype and works with any pySDC backend.
257 Raises
258 ------
259 NotImplementedError
260 If ``correction_precision`` was requested but the datatype cannot be built at another
261 precision, as is the case for datatypes not backed by a numpy array.
262 """
263 if self._work_dtype is None:
264 return value
265 try:
266 me = prob.dtype_u(self._work_init(prob))
267 me[:] = value if self._work_scale == 1.0 else value / self._work_scale
268 except (TypeError, NotImplementedError) as error:
269 raise NotImplementedError(
270 f'correction_precision is not supported for {prob.dtype_u.__name__}: it cannot be '
271 f'built at a different precision from the problem init tuple ({error})'
272 ) from error
273 return me
275 def _to_backend(self, prob, value):
276 """
277 Lift a possibly reduced-precision quantity back to backend precision.
279 A no-op when no reduced precision is in play, which keeps the default path datatype-agnostic.
280 """
281 if self._work_dtype is None:
282 return value
283 me = prob.dtype_u(prob.init)
284 me[:] = value
285 if self._work_scale != 1.0:
286 # widen first, scale second. The other order multiplies at the reduced precision, and a
287 # scale of 1e-8 then lands the result below float16's smallest subnormal on the way out
288 # -- the value is destroyed before it ever reaches the backend-precision array.
289 me *= self._work_scale
290 return me
292 def _coeff(self, value):
293 """
294 Cast a scalar coefficient so an accumulation stays at the precision it should.
296 The ``float`` matters. A coefficient taken out of a numpy array is an ``np.float64``, and
297 multiplying a reduced-precision level quantity by one of those upcasts the result to
298 ``float64`` under NEP 50 -- NumPy 2's rule -- which would quietly put the whole level back at
299 backend precision. A plain Python float is weak under both the old and the new rule, so the
300 array's own dtype wins.
301 """
302 if self._work_dtype is None:
303 return float(value)
304 return self._work_dtype.type(value)
306 def _residual_nodes(self):
307 r"""
308 Compute :math:`\varepsilon_m = u_0 + \tau_m + \Delta t (Q f^k)_m - u^k_m`.
310 This is the high-precision residual of iterative refinement. It is a difference of
311 :math:`\mathcal{O}(1)` quantities and is therefore always formed in backend precision --
312 unless a transfer handed one down, in which case that one is already exact and small, and
313 rebuilding it here is what the delta-form hierarchy exists to avoid.
315 Returns
316 -------
317 list
318 One residual per collocation node.
319 """
320 if self.eps_in is not None:
321 return self.eps_in
323 lvl = self.level
324 eps = self.integrate()
325 for m in range(self.coll.num_nodes):
326 eps[m] += lvl.u[0]
327 eps[m] -= lvl.u[m + 1]
328 if lvl.tau[m] is not None:
329 eps[m] += lvl.tau[m]
330 return eps
332 def _f_increment(self, prob, f_new, f_old, u_old, delta, t_node):
333 r"""
334 Form the right-hand side increment :math:`\Delta f = f(w+\delta) - f(w)`.
336 Formed by subtraction unless the problem can expand it analytically. The subtraction is a
337 cancellation of two :math:`\mathcal{O}(|f|)` quantities, so its absolute error is
338 :math:`\varepsilon |f|` in whatever precision the level stores ``f`` at, and :math:`|f|`
339 carries the operator norm. That is harmless while the level is at backend precision and
340 becomes the binding term as soon as it is not, which is why a problem living on a
341 reduced-precision level should provide ``eval_f_increment``.
343 Parameters
344 ----------
345 prob : pySDC.core.problem.Problem
346 The problem on this level.
347 f_new, f_old : dtype_f
348 ``f`` at :math:`w+\delta` and at :math:`w`, used by the subtraction fallback.
349 u_old : dtype_u
350 The base state :math:`w`.
351 delta : dtype_u
352 The correction :math:`\delta`.
353 t_node : float
354 Physical time of the collocation node.
356 Returns
357 -------
358 dtype_f
359 The increment, with the same splitting as ``eval_f``.
360 """
361 if hasattr(prob, 'eval_f_increment'):
362 increment = prob.eval_f_increment(u_old, delta, t_node)
363 else:
364 increment = prob.dtype_f(f_new)
365 increment -= f_old
366 # recorded for the residual recursion, so a transfer never recovers it by subtraction
367 self._dfs.append(total_increment(prob, increment))
368 return increment
370 def _solve_correction(self, rhs_corr, alpha, u_old, f_old, t_node, implicit_part=None):
371 r"""
372 Solve the node-local correction equation.
374 Parameters
375 ----------
376 rhs_corr : dtype_u
377 Right-hand side :math:`r` of the correction equation.
378 alpha : float
379 Implicit prefactor :math:`\alpha = \Delta t Q^\Delta_{mm}`.
380 u_old : dtype_u
381 Current nodal value :math:`u^k_m`, the base state of the correction.
382 f_old : dtype_f
383 ``f`` evaluated at ``u_old``; already stored on the level, so it costs nothing.
384 t_node : float
385 Physical time of the collocation node.
386 implicit_part : dtype_u, optional
387 The implicit component of ``f_old`` for IMEX problems. Defaults to ``f_old``.
389 Returns
390 -------
391 dtype_u
392 The correction :math:`\delta_m`.
393 """
394 prob = self.level.prob
395 f_impl_old = f_old if implicit_part is None else implicit_part
397 rhs_phys = self._to_backend(prob, rhs_corr)
399 if alpha == 0:
400 # explicit node: the correction is the residual itself
401 delta = rhs_phys
402 elif hasattr(prob, 'solve_system_delta'):
403 delta = prob.solve_system_delta(rhs_phys, alpha, u_old, f_old, t_node)
404 elif self._linear_implicit:
405 # f(w+d) - f(w) = A d, so solve_system already solves the correction equation once the
406 # affine part f(0, t) has been removed. f(0, t) vanishes for a homogeneous operator.
407 zero = prob.dtype_u(prob.init, val=0.0)
408 affine = prob.eval_f(zero, t_node)
409 rhs_phys -= alpha * (affine if implicit_part is None else affine.impl)
410 delta = prob.solve_system(rhs_phys, alpha, zero, t_node)
411 else:
412 # Fallback: substitute y = u_old + delta. Always correct, but the solver sees an O(1)
413 # unknown, so there is no precision benefit.
414 rhs_phys += u_old
415 rhs_phys -= alpha * f_impl_old
416 solution = prob.solve_system(rhs_phys, alpha, u_old, t_node)
417 delta = prob.dtype_u(solution)
418 delta -= u_old
419 # recorded so a transfer never has to recover the correction by subtraction
420 self._deltas.append(delta)
421 return delta
424class delta_implicit(DeltaFormMixin, generic_implicit):
425 """
426 Delta-form counterpart of :class:`generic_implicit`.
428 Mathematically identical to the standard sweep; see the module docstring for the sweeper
429 parameters ``correction_precision`` and ``linear_implicit``.
430 """
432 def _sweep_nodes(self):
433 """
434 Perform one delta-form sweep over all collocation nodes.
436 Returns
437 -------
438 None
439 """
440 lvl = self.level
441 prob = lvl.prob
442 assert lvl.status.unlocked
443 num_nodes = self.coll.num_nodes
444 self._delta_setup()
446 residual = self._residual_nodes()
447 self._set_work_scale(residual)
448 eps = [self._to_work(prob, value) for value in residual]
449 df = [None] * (num_nodes + 1)
451 for m in range(num_nodes):
452 t_node = lvl.time + lvl.dt * self.coll.nodes[m]
454 rhs_corr = type(eps[m])(eps[m])
455 for j in range(1, m + 1):
456 if self.QI[m + 1, j] != 0.0:
457 rhs_corr += self._coeff(lvl.dt * self.QI[m + 1, j]) * df[j]
459 alpha = lvl.dt * self.QI[m + 1, m + 1]
460 u_old = prob.dtype_u(lvl.u[m + 1])
461 f_old = prob.dtype_f(lvl.f[m + 1])
463 delta = self._solve_correction(rhs_corr, alpha, u_old, f_old, t_node)
465 lvl.u[m + 1] = u_old + self._to_backend(prob, self._to_work(prob, delta))
466 lvl.f[m + 1] = prob.eval_f(lvl.u[m + 1], t_node)
468 increment = self._f_increment(prob, lvl.f[m + 1], f_old, u_old, delta, t_node)
469 df[m + 1] = self._to_work(prob, increment)
471 lvl.status.updated = True
472 return None
475class delta_imex_1st_order(DeltaFormMixin, imex_1st_order):
476 """
477 Delta-form counterpart of :class:`imex_1st_order`.
479 The correction equation contains only differences of ``f``, never a Jacobian, so the
480 explicit/implicit splitting is untouched.
481 """
483 def _sweep_nodes(self):
484 """
485 Perform one delta-form IMEX sweep over all collocation nodes.
487 ``QE`` is strictly lower triangular, which :class:`imex_1st_order` already enforces, so the
488 explicit part never contributes to the node-local solve.
490 Returns
491 -------
492 None
493 """
494 lvl = self.level
495 prob = lvl.prob
496 assert lvl.status.unlocked
497 num_nodes = self.coll.num_nodes
498 self._delta_setup()
500 residual = self._residual_nodes()
501 self._set_work_scale(residual)
502 eps = [self._to_work(prob, value) for value in residual]
503 df_impl = [None] * (num_nodes + 1)
504 df_expl = [None] * (num_nodes + 1)
506 for m in range(num_nodes):
507 t_node = lvl.time + lvl.dt * self.coll.nodes[m]
509 rhs_corr = type(eps[m])(eps[m])
510 for j in range(1, m + 1):
511 if self.QI[m + 1, j] != 0.0:
512 rhs_corr += self._coeff(lvl.dt * self.QI[m + 1, j]) * df_impl[j]
513 if self.QE[m + 1, j] != 0.0:
514 rhs_corr += self._coeff(lvl.dt * self.QE[m + 1, j]) * df_expl[j]
516 alpha = lvl.dt * self.QI[m + 1, m + 1]
517 u_old = prob.dtype_u(lvl.u[m + 1])
518 f_old = prob.dtype_f(lvl.f[m + 1])
520 delta = self._solve_correction(rhs_corr, alpha, u_old, f_old, t_node, implicit_part=f_old.impl)
522 lvl.u[m + 1] = u_old + self._to_backend(prob, self._to_work(prob, delta))
523 lvl.f[m + 1] = prob.eval_f(lvl.u[m + 1], t_node)
525 increment = self._f_increment(prob, lvl.f[m + 1], f_old, u_old, delta, t_node)
526 df_impl[m + 1] = self._to_work(prob, prob.dtype_u(increment.impl))
527 df_expl[m + 1] = self._to_work(prob, prob.dtype_u(increment.expl))
529 lvl.status.updated = True
530 return None
533def total_increment(prob, increment):
534 """
535 The full right-hand side increment, recombining an IMEX splitting.
537 Parameters
538 ----------
539 prob : pySDC.core.problem.Problem
540 The problem the increment belongs to.
541 increment : dtype_f
542 The increment, possibly split into ``impl`` and ``expl``.
544 Returns
545 -------
546 dtype_u
547 The sum of the parts.
548 """
549 if not hasattr(increment, 'impl'):
550 return increment
551 return prob.dtype_u(increment.impl) + increment.expl