Coverage for pySDC/implementations/sweeper_classes/Runge_Kutta_Nystrom.py: 98%
89 statements
« prev ^ index » next coverage.py v7.6.1, created at 2024-09-19 09:13 +0000
« prev ^ index » next coverage.py v7.6.1, created at 2024-09-19 09:13 +0000
1import numpy as np
2import logging
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
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.
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.
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.
27 Please be aware that all fundamental parameters of the Sweeper are ignored. These include
29 - num_nodes
30 - collocation_class
31 - initial_guess
32 - QI
34 All of these variables are either determined by the RK rule, or are not part of an RK scheme.
36 Attribues:
37 butcher_tableau (ButcherTableau): Butcher tableau for the Runge-Kutta scheme that you want
38 """
40 def __init__(self, params):
41 """
42 Initialization routine for the custom sweeper
44 Args:
45 params: parameters for the sweeper
46 """
47 # set up logger
48 self.logger = logging.getLogger('sweeper')
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)
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')
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
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 )
72 self.params = _Pars(params)
74 self.coll = params['butcher_tableau']
75 self.coll_bar = params['butcher_tableau_bar']
77 # This will be set as soon as the sweeper is instantiated at the level
78 self.__level = None
80 self.parallelizable = False
81 self.QI = self.coll.Qmat
82 self.Qx = self.coll_bar.Qmat
84 def get_full_f(self, f):
85 """
86 Test the right hand side funtion is the correct type
88 Args:
89 f (dtype_f): Right hand side at a single node
91 Returns:
92 mesh: Full right hand side as a mesh
93 """
95 if type(f) in [particles, fields, acceleration]:
96 return f
97 else:
98 raise NotImplementedError(f'Type \"{type(f)}\" not implemented')
100 def update_nodes(self):
101 """
102 Update the u- and f-values at the collocation nodes
104 Returns:
105 None
106 """
108 # get current level and problem description
109 L = self.level
110 P = L.prob
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!"
116 # get number of collocation nodes for easier access
117 M = self.coll.num_nodes
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
124 for j in range(1, m + 1):
125 # build RHS from f-terms (containing the E field) and the B field
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 """
131 Implicit part only works for Velocity-Verlet scheme
132 Boris solver for the implicit part
134 """
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])
141 else:
142 rhs.vel += L.dt * self.QI[m + 1, j] * self.get_full_f(f)
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])
155 # indicate presence of new values at this level
157 L.status.updated = True
159 return None
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]
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 """
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
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)
192 super(RKN, self).__init__(params)
195class Velocity_Verlet(RungeKuttaNystrom):
196 """
197 Velocity-Verlet scheme
198 https://de.wikipedia.org/wiki/Verlet-Algorithmus
199 """
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
212 super(Velocity_Verlet, self).__init__(params)