Coverage for pySDC/implementations/transfer_classes/BaseTransferDelta.py: 98%

1r""" 

2Space-time transfer for the delta-form hierarchy. 

3 

4:class:`BaseTransfer` makes a coarse level rebuild its own residual from :math:`\mathcal{O}(1)` 

5coarse data and recovers the coarse-grid correction as :math:`u_G - u_G^{\mathrm{old}}`. Both have 

6exact algebraic replacements -- :math:`\varepsilon_G = R\varepsilon_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. 

10 

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\varepsilon_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. 

18 

19Identical to :class:`BaseTransfer` up to round-off at backend precision, verified for two, three and 

20four 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`. 

23 

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""" 

30 

31from pySDC.core.base_transfer import BaseTransfer 

32from pySDC.core.errors import UnlockError 

33 

34 

35class delta_transfer(BaseTransfer): 

36 """ 

37 Space-time transfer that passes a residual down and an accumulated correction up. 

38 

39 Drop-in for :class:`BaseTransfer`, and identical to it up to round-off when every level is at 

40 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 """ 

43 

44 def coarse_reads_tau(self): 

45 """ 

46 Whether the coarse level reads the FAS :math:`\\tau` at all. 

47 

48 The sweep does not -- it is handed the restricted fine residual, which is what :math:`\\tau` 

49 exists to produce -- and neither does :meth:`DeltaFormMixin.compute_residual`. The one 

50 remaining reader is :meth:`compute_end_point`, and only when the end point comes from the 

51 quadrature update rather than a copy of the last node. When nothing reads it, building it is 

52 a coarse ``integrate()`` and a restriction spent on a quantity that is then discarded. 

53 

54 Returns 

55 ------- 

56 bool 

57 Whether :math:`\\tau` has to be built. 

58 """ 

59 SG = self.coarse.sweep 

60 return not (SG.coll.right_is_node and not SG.params.do_coll_update) 

61 

62 def restrict_state(self): 

63 """ 

64 Restrict ``u`` and re-evaluate ``f``, which is :meth:`BaseTransfer.restrict` without 

65 :math:`\\tau`. 

66 

67 Returns 

68 ------- 

69 None 

70 

71 Raises 

72 ------ 

73 UnlockError 

74 If the fine level has not been unlocked yet. 

75 """ 

76 F, G = self.fine, self.coarse 

77 SF, SG, PG = F.sweep, G.sweep, G.prob 

78 if not F.status.unlocked: 

79 raise UnlockError('fine level is still locked, cannot use data from there') 

80 

81 tmp_u = [self.space_transfer.restrict(F.u[m]) for m in range(1, SF.coll.num_nodes + 1)] 

82 G.u[0] = self.space_transfer.restrict(F.u[0]) 

83 for n in range(1, SG.coll.num_nodes + 1): 

84 # float(), not the raw np.float64 entry: see DeltaFormMixin._coeff 

85 G.u[n] = float(self.Rcoll[n - 1, 0]) * tmp_u[0] 

86 for m in range(1, SF.coll.num_nodes): 

87 G.u[n] += float(self.Rcoll[n - 1, m]) * tmp_u[m] 

88 

89 G.f[0] = PG.eval_f(G.u[0], G.time) 

90 for m in range(1, SG.coll.num_nodes + 1): 

91 G.f[m] = PG.eval_f(G.u[m], G.time + G.dt * SG.coll.nodes[m - 1]) 

92 G.uold[m] = PG.dtype_u(G.u[m]) 

93 G.fold[m] = PG.dtype_f(G.f[m]) 

94 

95 G.status.unlocked = True 

96 return None 

97 

98 def restrict(self): 

99 """ 

100 Restrict the state, then hand the coarse level the restricted fine residual. 

101 

102 Returns 

103 ------- 

104 None 

105 """ 

106 SF, SG = self.fine.sweep, self.coarse.sweep 

107 SF._delta_setup() 

108 SG._delta_setup() 

109 eps_F = SF._residual_nodes() 

110 super().restrict() if self.coarse_reads_tau() else self.restrict_state() 

111 

112 tmp = [self.space_transfer.restrict(eps) for eps in eps_F] 

113 eps_G = [] 

114 for n in range(SG.coll.num_nodes): 

115 acc = float(self.Rcoll[n, 0]) * tmp[0] 

116 for m in range(1, SF.coll.num_nodes): 

117 acc += float(self.Rcoll[n, m]) * tmp[m] 

118 eps_G.append(SG._to_work(self.coarse.prob, acc)) 

119 SG.eps_in = eps_G 

120 SG.delta_acc = None 

121 # after the rounding, so the reference is the value the level actually holds and a level 

122 # that receives nothing shifts its residual by exactly zero 

123 self.coarse.u0_reference = self.coarse.prob.dtype_u(self.coarse.u[0]) 

124 return None 

125 

126 def prolong(self): 

127 """ 

128 Add the coarse level's accumulated correction to the fine level. 

129 

130 If the fine level is itself a coarse level of something above it, its residual and its own 

131 accumulated correction are advanced by the same update, so neither has to be rebuilt from 

132 :math:`\\mathcal{O}(1)` state later in the V-cycle. 

133 

134 Returns 

135 ------- 

136 None 

137 """ 

138 if self.coarse.sweep.delta_acc is None: 

139 return super().prolong() 

140 

141 F, PF = self.fine, self.fine.prob 

142 SF, SG = self.fine.sweep, self.coarse.sweep 

143 SF._delta_setup() 

144 tmp = [self.space_transfer.prolong(delta) for delta in self.coarse.sweep.delta_acc] 

145 

146 corr = [] 

147 for n in range(1, SF.coll.num_nodes + 1): 

148 c = float(self.Pcoll[n - 1, 0]) * tmp[0] 

149 for m in range(1, SG.coll.num_nodes): 

150 c += float(self.Pcoll[n - 1, m]) * tmp[m] 

151 # Quantise the correction at the correction precision, then bring it back to backend 

152 # units before applying it. `_to_work` alone also divides by the sweep's scale, and 

153 # adding *that* to the nodal value would be wrong by a factor of the scale. 

154 c = SF._to_backend(PF, SF._to_work(PF, c)) 

155 corr.append(SF._to_work(PF, c)) 

156 

157 t_node = F.time + F.dt * SF.coll.nodes[n - 1] 

158 u_old, f_old = PF.dtype_u(F.u[n]), PF.dtype_f(F.f[n]) 

159 F.u[n] += c 

160 F.f[n] = PF.eval_f(F.u[n], t_node) 

161 if SF.eps_in is not None: 

162 # records into SF._dfs, which is what advance_residual reads below 

163 SF._f_increment(PF, F.f[n], f_old, u_old, c, t_node) 

164 

165 if SF.eps_in is not None: 

166 SF.eps_in = SF.advance_residual(SF.eps_in, corr, SF._dfs) 

167 SF.accumulate(corr) 

168 return None