Coverage for pySDC/implementations/transfer_classes/BaseTransferDelta.py: 99%
69 statements
« prev ^ index » next coverage.py v7.17.0, created at 2026-11-02 12:30 +0000
1r"""
2Space-time transfer for the delta-form hierarchy.
4:class:`~pySDC.core.base_transfer.BaseTransfer` makes a coarse level rebuild its own residual from
5:math:`\mathcal{O}(1)` coarse data and recovers the coarse-grid correction as :math:`u_G - u_G^{\mathrm{old}}`.
6Both have exact algebraic replacements -- :math:`r_G = R\,r_F`, which is what the FAS
7:math:`\tau` is *for*, and the correction is the sum of the sweep's own increments -- and this
8transfer uses those instead. It is what decides whether a level rebuilds its residual or is handed
9one: the delta-form sweepers do both, and this is what hands one down.
11The FAS :math:`\tau` is **substituted, not discarded**. Writing
12:math:`\tau = R(\Delta t\,Q_F f_F) - \Delta t\,Q_G f_G` and putting it into the coarse residual
13:math:`\Delta t (Q_G f_G) + u_G[0] - u_G[m] + \tau` cancels the :math:`\Delta t\,Q_G f_G` terms
14identically and leaves :math:`R\,r_F`. So a coarse level handed that residual is solving
15precisely the FAS-corrected problem -- and must *not* also add :math:`\tau`, which would count it
16twice. What is skipped is materialising :math:`\tau` as an array, since the quantity it exists to
17produce arrives directly.
19Identical to :class:`~pySDC.core.base_transfer.BaseTransfer` up to round-off at backend precision, verified
20for two, three and four levels and for PFASST -- and against a control with :math:`\tau` genuinely zeroed, which does
21not converge at all. It is also slightly cheaper, because nothing then reads :math:`\tau` except
22:meth:`compute_end_point` in the quadrature-update case: see :meth:`delta_transfer.restrict`.
24What the reformulation buys is that every coarse-level quantity becomes proportional to the fine
25residual rather than to :math:`|u|`. A coarse level is a preconditioner whose returned correction
26tends to zero, so once nothing on it is rebuilt out of :math:`\mathcal{O}(1)` state, its arithmetic
27error is proportional to that correction rather than to the solution -- which is what lets it run
28below backend precision without capping the accuracy of the iteration.
29"""
31from pySDC.core.base_transfer import BaseTransfer
32from pySDC.core.errors import UnlockError
35class delta_transfer(BaseTransfer):
36 """
37 Space-time transfer that passes a residual down and an accumulated correction up.
39 Drop-in for :class:`~pySDC.core.base_transfer.BaseTransfer`, and identical to it up to round-off when every
40 level is at backend precision. Falls back to the stock prolongation when the coarse level has banked no
41 corrections, so a hierarchy mixing delta-form and stock sweepers still runs.
42 """
44 def __init__(self, fine_level, coarse_level, *args, **kwargs):
45 super().__init__(fine_level, coarse_level, *args, **kwargs)
46 # the levels are connected fine to coarse, so the transfer onto the fine level can only
47 # find this one through the level's sweeper, and only once it is asked: see coarse_reads_tau
48 fine_level.sweep.coarser_transfer = self
50 def coarse_reads_tau(self):
51 """
52 Whether the coarse level's FAS :math:`\\tau` is read at all.
54 The sweep does not -- it is handed the restricted fine residual, which is what :math:`\\tau`
55 exists to produce -- and neither does :meth:`DeltaFormMixin.compute_residual`. The one
56 remaining reader on the level is :meth:`compute_end_point`, and only when the end point
57 comes from the quadrature update rather than a copy of the last node. The other is the
58 restriction to the next coarser level, which folds this :math:`\\tau` into its own when that
59 one is built, so the question covers every level below. When nothing reads it, building it
60 is a coarse ``integrate()`` and a restriction spent on a quantity that is then discarded.
62 Returns
63 -------
64 bool
65 Whether :math:`\\tau` has to be built.
66 """
67 SG = self.coarse.sweep
68 below = getattr(SG, 'coarser_transfer', None)
69 return not (SG.coll.right_is_node and not SG.params.do_coll_update) or (
70 below is not None and below.coarse_reads_tau()
71 )
73 def restrict_state(self):
74 """
75 Restrict ``u`` and re-evaluate ``f``, which is :meth:`BaseTransfer.restrict` without
76 :math:`\\tau`.
78 Returns
79 -------
80 None
82 Raises
83 ------
84 UnlockError
85 If the fine level has not been unlocked yet.
86 """
87 F, G = self.fine, self.coarse
88 SF, SG, PG = F.sweep, G.sweep, G.prob
89 if not F.status.unlocked:
90 raise UnlockError('fine level is still locked, cannot use data from there')
92 tmp_u = [self.space_transfer.restrict(F.u[m]) for m in range(1, SF.coll.num_nodes + 1)]
93 G.u[0] = self.space_transfer.restrict(F.u[0])
94 for n in range(1, SG.coll.num_nodes + 1):
95 # float(), not the raw np.float64 entry: see DeltaFormMixin._coeff
96 G.u[n] = float(self.Rcoll[n - 1, 0]) * tmp_u[0]
97 for m in range(1, SF.coll.num_nodes):
98 G.u[n] += float(self.Rcoll[n - 1, m]) * tmp_u[m]
100 G.f[0] = PG.eval_f(G.u[0], G.time)
101 for m in range(1, SG.coll.num_nodes + 1):
102 G.f[m] = PG.eval_f(G.u[m], G.time + G.dt * SG.coll.nodes[m - 1])
103 G.uold[m] = PG.dtype_u(G.u[m])
104 G.fold[m] = PG.dtype_f(G.f[m])
105 PG.use_restricted_boundary_values(G.u[1:], [G.time + G.dt * tau for tau in SG.coll.nodes])
107 G.status.unlocked = True
108 return None
110 def restrict(self):
111 """
112 Restrict the state, then hand the coarse level the restricted fine residual.
114 Returns
115 -------
116 None
117 """
118 SF, SG = self.fine.sweep, self.coarse.sweep
119 SF._delta_setup()
120 SG._delta_setup()
121 eps_F = SF._residual_nodes()
122 super().restrict() if self.coarse_reads_tau() else self.restrict_state()
124 tmp = [self.space_transfer.restrict(eps) for eps in eps_F]
125 eps_G = []
126 for n in range(SG.coll.num_nodes):
127 acc = float(self.Rcoll[n, 0]) * tmp[0]
128 for m in range(1, SF.coll.num_nodes):
129 acc += float(self.Rcoll[n, m]) * tmp[m]
130 eps_G.append(SG._to_work(self.coarse.prob, acc))
131 SG.eps_in = eps_G
132 SG.delta_acc = None
133 # after the rounding, so the reference is the value the level actually holds and a level
134 # that receives nothing shifts its residual by exactly zero
135 self.coarse.u0_reference = self.coarse.prob.dtype_u(self.coarse.u[0])
136 return None
138 def prolong(self):
139 """
140 Add the coarse level's accumulated correction to the fine level.
142 If the fine level is itself a coarse level of something above it, its residual and its own
143 accumulated correction are advanced by the same update, so neither has to be rebuilt from
144 :math:`\\mathcal{O}(1)` state later in the V-cycle.
146 Returns
147 -------
148 None
149 """
150 if self.coarse.sweep.delta_acc is None:
151 return super().prolong()
153 F, PF = self.fine, self.fine.prob
154 SF, SG = self.fine.sweep, self.coarse.sweep
155 SF._delta_setup()
156 tmp = [self.space_transfer.prolong(delta) for delta in self.coarse.sweep.delta_acc]
158 corr = []
159 for n in range(1, SF.coll.num_nodes + 1):
160 c = float(self.Pcoll[n - 1, 0]) * tmp[0]
161 for m in range(1, SG.coll.num_nodes):
162 c += float(self.Pcoll[n - 1, m]) * tmp[m]
163 # Quantise the correction at the correction precision, then bring it back to backend
164 # units before applying it. `_to_work` alone also divides by the sweep's scale, and
165 # adding *that* to the nodal value would be wrong by a factor of the scale.
166 c = SF._to_backend(PF, SF._to_work(PF, c))
167 corr.append(SF._to_work(PF, c))
169 t_node = F.time + F.dt * SF.coll.nodes[n - 1]
170 u_old, f_old = PF.dtype_u(F.u[n]), PF.dtype_f(F.f[n])
171 F.u[n] += c
172 F.f[n] = PF.eval_f(F.u[n], t_node)
173 if SF.eps_in is not None:
174 # records into SF._dfs, which is what advance_residual reads below
175 SF._f_increment(PF, F.f[n], f_old, u_old, c, t_node)
177 if SF.eps_in is not None:
178 SF.eps_in = SF.advance_residual(SF.eps_in, corr, SF._dfs)
179 SF.accumulate(corr)
180 return None