CMS 3D CMS Logo

CMSTDormandPrince45.h
Go to the documentation of this file.
1 //
2 // copied from G4TDormandPrince45
3 //
4 // Class desription:
5 //
6 // An implementation of the 5th order embedded RK method from the paper:
7 // J. R. Dormand and P. J. Prince, "A family of embedded Runge-Kutta formulae"
8 // Journal of computational and applied Math., vol.6, no.1, pp.19-26, 1980.
9 //
10 // DormandPrince7 - 5(4) embedded RK method
11 //
12 //
13 // Created: Somnath Banerjee, Google Summer of Code 2015, 25 May 2015
14 // Supervision: John Apostolakis, CERN
15 //
16 // Modified by Vincenzo Innocente on 7/2022 to account for constant momentum magnitude
17 //
18 // --------------------------------------------------------------------
19 #ifndef CMSTDORMAND_PRINCE_45_HH
20 #define CMSTDORMAND_PRINCE_45_HH
21 
22 #include <cassert>
23 #include "G4MagIntegratorStepper.hh"
24 #include "G4FieldUtils.hh"
25 
26 template <class T_Equation, unsigned int N = 6>
27 class CMSTDormandPrince45 : public G4MagIntegratorStepper {
28 public:
29  CMSTDormandPrince45(T_Equation* equation);
30  CMSTDormandPrince45(T_Equation* equation, G4int numVar); // must have numVar == N
31 
32  inline void StepWithError(
33  const G4double yInput[], const G4double dydx[], G4double hstep, G4double yOutput[], G4double yError[]);
34 
35  void Stepper(
36  const G4double yInput[], const G4double dydx[], G4double hstep, G4double yOutput[], G4double yError[]) final;
37 
38  inline void StepWithFinalDerivate(const G4double yInput[],
39  const G4double dydx[],
40  G4double hstep,
41  G4double yOutput[],
42  G4double yError[],
43  G4double dydxOutput[]);
44 
45  inline void SetupInterpolation() {}
46 
47  void Interpolate(G4double tau, G4double yOut[]) const { Interpolate4thOrder(yOut, tau); }
48  // For calculating the output at the tau fraction of Step
49 
50  G4double DistChord() const final;
51 
52  G4int IntegratorOrder() const override { return 4; }
53 
54  const field_utils::ShortState<N>& GetYOut() const { return fyOut; }
55 
56  void Interpolate4thOrder(G4double yOut[], G4double tau) const;
57 
59  void Interpolate5thOrder(G4double yOut[], G4double tau) const;
60 
61  // __attribute__((always_inline))
62  void RightHandSideInl(const G4double y[], G4double inv_momentum_magnitude, G4double dydx[]) {
63  fEquation_Rhs->T_Equation::TRightHandSide(y, inv_momentum_magnitude, dydx);
64  }
65 
66  inline void Stepper(const G4double yInput[],
67  const G4double dydx[],
68  G4double hstep,
69  G4double yOutput[],
70  G4double yError[],
71  G4double dydxOutput[]) {
72  StepWithFinalDerivate(yInput, dydx, hstep, yOutput, yError, dydxOutput);
73  }
74 
75  T_Equation* GetSpecificEquation() { return fEquation_Rhs; }
76 
77  static constexpr int N8 = N > 8 ? N : 8; // y[
78 
79 private:
80  field_utils::ShortState<N> ak2, ak3, ak4, ak5, ak6, ak7, ak8, ak9;
81  field_utils::ShortState<N8> fyIn;
82  field_utils::ShortState<N> fyOut, fdydxIn;
83 
84  // - Simpler :
85  // field_utils::State ak2, ak3, ak4, ak5, ak6, ak7, ak8, ak9;
86  // field_utils::State fyIn, fyOut, fdydxIn;
87 
88  G4double fLastStepLength = -1.0;
89  T_Equation* fEquation_Rhs;
90 };
91 
92 // G4TDormandPrince745 implementation -- borrowed from G4DormandPrince745
93 //
94 // DormandPrince7 - 5(4) non-FSAL
95 // definition of the stepper() method that evaluates one step in
96 // field propagation.
97 // The coefficients and the algorithm have been adapted from
98 //
99 // J. R. Dormand and P. J. Prince, "A family of embedded Runge-Kutta formulae"
100 // Journal of computational and applied Math., vol.6, no.1, pp.19-26, 1980.
101 //
102 // The Butcher table of the Dormand-Prince-7-4-5 method is as follows :
103 //
104 // 0 |
105 // 1/5 | 1/5
106 // 3/10| 3/40 9/40
107 // 4/5 | 44/45 56/15 32/9
108 // 8/9 | 19372/6561 25360/2187 64448/6561 212/729
109 // 1 | 9017/3168 355/33 46732/5247 49/176 5103/18656
110 // 1 | 35/384 0 500/1113 125/192 2187/6784 11/84
111 // ------------------------------------------------------------------------
112 // 35/384 0 500/1113 125/192 2187/6784 11/84 0
113 // 5179/57600 0 7571/16695 393/640 92097/339200 187/2100 1/40
114 //
115 // Created: Somnath Banerjee, Google Summer of Code 2015, 25 May 2015
116 // Supervision: John Apostolakis, CERN
117 // --------------------------------------------------------------------
118 
119 #include "G4LineSection.hh"
120 
121 #include <cstring>
122 
123 // using namespace field_utils;
124 
126 // Constructor
127 //
128 template <class T_Equation, unsigned int N>
130  : G4MagIntegratorStepper(dynamic_cast<G4EquationOfMotion*>(equation), N), fEquation_Rhs(equation) {
131  // assert( dynamic_cast<G4EquationOfMotion*>(equation) != nullptr );
132  if (dynamic_cast<G4EquationOfMotion*>(equation) == nullptr) {
133  G4Exception("G4TDormandPrince745CMS: constructor",
134  "GeomField0001",
135  FatalException,
136  "T_Equation is not an G4EquationOfMotion.");
137  }
138 
139  /***
140  assert( equation->GetNumberOfVariables == N );
141  if( equation->GetNumberOfVariables != N ){
142  G4ExceptionDescription msg;
143  msg << "Equation has an incompatible number of variables." ;
144  msg << " template N = " << N << " equation-Nvar= "
145  << equation->GetNumberOfVariables;
146  G4Exception("G4TCashKarpRKF45: constructor", "GeomField0001",
147  FatalException, msg );
148  } ****/
149 }
150 
151 template <class T_Equation, unsigned int N>
153  : CMSTDormandPrince45<T_Equation, N>(equation) {
154  if (numVar != G4int(N)) {
155  G4ExceptionDescription msg;
156  msg << "Equation has an incompatible number of variables.";
157  msg << " template N = " << N << " argument numVar = " << numVar;
158  // << " equation-Nvar= " << equation->GetNumberOfVariables(); // --> Expected later
159  G4Exception("G4TCashKarpRKF45CMS: constructor", "GeomField0001", FatalErrorInArgument, msg);
160  }
161  assert(numVar == N);
162 }
163 
164 template <class T_Equation, unsigned int N>
166  const G4double dydx[],
167  G4double hstep,
168  G4double yOutput[],
169  G4double yError[],
170  G4double dydxOutput[]) {
171  StepWithError(yInput, dydx, hstep, yOutput, yError);
172  field_utils::copy(dydxOutput, ak7, N);
173 }
174 
175 // Stepper
176 //
177 // Passing in the value of yInput[],the first time dydx[] and Step length
178 // Giving back yOut and yErr arrays for output and error respectively
179 //
180 
181 template <class T_Equation, unsigned int N>
183  const G4double yInput[], const G4double dydx[], G4double hstep, G4double yOut[], G4double yErr[]) {
184  // The parameters of the Butcher tableu
185  //
186  constexpr G4double b21 = 0.2, b31 = 3.0 / 40.0, b32 = 9.0 / 40.0, b41 = 44.0 / 45.0, b42 = -56.0 / 15.0,
187  b43 = 32.0 / 9.0,
188 
189  b51 = 19372.0 / 6561.0, b52 = -25360.0 / 2187.0, b53 = 64448.0 / 6561.0, b54 = -212.0 / 729.0,
190 
191  b61 = 9017.0 / 3168.0, b62 = -355.0 / 33.0, b63 = 46732.0 / 5247.0, b64 = 49.0 / 176.0,
192  b65 = -5103.0 / 18656.0,
193 
194  b71 = 35.0 / 384.0, b72 = 0., b73 = 500.0 / 1113.0, b74 = 125.0 / 192.0, b75 = -2187.0 / 6784.0,
195  b76 = 11.0 / 84.0,
196 
197  //Sum of columns, sum(bij) = ei
198  // e1 = 0. ,
199  // e2 = 1.0/5.0 ,
200  // e3 = 3.0/10.0 ,
201  // e4 = 4.0/5.0 ,
202  // e5 = 8.0/9.0 ,
203  // e6 = 1.0 ,
204  // e7 = 1.0 ,
205 
206  // Difference between the higher and the lower order method coeff. :
207  // b7j are the coefficients of higher order
208 
209  dc1 = -(b71 - 5179.0 / 57600.0), dc2 = -(b72 - .0), dc3 = -(b73 - 7571.0 / 16695.0), dc4 = -(b74 - 393.0 / 640.0),
210  dc5 = -(b75 + 92097.0 / 339200.0), dc6 = -(b76 - 187.0 / 2100.0), dc7 = -(-1.0 / 40.0);
211 
212  // const G4int numberOfVariables = GetNumberOfVariables();
213  // The number of variables to be integrated over
214  field_utils::ShortState<N8> yTemp;
215 
216  yOut[7] = yTemp[7] = fyIn[7] = yInput[7]; // Pass along the time - used in RightHandSide
217 
218  // Saving yInput because yInput and yOut can be aliases for same array
219  //
220  for (unsigned int i = 0; i < N; ++i) {
221  fyIn[i] = yInput[i];
222  yTemp[i] = yInput[i] + b21 * hstep * dydx[i];
223  }
224  G4double momentum_mag_square = yTemp[3] * yTemp[3] + yTemp[4] * yTemp[4] + yTemp[5] * yTemp[5];
225  G4double inv_momentum_magnitude = 1.0 / std::sqrt(momentum_mag_square);
226  RightHandSideInl(yTemp, inv_momentum_magnitude, ak2); // 2nd stage
227 
228  for (unsigned int i = 0; i < N; ++i) {
229  yTemp[i] = fyIn[i] + hstep * (b31 * dydx[i] + b32 * ak2[i]);
230  }
231  RightHandSideInl(yTemp, inv_momentum_magnitude, ak3); // 3rd stage
232 
233  for (unsigned int i = 0; i < N; ++i) {
234  yTemp[i] = fyIn[i] + hstep * (b41 * dydx[i] + b42 * ak2[i] + b43 * ak3[i]);
235  }
236  RightHandSideInl(yTemp, inv_momentum_magnitude, ak4); // 4th stage
237 
238  for (unsigned int i = 0; i < N; ++i) {
239  yTemp[i] = fyIn[i] + hstep * (b51 * dydx[i] + b52 * ak2[i] + b53 * ak3[i] + b54 * ak4[i]);
240  }
241  RightHandSideInl(yTemp, inv_momentum_magnitude, ak5); // 5th stage
242 
243  for (unsigned int i = 0; i < N; ++i) {
244  yTemp[i] = fyIn[i] + hstep * (b61 * dydx[i] + b62 * ak2[i] + b63 * ak3[i] + b64 * ak4[i] + b65 * ak5[i]);
245  }
246  RightHandSideInl(yTemp, inv_momentum_magnitude, ak6); // 6th stage
247 
248  for (unsigned int i = 0; i < N; ++i) {
249  yOut[i] =
250  fyIn[i] + hstep * (b71 * dydx[i] + b72 * ak2[i] + b73 * ak3[i] + b74 * ak4[i] + b75 * ak5[i] + b76 * ak6[i]);
251  }
252  RightHandSideInl(yOut, inv_momentum_magnitude, ak7); // 7th and Final stage
253 
254  for (unsigned int i = 0; i < N; ++i) {
255  yErr[i] = hstep * (dc1 * dydx[i] + dc2 * ak2[i] + dc3 * ak3[i] + dc4 * ak4[i] + dc5 * ak5[i] + dc6 * ak6[i] +
256  dc7 * ak7[i]) +
257  1.5e-18;
258 
259  // Store Input and Final values, for possible use in calculating chord
260  //
261  fyOut[i] = yOut[i];
262  fdydxIn[i] = dydx[i];
263  }
264 
265  fLastStepLength = hstep;
266 }
267 
268 template <class T_Equation, unsigned int N>
270  const G4double yInput[], const G4double dydx[], G4double Step, G4double yOutput[], G4double yError[]) {
271  assert(yOutput != yInput);
272  assert(yError != yInput);
273 
274  StepWithError(yInput, dydx, Step, yOutput, yError);
275 }
276 
277 template <class T_Equation, unsigned int N>
279  // Coefficients were taken from Some Practical Runge-Kutta Formulas
280  // by Lawrence F. Shampine, page 149, c*
281  //
282  const G4double hf1 = 6025192743.0 / 30085553152.0, hf3 = 51252292925.0 / 65400821598.0,
283  hf4 = -2691868925.0 / 45128329728.0, hf5 = 187940372067.0 / 1594534317056.0,
284  hf6 = -1776094331.0 / 19743644256.0, hf7 = 11237099.0 / 235043384.0;
285 
286  G4ThreeVector mid;
287 
288  for (unsigned int i = 0; i < 3; ++i) {
289  mid[i] =
290  fyIn[i] + 0.5 * fLastStepLength *
291  (hf1 * fdydxIn[i] + hf3 * ak3[i] + hf4 * ak4[i] + hf5 * ak5[i] + hf6 * ak6[i] + hf7 * ak7[i]);
292  }
293 
294  const G4ThreeVector begin = makeVector(fyIn, field_utils::Value3D::Position);
295  const G4ThreeVector end = makeVector(fyOut, field_utils::Value3D::Position);
296 
297  return G4LineSection::Distline(mid, begin, end);
298 }
299 
300 // The lower (4th) order interpolant given by Dormand and Prince:
301 // J. R. Dormand and P. J. Prince, "Runge-Kutta triples"
302 // Computers & Mathematics with Applications, vol. 12, no. 9,
303 // pp. 1007-1017, 1986.
304 //
305 template <class T_Equation, unsigned int N>
306 void CMSTDormandPrince45<T_Equation, N>::Interpolate4thOrder(G4double yOut[], G4double tau) const {
307  // const G4int numberOfVariables = GetNumberOfVariables();
308 
309  const G4double tau2 = tau * tau, tau3 = tau * tau2, tau4 = tau2 * tau2;
310 
311  const G4double bf1 =
312  1.0 / 11282082432.0 *
313  (157015080.0 * tau4 - 13107642775.0 * tau3 + 34969693132.0 * tau2 - 32272833064.0 * tau + 11282082432.0);
314 
315  const G4double bf3 =
316  -100.0 / 32700410799.0 * tau * (15701508.0 * tau3 - 914128567.0 * tau2 + 2074956840.0 * tau - 1323431896.0);
317 
318  const G4double bf4 =
319  25.0 / 5641041216.0 * tau * (94209048.0 * tau3 - 1518414297.0 * tau2 + 2460397220.0 * tau - 889289856.0);
320 
321  const G4double bf5 =
322  -2187.0 / 199316789632.0 * tau * (52338360.0 * tau3 - 451824525.0 * tau2 + 687873124.0 * tau - 259006536.0);
323 
324  const G4double bf6 =
325  11.0 / 2467955532.0 * tau * (106151040.0 * tau3 - 661884105.0 * tau2 + 946554244.0 * tau - 361440756.0);
326 
327  const G4double bf7 = 1.0 / 29380423.0 * tau * (1.0 - tau) * (8293050.0 * tau2 - 82437520.0 * tau + 44764047.0);
328 
329  for (unsigned int i = 0; i < N; ++i) {
330  yOut[i] =
331  fyIn[i] + fLastStepLength * tau *
332  (bf1 * fdydxIn[i] + bf3 * ak3[i] + bf4 * ak4[i] + bf5 * ak5[i] + bf6 * ak6[i] + bf7 * ak7[i]);
333  }
334 }
335 
336 // Following interpolant of order 5 was given by Baker,Dormand,Gilmore, Prince :
337 // T. S. Baker, J. R. Dormand, J. P. Gilmore, and P. J. Prince,
338 // "Continuous approximation with embedded Runge-Kutta methods"
339 // Applied Numerical Mathematics, vol. 22, no. 1, pp. 51-62, 1996.
340 //
341 // Calculating the extra stages for the interpolant
342 //
343 template <class T_Equation, unsigned int N>
345  // Coefficients for the additional stages
346  //
347  const G4double b81 = 6245.0 / 62208.0, b82 = 0.0, b83 = 8875.0 / 103032.0, b84 = -125.0 / 1728.0,
348  b85 = 801.0 / 13568.0, b86 = -13519.0 / 368064.0, b87 = 11105.0 / 368064.0,
349 
350  b91 = 632855.0 / 4478976.0, b92 = 0.0, b93 = 4146875.0 / 6491016.0, b94 = 5490625.0 / 14183424.0,
351  b95 = -15975.0 / 108544.0, b96 = 8295925.0 / 220286304.0, b97 = -1779595.0 / 62938944.0,
352  b98 = -805.0 / 4104.0;
353 
354  // const G4int numberOfVariables = GetNumberOfVariables();
355  field_utils::ShortState<N> yTemp;
356 
357  // Evaluate the extra stages
358  //
359  for (unsigned int i = 0; i < N; ++i) {
360  yTemp[i] = fyIn[i] + fLastStepLength * (b81 * fdydxIn[i] + b82 * ak2[i] + b83 * ak3[i] + b84 * ak4[i] +
361  b85 * ak5[i] + b86 * ak6[i] + b87 * ak7[i]);
362  }
363  G4double momentum_mag_square = yTemp[3] * yTemp[3] + yTemp[4] * yTemp[4] + yTemp[5] * yTemp[5];
364  G4double inv_momentum_magnitude = 1.0 / std::sqrt(momentum_mag_square);
365  RightHandSideInl(yTemp, inv_momentum_magnitude, ak8); // 8th Stage
366 
367  for (unsigned int i = 0; i < N; ++i) {
368  yTemp[i] = fyIn[i] + fLastStepLength * (b91 * fdydxIn[i] + b92 * ak2[i] + b93 * ak3[i] + b94 * ak4[i] +
369  b95 * ak5[i] + b96 * ak6[i] + b97 * ak7[i] + b98 * ak8[i]);
370  }
371  RightHandSideInl(yTemp, inv_momentum_magnitude, ak9); // 9th Stage
372 }
373 
374 // Calculating the interpolated result yOut with the coefficients
375 //
376 template <class T_Equation, unsigned int N>
377 void CMSTDormandPrince45<T_Equation, N>::Interpolate5thOrder(G4double yOut[], G4double tau) const {
378  // Define the coefficients for the polynomials
379  //
380  G4double bi[10][5];
381 
382  // COEFFICIENTS OF bi[1]
383  bi[1][0] = 1.0, bi[1][1] = -38039.0 / 7040.0, bi[1][2] = 125923.0 / 10560.0, bi[1][3] = -19683.0 / 1760.0,
384  bi[1][4] = 3303.0 / 880.0,
385  // --------------------------------------------------------
386  //
387  // COEFFICIENTS OF bi[2]
388  bi[2][0] = 0.0, bi[2][1] = 0.0, bi[2][2] = 0.0, bi[2][3] = 0.0, bi[2][4] = 0.0,
389  // --------------------------------------------------------
390  //
391  // COEFFICIENTS OF bi[3]
392  bi[3][0] = 0.0, bi[3][1] = -12500.0 / 4081.0, bi[3][2] = 205000.0 / 12243.0, bi[3][3] = -90000.0 / 4081.0,
393  bi[3][4] = 36000.0 / 4081.0,
394  // --------------------------------------------------------
395  //
396  // COEFFICIENTS OF bi[4]
397  bi[4][0] = 0.0, bi[4][1] = -3125.0 / 704.0, bi[4][2] = 25625.0 / 1056.0, bi[4][3] = -5625.0 / 176.0,
398  bi[4][4] = 1125.0 / 88.0,
399  // --------------------------------------------------------
400  //
401  // COEFFICIENTS OF bi[5]
402  bi[5][0] = 0.0, bi[5][1] = 164025.0 / 74624.0, bi[5][2] = -448335.0 / 37312.0, bi[5][3] = 295245.0 / 18656.0,
403  bi[5][4] = -59049.0 / 9328.0,
404  // --------------------------------------------------------
405  //
406  // COEFFICIENTS OF bi[6]
407  bi[6][0] = 0.0, bi[6][1] = -25.0 / 28.0, bi[6][2] = 205.0 / 42.0, bi[6][3] = -45.0 / 7.0, bi[6][4] = 18.0 / 7.0,
408  // --------------------------------------------------------
409  //
410  // COEFFICIENTS OF bi[7]
411  bi[7][0] = 0.0, bi[7][1] = -2.0 / 11.0, bi[7][2] = 73.0 / 55.0, bi[7][3] = -171.0 / 55.0, bi[7][4] = 108.0 / 55.0,
412  // --------------------------------------------------------
413  //
414  // COEFFICIENTS OF bi[8]
415  bi[8][0] = 0.0, bi[8][1] = 189.0 / 22.0, bi[8][2] = -1593.0 / 55.0, bi[8][3] = 3537.0 / 110.0,
416  bi[8][4] = -648.0 / 55.0,
417  // --------------------------------------------------------
418  //
419  // COEFFICIENTS OF bi[9]
420  bi[9][0] = 0.0, bi[9][1] = 351.0 / 110.0, bi[9][2] = -999.0 / 55.0, bi[9][3] = 2943.0 / 110.0,
421  bi[9][4] = -648.0 / 55.0;
422  // --------------------------------------------------------
423 
424  // Calculating the polynomials
425 
426  G4double b[10];
427  std::memset(b, 0.0, sizeof(b));
428 
429  G4double tauPower = 1.0;
430  for (G4int j = 0; j <= 4; ++j) {
431  for (G4int iStage = 1; iStage <= 9; ++iStage) {
432  b[iStage] += bi[iStage][j] * tauPower;
433  }
434  tauPower *= tau;
435  }
436 
437  // const G4int numberOfVariables = GetNumberOfVariables();
438  const G4double stepLen = fLastStepLength * tau;
439  for (G4int i = 0; i < N; ++i) {
440  yOut[i] = fyIn[i] + stepLen * (b[1] * fdydxIn[i] + b[2] * ak2[i] + b[3] * ak3[i] + b[4] * ak4[i] + b[5] * ak5[i] +
441  b[6] * ak6[i] + b[7] * ak7[i] + b[8] * ak8[i] + b[9] * ak9[i]);
442  }
443 }
444 
445 #endif
void Stepper(const G4double yInput[], const G4double dydx[], G4double hstep, G4double yOutput[], G4double yError[], G4double dydxOutput[])
static const signed char b64[0x100]
Definition: base64.cc:208
field_utils::ShortState< N > ak7
const field_utils::ShortState< N > & GetYOut() const
assert(be >=bs)
field_utils::ShortState< N > ak3
field_utils::ShortState< N > fyOut
void RightHandSideInl(const G4double y[], G4double inv_momentum_magnitude, G4double dydx[])
T_Equation * GetSpecificEquation()
field_utils::ShortState< N > ak9
G4double DistChord() const final
void Interpolate(G4double tau, G4double yOut[]) const
field_utils::ShortState< N > ak5
static constexpr int N8
void StepWithFinalDerivate(const G4double yInput[], const G4double dydx[], G4double hstep, G4double yOutput[], G4double yError[], G4double dydxOutput[])
T sqrt(T t)
Definition: SSEVec.h:19
CMSTDormandPrince45(T_Equation *equation)
field_utils::ShortState< N > ak8
G4int IntegratorOrder() const override
void StepWithError(const G4double yInput[], const G4double dydx[], G4double hstep, G4double yOutput[], G4double yError[])
void Interpolate5thOrder(G4double yOut[], G4double tau) const
field_utils::ShortState< N > ak4
#define N
Definition: blowfish.cc:9
double b
Definition: hdecay.h:120
tuple msg
Definition: mps_check.py:286
void Interpolate4thOrder(G4double yOut[], G4double tau) const
field_utils::ShortState< N > fdydxIn
field_utils::ShortState< N8 > fyIn
field_utils::ShortState< N > ak2
field_utils::ShortState< N > ak6
void Stepper(const G4double yInput[], const G4double dydx[], G4double hstep, G4double yOutput[], G4double yError[]) final