Coverage for pySDC/implementations/sweeper_classes/Runge_Kutta_Nystrom.py: 96%

1import numpy as np 

2import logging 

3 

4from pySDC.core.sweeper import Sweeper, _Pars 

5from pySDC.core.errors import ParameterError 

6from pySDC.implementations.datatype_classes.particles import particles, fields, acceleration 

7from pySDC.implementations.sweeper_classes.Runge_Kutta import RungeKutta 

8 

9 

10class ButcherTableauNoCollUpdate(object): 

11 """Version of Butcher Tableau that does not need a collocation update because the weights are put in the last line of Q""" 

12 

13 def __init__(self, weights, nodes, matrix): 

14 """ 

15 Initialization routine to get a quadrature matrix out of a Butcher tableau 

16 

17 Args: 

18 weights (numpy.ndarray): Butcher tableau weights 

19 nodes (numpy.ndarray): Butcher tableau nodes 

20 matrix (numpy.ndarray): Butcher tableau entries 

21 """ 

22 # check if the arguments have the correct form 

23 if type(matrix) != np.ndarray: 

24 raise ParameterError('Runge-Kutta matrix needs to be supplied as a numpy array!') 

25 elif len(np.unique(matrix.shape)) != 1 or len(matrix.shape) != 2: 

26 raise ParameterError('Runge-Kutta matrix needs to be a square 2D numpy array!') 

27 

28 if type(weights) != np.ndarray: 

29 raise ParameterError('Weights need to be supplied as a numpy array!') 

30 elif len(weights.shape) != 1: 

31 raise ParameterError(f'Incompatible dimension of weights! Need 1, got {len(weights.shape)}') 

32 elif len(weights) != matrix.shape[0]: 

33 raise ParameterError(f'Incompatible number of weights! Need {matrix.shape[0]}, got {len(weights)}') 

34 

35 if type(nodes) != np.ndarray: 

36 raise ParameterError('Nodes need to be supplied as a numpy array!') 

37 elif len(nodes.shape) != 1: 

38 raise ParameterError(f'Incompatible dimension of nodes! Need 1, got {len(nodes.shape)}') 

39 elif len(nodes) != matrix.shape[0]: 

40 raise ParameterError(f'Incompatible number of nodes! Need {matrix.shape[0]}, got {len(nodes)}') 

41 

42 self.globally_stiffly_accurate = np.allclose(matrix[-1], weights) 

43 

44 self.tleft = 0.0 

45 self.tright = 1.0 

46 self.num_solution_stages = 0 if self.globally_stiffly_accurate else 1 

47 self.num_nodes = matrix.shape[0] + self.num_solution_stages 

48 self.weights = weights 

49 

50 if self.globally_stiffly_accurate: 

51 # For globally stiffly accurate methods, the last row of the Butcher tableau is the same as the weights. 

52 self.nodes = np.append([0], nodes) 

53 self.Qmat = np.zeros([self.num_nodes + 1, self.num_nodes + 1]) 

54 self.Qmat[1:, 1:] = matrix 

55 else: 

56 self.nodes = np.append(np.append([0], nodes), [1]) 

57 self.Qmat = np.zeros([self.num_nodes + 1, self.num_nodes + 1]) 

58 self.Qmat[1:-1, 1:-1] = matrix 

59 self.Qmat[-1, 1:-1] = weights # this is for computing the solution to the step from the previous stages 

60 

61 self.left_is_node = True 

62 self.right_is_node = self.nodes[-1] == self.tright 

63 

64 # compute distances between the nodes 

65 if self.num_nodes > 1: 

66 self.delta_m = self.nodes[1:] - self.nodes[:-1] 

67 else: 

68 self.delta_m = np.zeros(1) 

69 self.delta_m[0] = self.nodes[0] - self.tleft 

70 

71 # check if the RK scheme is implicit 

72 self.implicit = any(matrix[i, i] != 0 for i in range(self.num_nodes - self.num_solution_stages)) 

73 

74 

75class RungeKuttaNystrom(RungeKutta): 

76 """ 

77 Runge-Kutta scheme that fits the interface of a sweeper. 

78 Actually, the sweeper idea fits the Runge-Kutta idea when using only lower triangular rules, where solutions 

79 at the nodes are successively computed from earlier nodes. However, we only perform a single iteration of this. 

80 

81 We have two choices to realise a Runge-Kutta sweeper: We can choose Q = Q_Delta = <Butcher tableau>, but in this 

82 implementation, that would lead to a lot of wasted FLOPS from integrating with Q and then with Q_Delta and 

83 subtracting the two. For that reason, we built this new sweeper, which does not have a preconditioner. 

84 

85 This class only supports lower triangular Butcher tableaux such that the system can be solved with forward 

86 substitution. In this way, we don't get the maximum order that we could for the number of stages, but computing the 

87 stages is much cheaper. In particular, if the Butcher tableaux is strictly lower triangular, we get an explicit 

88 method, which does not require us to solve a system of equations to compute the stages. 

89 

90 Please be aware that all fundamental parameters of the Sweeper are ignored. These include 

91 

92 - num_nodes 

93 - collocation_class 

94 - initial_guess 

95 - QI 

96 

97 All of these variables are either determined by the RK rule, or are not part of an RK scheme. 

98 

99 Attribues: 

100 butcher_tableau (ButcherTableauNoCollUpdate): Butcher tableau for the Runge-Kutta scheme that you want 

101 """ 

102 

103 ButcherTableauClass = ButcherTableauNoCollUpdate 

104 weights_bar = None 

105 matrix_bar = None 

106 

107 def __init__(self, params): 

108 """ 

109 Initialization routine for the custom sweeper 

110 

111 Args: 

112 params: parameters for the sweeper 

113 """ 

114 super().__init__(params) 

115 self.coll_bar = self.get_Butcher_tableau_bar() 

116 self.Qx = self.coll_bar.Qmat 

117 

118 @classmethod 

119 def get_Butcher_tableau_bar(cls): 

120 return cls.ButcherTableauClass(cls.weights_bar, cls.nodes, cls.matrix_bar) 

121 

122 def get_full_f(self, f): 

123 """ 

124 Test the right hand side function is the correct type 

125 

126 Args: 

127 f (dtype_f): Right hand side at a single node 

128 

129 Returns: 

130 mesh: Full right hand side as a mesh 

131 """ 

132 

133 if type(f) in [particles, fields, acceleration]: 

134 return f 

135 else: 

136 raise NotImplementedError(f'Type \"{type(f)}\" not implemented') 

137 

138 def update_nodes(self): 

139 """ 

140 Update the u- and f-values at the collocation nodes 

141 

142 Returns: 

143 None 

144 """ 

145 

146 # get current level and problem description 

147 L = self.level 

148 P = L.prob 

149 

150 # only if the level has been touched before 

151 assert L.status.unlocked 

152 assert L.status.sweep <= 1, "RK schemes are direct solvers. Please perform only 1 iteration!" 

153 

154 # get number of collocation nodes for easier access 

155 M = self.coll.num_nodes 

156 

157 for m in range(0, M): 

158 # build rhs, consisting of the known values from above and new values from previous nodes (at k+1) 

159 rhs = P.dtype_u(L.u[0]) 

160 rhs.pos += L.dt * self.coll.nodes[m + 1] * L.u[0].vel 

161 

162 for j in range(1, m + 1): 

163 # build RHS from f-terms (containing the E field) and the B field 

164 

165 f = P.build_f(L.f[j], L.u[j], L.time + L.dt * self.coll.nodes[j]) 

166 rhs.pos += L.dt**2 * self.Qx[m + 1, j] * self.get_full_f(f) 

167 """ 

168 

169 Implicit part only works for Velocity-Verlet scheme 

170 Boris solver for the implicit part 

171 

172 """ 

173 

174 if self.coll.implicit: 

175 ck = rhs.vel * 0.0 

176 L.f[3] = P.eval_f(rhs, L.time + L.dt) 

177 rhs.vel = P.boris_solver(ck, L.dt, L.f[0], L.f[3], L.u[0]) 

178 

179 else: 

180 rhs.vel += L.dt * self.QI[m + 1, j] * self.get_full_f(f) 

181 

182 # implicit solve with prefactor stemming from the diagonal of Qd 

183 L.u[m + 1] = rhs 

184 # update function values 

185 if self.coll.implicit: 

186 # That is why it only works for the Velocity-Verlet scheme 

187 L.f[0] = P.eval_f(L.u[0], L.time) 

188 L.f[m + 1] = P.dtype_f(L.f[0]) 

189 else: 

190 if m != self.coll.num_nodes - 1: 

191 L.f[m + 1] = P.eval_f(L.u[m + 1], L.time + L.dt * self.coll.nodes[m]) 

192 

193 # indicate presence of new values at this level 

194 

195 L.status.updated = True 

196 

197 return None 

198 

199 def compute_end_point(self): 

200 """ 

201 In this Runge-Kutta implementation, the solution to the step is always stored in the last node 

202 """ 

203 self.level.uend = self.level.u[-1] 

204 

205 

206class RKN(RungeKuttaNystrom): 

207 """ 

208 Runge-Kutta-Nystrom method 

209 https://link.springer.com/book/10.1007/978-3-540-78862-1 

210 page: 284 

211 Chapter: II.14 Numerical methods for Second order differential equations 

212 """ 

213 

214 nodes = np.array([0.0, 0.5, 0.5, 1]) 

215 weights = np.array([1.0, 2.0, 2.0, 1.0]) / 6.0 

216 matrix = np.zeros([4, 4]) 

217 matrix[1, 0] = 0.5 

218 matrix[2, 1] = 0.5 

219 matrix[3, 2] = 1.0 

220 

221 weights_bar = np.array([1.0, 1.0, 1.0, 0]) / 6.0 

222 matrix_bar = np.zeros([4, 4]) 

223 matrix_bar[1, 0] = 1 / 8 

224 matrix_bar[2, 0] = 1 / 8 

225 matrix_bar[3, 2] = 1 / 2 

226 

227 

228class Velocity_Verlet(RungeKuttaNystrom): 

229 """ 

230 Velocity-Verlet scheme 

231 https://de.wikipedia.org/wiki/Verlet-Algorithmus 

232 """ 

233 

234 nodes = np.array([1.0, 1.0]) 

235 weights = np.array([1 / 2, 0]) 

236 matrix = np.zeros([2, 2]) 

237 matrix[1, 1] = 1 

238 weights_bar = np.array([1 / 2, 0]) 

239 matrix_bar = np.zeros([2, 2])