ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4DormandPrince745.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4DormandPrince745.cc
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // G4DormandPrince745 implementation
27 //
28 // DormandPrince7 - 5(4) non-FSAL
29 // definition of the stepper() method that evaluates one step in
30 // field propagation.
31 // The coefficients and the algorithm have been adapted from
32 //
33 // J. R. Dormand and P. J. Prince, "A family of embedded Runge-Kutta formulae"
34 // Journal of computational and applied Math., vol.6, no.1, pp.19-26, 1980.
35 //
36 // The Butcher table of the Dormand-Prince-7-4-5 method is as follows :
37 //
38 // 0 |
39 // 1/5 | 1/5
40 // 3/10| 3/40 9/40
41 // 4/5 | 44/45 56/15 32/9
42 // 8/9 | 19372/6561 25360/2187 64448/6561 212/729
43 // 1 | 9017/3168 355/33 46732/5247 49/176 5103/18656
44 // 1 | 35/384 0 500/1113 125/192 2187/6784 11/84
45 // ------------------------------------------------------------------------
46 // 35/384 0 500/1113 125/192 2187/6784 11/84 0
47 // 5179/57600 0 7571/16695 393/640 92097/339200 187/2100 1/40
48 //
49 // Created: Somnath Banerjee, Google Summer of Code 2015, 25 May 2015
50 // Supervision: John Apostolakis, CERN
51 // --------------------------------------------------------------------
52 
53 #include "G4DormandPrince745.hh"
54 #include "G4LineSection.hh"
55 
56 #include <cstring>
57 
58 using namespace field_utils;
59 
61  G4int noIntegrationVariables)
62  : G4MagIntegratorStepper(equation, noIntegrationVariables)
63 {
64 }
65 
67  const G4double dydx[],
68  G4double hstep,
69  G4double yOutput[],
70  G4double yError[],
71  G4double dydxOutput[])
72 {
73  Stepper(yInput, dydx, hstep, yOutput, yError);
74  copy(dydxOutput, ak7);
75 }
76 
77 // Stepper
78 //
79 // Passing in the value of yInput[],the first time dydx[] and Step length
80 // Giving back yOut and yErr arrays for output and error respectively
81 //
83  const G4double dydx[],
84  G4double hstep,
85  G4double yOut[],
86  G4double yErr[])
87 {
88  // The various constants defined on the basis of butcher tableu
89  //
90  const G4double b21 = 0.2,
91  b31 = 3.0 / 40.0, b32 = 9.0 / 40.0,
92  b41 = 44.0 / 45.0, b42 = -56.0 / 15.0, b43 = 32.0/9.0,
93 
94  b51 = 19372.0 / 6561.0, b52 = -25360.0 / 2187.0, b53 = 64448.0 / 6561.0,
95  b54 = -212.0 / 729.0,
96 
97  b61 = 9017.0 / 3168.0 , b62 = -355.0 / 33.0,
98  b63 = 46732.0 / 5247.0, b64 = 49.0 / 176.0,
99  b65 = -5103.0 / 18656.0,
100 
101  b71 = 35.0 / 384.0, b72 = 0.,
102  b73 = 500.0 / 1113.0, b74 = 125.0 / 192.0,
103  b75 = -2187.0 / 6784.0, b76 = 11.0 / 84.0,
104 
105  //Sum of columns, sum(bij) = ei
106  // e1 = 0. ,
107  // e2 = 1.0/5.0 ,
108  // e3 = 3.0/10.0 ,
109  // e4 = 4.0/5.0 ,
110  // e5 = 8.0/9.0 ,
111  // e6 = 1.0 ,
112  // e7 = 1.0 ,
113 
114  // Difference between the higher and the lower order method coeff. :
115  // b7j are the coefficients of higher order
116 
117  dc1 = -(b71 - 5179.0 / 57600.0),
118  dc2 = -(b72 - .0),
119  dc3 = -(b73 - 7571.0 / 16695.0),
120  dc4 = -(b74 - 393.0 / 640.0),
121  dc5 = -(b75 + 92097.0 / 339200.0),
122  dc6 = -(b76 - 187.0 / 2100.0),
123  dc7 = -(- 1.0 / 40.0);
124 
125  const G4int numberOfVariables = GetNumberOfVariables();
126  State yTemp;
127 
128  // The number of variables to be integrated over
129  //
130  yOut[7] = yTemp[7] = yInput[7];
131 
132  // Saving yInput because yInput and yOut can be aliases for same array
133  //
134  for(G4int i = 0; i < numberOfVariables; ++i)
135  {
136  fyIn[i] = yInput[i];
137  }
138  // RightHandSide(yIn, dydx); // Not done! 1st stage
139 
140  for(G4int i = 0; i < numberOfVariables; ++i)
141  {
142  yTemp[i] = fyIn[i] + b21 * hstep * dydx[i];
143  }
144  RightHandSide(yTemp, ak2); // 2nd stage
145 
146  for(G4int i = 0; i < numberOfVariables; ++i)
147  {
148  yTemp[i] = fyIn[i] + hstep * (b31 * dydx[i] + b32 * ak2[i]);
149  }
150  RightHandSide(yTemp, ak3); // 3rd stage
151 
152  for(G4int i = 0; i < numberOfVariables; ++i)
153  {
154  yTemp[i] = fyIn[i] + hstep * (
155  b41 * dydx[i] + b42 * ak2[i] + b43 * ak3[i]);
156  }
157  RightHandSide(yTemp, ak4); // 4th stage
158 
159  for(G4int i = 0; i < numberOfVariables; ++i)
160  {
161  yTemp[i] = fyIn[i] + hstep * (
162  b51 * dydx[i] + b52 * ak2[i] + b53 * ak3[i] + b54 * ak4[i]);
163  }
164  RightHandSide(yTemp, ak5); // 5th stage
165 
166  for(G4int i = 0; i < numberOfVariables; ++i)
167  {
168  yTemp[i] = fyIn[i] + hstep * (
169  b61 * dydx[i] + b62 * ak2[i] +
170  b63 * ak3[i] + b64 * ak4[i] + b65 * ak5[i]);
171  }
172  RightHandSide(yTemp, ak6); // 6th stage
173 
174  for(G4int i = 0; i < numberOfVariables; ++i)
175  {
176  yOut[i] = fyIn[i] + hstep * (
177  b71 * dydx[i] + b72 * ak2[i] + b73 * ak3[i] +
178  b74 * ak4[i] + b75 * ak5[i] + b76 * ak6[i]);
179  }
180  RightHandSide(yOut, ak7); // 7th and Final stage
181 
182  for(G4int i = 0; i < numberOfVariables; ++i)
183  {
184  yErr[i] = hstep * (
185  dc1 * dydx[i] + dc2 * ak2[i] +
186  dc3 * ak3[i] + dc4 * ak4[i] +
187  dc5 * ak5[i] + dc6 * ak6[i] + dc7 * ak7[i]
188  ) + 1.5e-18;
189 
190  // Store Input and Final values, for possible use in calculating chord
191  //
192  fyOut[i] = yOut[i];
193  fdydxIn[i] = dydx[i];
194  }
195 
196  fLastStepLength = hstep;
197 }
198 
200 {
201  // Coefficients were taken from Some Practical Runge-Kutta Formulas
202  // by Lawrence F. Shampine, page 149, c*
203  //
204  const G4double hf1 = 6025192743.0 / 30085553152.0,
205  hf3 = 51252292925.0 / 65400821598.0,
206  hf4 = - 2691868925.0 / 45128329728.0,
207  hf5 = 187940372067.0 / 1594534317056.0,
208  hf6 = - 1776094331.0 / 19743644256.0,
209  hf7 = 11237099.0 / 235043384.0;
210 
211  G4ThreeVector mid;
212 
213  for(G4int i = 0; i < 3; ++i)
214  {
215  mid[i] = fyIn[i] + 0.5 * fLastStepLength * (
216  hf1 * fdydxIn[i] + hf3 * ak3[i] +
217  hf4 * ak4[i] + hf5 * ak5[i] + hf6 * ak6[i] + hf7 * ak7[i]);
218  }
219 
220  const G4ThreeVector begin = makeVector(fyIn, Value3D::Position);
221  const G4ThreeVector end = makeVector(fyOut, Value3D::Position);
222 
223  return G4LineSection::Distline(mid, begin, end);
224 }
225 
226 // The lower (4th) order interpolant given by Dormand and Prince:
227 // J. R. Dormand and P. J. Prince, "Runge-Kutta triples"
228 // Computers & Mathematics with Applications, vol. 12, no. 9,
229 // pp. 1007-1017, 1986.
230 //
233 {
234  const G4int numberOfVariables = GetNumberOfVariables();
235 
236  const G4double tau2 = tau * tau,
237  tau3 = tau * tau2,
238  tau4 = tau2 * tau2;
239 
240  const G4double bf1 = 1.0 / 11282082432.0 * (
241  157015080.0 * tau4 - 13107642775.0 * tau3 + 34969693132.0 * tau2 -
242  32272833064.0 * tau + 11282082432.0);
243 
244  const G4double bf3 = - 100.0 / 32700410799.0 * tau * (
245  15701508.0 * tau3 - 914128567.0 * tau2 + 2074956840.0 * tau -
246  1323431896.0);
247 
248  const G4double bf4 = 25.0 / 5641041216.0 * tau * (
249  94209048.0 * tau3 - 1518414297.0 * tau2 + 2460397220.0 * tau -
250  889289856.0);
251 
252  const G4double bf5 = - 2187.0 / 199316789632.0 * tau * (
253  52338360.0 * tau3 - 451824525.0 * tau2 + 687873124.0 * tau -
254  259006536.0);
255 
256  const G4double bf6 = 11.0 / 2467955532.0 * tau * (
257  106151040.0 * tau3 - 661884105.0 * tau2 +
258  946554244.0 * tau - 361440756.0);
259 
260  const G4double bf7 = 1.0 / 29380423.0 * tau * (1.0 - tau) * (
261  8293050.0 * tau2 - 82437520.0 * tau + 44764047.0);
262 
263  for(G4int i = 0; i < numberOfVariables; ++i)
264  {
265  yOut[i] = fyIn[i] + fLastStepLength * tau * (
266  bf1 * fdydxIn[i] + bf3 * ak3[i] + bf4 * ak4[i] +
267  bf5 * ak5[i] + bf6 * ak6[i] + bf7 * ak7[i]);
268  }
269 }
270 
271 // Following interpolant of order 5 was given by Baker,Dormand,Gilmore, Prince :
272 // T. S. Baker, J. R. Dormand, J. P. Gilmore, and P. J. Prince,
273 // "Continuous approximation with embedded Runge-Kutta methods"
274 // Applied Numerical Mathematics, vol. 22, no. 1, pp. 51-62, 1996.
275 //
276 // Calculating the extra stages for the interpolant
277 //
279 {
280  // Coefficients for the additional stages
281  //
282  const G4double b81 = 6245.0 / 62208.0,
283  b82 = 0.0,
284  b83 = 8875.0 / 103032.0,
285  b84 = -125.0 / 1728.0,
286  b85 = 801.0 / 13568.0,
287  b86 = -13519.0 / 368064.0,
288  b87 = 11105.0 / 368064.0,
289 
290  b91 = 632855.0 / 4478976.0,
291  b92 = 0.0,
292  b93 = 4146875.0 / 6491016.0,
293  b94 = 5490625.0 /14183424.0,
294  b95 = -15975.0 / 108544.0,
295  b96 = 8295925.0 / 220286304.0,
296  b97 = -1779595.0 / 62938944.0,
297  b98 = -805.0 / 4104.0;
298 
299  const G4int numberOfVariables = GetNumberOfVariables();
300  State yTemp;
301 
302  // Evaluate the extra stages
303  //
304  for(G4int i = 0; i < numberOfVariables; ++i)
305  {
306  yTemp[i] = fyIn[i] + fLastStepLength * (
307  b81 * fdydxIn[i] + b82 * ak2[i] + b83 * ak3[i] +
308  b84 * ak4[i] + b85 * ak5[i] + b86 * ak6[i] +
309  b87 * ak7[i]
310  );
311  }
312  RightHandSide(yTemp, ak8); // 8th Stage
313 
314  for(G4int i = 0; i < numberOfVariables; ++i)
315  {
316  yTemp[i] = fyIn[i] + fLastStepLength * (
317  b91 * fdydxIn[i] + b92 * ak2[i] + b93 * ak3[i] +
318  b94 * ak4[i] + b95 * ak5[i] + b96 * ak6[i] +
319  b97 * ak7[i] + b98 * ak8[i]
320  );
321  }
322  RightHandSide(yTemp, ak9); // 9th Stage
323 }
324 
325 // Calculating the interpolated result yOut with the coefficients
326 //
329 {
330  // Define the coefficients for the polynomials
331  //
332  G4double bi[10][5];
333 
334  // COEFFICIENTS OF bi[1]
335  bi[1][0] = 1.0,
336  bi[1][1] = -38039.0 / 7040.0,
337  bi[1][2] = 125923.0 / 10560.0,
338  bi[1][3] = -19683.0 / 1760.0,
339  bi[1][4] = 3303.0 / 880.0,
340  // --------------------------------------------------------
341  //
342  // COEFFICIENTS OF bi[2]
343  bi[2][0] = 0.0,
344  bi[2][1] = 0.0,
345  bi[2][2] = 0.0,
346  bi[2][3] = 0.0,
347  bi[2][4] = 0.0,
348  // --------------------------------------------------------
349  //
350  // COEFFICIENTS OF bi[3]
351  bi[3][0] = 0.0,
352  bi[3][1] = -12500.0 / 4081.0,
353  bi[3][2] = 205000.0 / 12243.0,
354  bi[3][3] = -90000.0 / 4081.0,
355  bi[3][4] = 36000.0 / 4081.0,
356  // --------------------------------------------------------
357  //
358  // COEFFICIENTS OF bi[4]
359  bi[4][0] = 0.0,
360  bi[4][1] = -3125.0 / 704.0,
361  bi[4][2] = 25625.0 / 1056.0,
362  bi[4][3] = -5625.0 / 176.0,
363  bi[4][4] = 1125.0 / 88.0,
364  // --------------------------------------------------------
365  //
366  // COEFFICIENTS OF bi[5]
367  bi[5][0] = 0.0,
368  bi[5][1] = 164025.0 / 74624.0,
369  bi[5][2] = -448335.0 / 37312.0,
370  bi[5][3] = 295245.0 / 18656.0,
371  bi[5][4] = -59049.0 / 9328.0,
372  // --------------------------------------------------------
373  //
374  // COEFFICIENTS OF bi[6]
375  bi[6][0] = 0.0,
376  bi[6][1] = -25.0 / 28.0,
377  bi[6][2] = 205.0 / 42.0,
378  bi[6][3] = -45.0 / 7.0,
379  bi[6][4] = 18.0 / 7.0,
380  // --------------------------------------------------------
381  //
382  // COEFFICIENTS OF bi[7]
383  bi[7][0] = 0.0,
384  bi[7][1] = -2.0 / 11.0,
385  bi[7][2] = 73.0 / 55.0,
386  bi[7][3] = -171.0 / 55.0,
387  bi[7][4] = 108.0 / 55.0,
388  // --------------------------------------------------------
389  //
390  // COEFFICIENTS OF bi[8]
391  bi[8][0] = 0.0,
392  bi[8][1] = 189.0 / 22.0,
393  bi[8][2] = -1593.0 / 55.0,
394  bi[8][3] = 3537.0 / 110.0,
395  bi[8][4] = -648.0 / 55.0,
396  // --------------------------------------------------------
397  //
398  // COEFFICIENTS OF bi[9]
399  bi[9][0] = 0.0,
400  bi[9][1] = 351.0 / 110.0,
401  bi[9][2] = -999.0 / 55.0,
402  bi[9][3] = 2943.0 / 110.0,
403  bi[9][4] = -648.0 / 55.0;
404  // --------------------------------------------------------
405 
406  // Calculating the polynomials
407 
408  G4double b[10];
409  std::memset(b, 0.0, sizeof(b));
410 
411  G4double tauPower = 1.0;
412  for(G4int j = 0; j <= 4; ++j)
413  {
414  for(G4int iStage = 1; iStage <= 9; ++iStage)
415  {
416  b[iStage] += bi[iStage][j] * tauPower;
417  }
418  tauPower *= tau;
419  }
420 
421  const G4int numberOfVariables = GetNumberOfVariables();
422  const G4double stepLen = fLastStepLength * tau;
423  for(G4int i = 0; i < numberOfVariables; ++i)
424  {
425  yOut[i] = fyIn[i] + stepLen * (
426  b[1] * fdydxIn[i] + b[2] * ak2[i] + b[3] * ak3[i] +
427  b[4] * ak4[i] + b[5] * ak5[i] + b[6] * ak6[i] +
428  b[7] * ak7[i] + b[8] * ak8[i] + b[9] * ak9[i]
429  );
430  }
431 }