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

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 ButcherTableau 

8from copy import deepcopy 

9from pySDC.implementations.sweeper_classes.Runge_Kutta import RungeKutta 

10 

11 

12class RungeKuttaNystrom(RungeKutta): 

13 """ 

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

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

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

17 

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

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

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

21 

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

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

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

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

26 

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

28 

29 - num_nodes 

30 - collocation_class 

31 - initial_guess 

32 - QI 

33 

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

35 

36 Attribues: 

37 butcher_tableau (ButcherTableau): Butcher tableau for the Runge-Kutta scheme that you want 

38 """ 

39 

40 def __init__(self, params): 

41 """ 

42 Initialization routine for the custom sweeper 

43 

44 Args: 

45 params: parameters for the sweeper 

46 """ 

47 # set up logger 

48 self.logger = logging.getLogger('sweeper') 

49 

50 essential_keys = ['butcher_tableau'] 

51 for key in essential_keys: 

52 if key not in params: 

53 msg = 'need %s to instantiate step, only got %s' % (key, str(params.keys())) 

54 self.logger.error(msg) 

55 raise ParameterError(msg) 

56 

57 # check if some parameters are set which only apply to actual sweepers 

58 for key in ['initial_guess', 'collocation_class', 'num_nodes']: 

59 if key in params: 

60 self.logger.warning(f'"{key}" will be ignored by Runge-Kutta sweeper') 

61 

62 # set parameters to their actual values 

63 params['initial_guess'] = 'zero' 

64 params['collocation_class'] = type(params['butcher_tableau']) 

65 params['num_nodes'] = params['butcher_tableau'].num_nodes 

66 

67 # disable residual computation by default 

68 params['skip_residual_computation'] = params.get( 

69 'skip_residual_computation', ('IT_CHECK', 'IT_FINE', 'IT_COARSE', 'IT_UP', 'IT_DOWN') 

70 ) 

71 

72 self.params = _Pars(params) 

73 

74 self.coll = params['butcher_tableau'] 

75 self.coll_bar = params['butcher_tableau_bar'] 

76 

77 # This will be set as soon as the sweeper is instantiated at the level 

78 self.__level = None 

79 

80 self.parallelizable = False 

81 self.QI = self.coll.Qmat 

82 self.Qx = self.coll_bar.Qmat 

83 

84 def get_full_f(self, f): 

85 """ 

86 Test the right hand side funtion is the correct type 

87 

88 Args: 

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

90 

91 Returns: 

92 mesh: Full right hand side as a mesh 

93 """ 

94 

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

96 return f 

97 else: 

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

99 

100 def update_nodes(self): 

101 """ 

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

103 

104 Returns: 

105 None 

106 """ 

107 

108 # get current level and problem description 

109 L = self.level 

110 P = L.prob 

111 

112 # only if the level has been touched before 

113 assert L.status.unlocked 

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

115 

116 # get number of collocation nodes for easier access 

117 M = self.coll.num_nodes 

118 

119 for m in range(0, M): 

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

121 rhs = deepcopy(L.u[0]) 

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

123 

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

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

126 

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

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

129 """ 

130 

131 Implicit part only works for Velocity-Verlet scheme 

132 Boris solver for the implicit part 

133 

134 """ 

135 

136 if self.coll.implicit: 

137 ck = rhs.vel * 0.0 

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

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

140 

141 else: 

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

143 

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

145 L.u[m + 1] = rhs 

146 # update function values 

147 if self.coll.implicit: 

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

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

150 L.f[m + 1] = deepcopy(L.f[0]) 

151 else: 

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

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

154 

155 # indicate presence of new values at this level 

156 

157 L.status.updated = True 

158 

159 return None 

160 

161 def compute_end_point(self): 

162 """ 

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

164 """ 

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

166 

167 

168class RKN(RungeKuttaNystrom): 

169 """ 

170 Runge-Kutta-Nystrom method 

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

172 page: 284 

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

174 """ 

175 

176 def __init__(self, params): 

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

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

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

180 matrix[1, 0] = 0.5 

181 matrix[2, 1] = 0.5 

182 matrix[3, 2] = 1.0 

183 

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

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

186 matrix_bar[1, 0] = 1 / 8 

187 matrix_bar[2, 0] = 1 / 8 

188 matrix_bar[3, 2] = 1 / 2 

189 params['butcher_tableau'] = ButcherTableau(weights, nodes, matrix) 

190 params['butcher_tableau_bar'] = ButcherTableau(weights_bar, nodes, matrix_bar) 

191 

192 super(RKN, self).__init__(params) 

193 

194 

195class Velocity_Verlet(RungeKuttaNystrom): 

196 """ 

197 Velocity-Verlet scheme 

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

199 """ 

200 

201 def __init__(self, params): 

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

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

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

205 matrix[1, 1] = 1 

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

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

208 params['butcher_tableau'] = ButcherTableau(weights, nodes, matrix) 

209 params['butcher_tableau_bar'] = ButcherTableau(weights_bar, nodes, matrix_bar) 

210 params['Velocity_verlet'] = True 

211 

212 super(Velocity_Verlet, self).__init__(params)