ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4DormandPrinceRK78.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4DormandPrinceRK78.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 // G4DormandPrinceRK78 implementation
27 //
28 // Dormand-Prince 8(7)13M non-FSAL, based on RK scheme from:
29 // P.J. Prince, J.R. Dormand, "High order embedded Runge-Kutta formulae",
30 // Journal of Computational and Applied Mathematics, Volume 7, Issue 1, 1981,
31 // Pages 67-75, ISSN 0377-0427, DOI: 10.1016/0771-050X(81)90010-3
32 //
33 // Created: Somnath Banerjee, Google Summer of Code 2015, 28 June 2015
34 // Supervision: John Apostolakis, CERN
35 // --------------------------------------------------------------------
36 
37 #include "G4DormandPrinceRK78.hh"
38 #include "G4LineSection.hh"
39 
40 // Constructor
41 //
42 G4DormandPrinceRK78::G4DormandPrinceRK78(G4EquationOfMotion* EqRhs,
43  G4int noIntegrationVariables,
44  G4bool primary)
45  : G4MagIntegratorStepper(EqRhs, noIntegrationVariables)
46 {
47  const G4int numberOfVariables = noIntegrationVariables;
48 
49  // New Chunk of memory being created for use by the stepper
50 
51  // aki - for storing intermediate RHS
52  //
53  ak2 = new G4double[numberOfVariables];
54  ak3 = new G4double[numberOfVariables];
55  ak4 = new G4double[numberOfVariables];
56  ak5 = new G4double[numberOfVariables];
57  ak6 = new G4double[numberOfVariables];
58  ak7 = new G4double[numberOfVariables];
59  ak8 = new G4double[numberOfVariables];
60  ak9 = new G4double[numberOfVariables];
61  ak10 = new G4double[numberOfVariables];
62  ak11 = new G4double[numberOfVariables];
63  ak12 = new G4double[numberOfVariables];
64  ak13 = new G4double[numberOfVariables];
65 
66  const G4int numStateVars = std::max(noIntegrationVariables, 8);
67  yTemp = new G4double[numStateVars];
68  yIn = new G4double[numStateVars] ;
69 
70  fLastInitialVector = new G4double[numStateVars] ;
71  fLastFinalVector = new G4double[numStateVars] ;
72 
73  fLastDyDx = new G4double[numStateVars];
74 
75  fMidVector = new G4double[numStateVars];
76  fMidError = new G4double[numStateVars];
77 
78  if( primary )
79  {
80  fAuxStepper = new G4DormandPrinceRK78(EqRhs, numberOfVariables, !primary);
81  }
82 }
83 
84 // Destructor
85 //
86 G4DormandPrinceRK78::~G4DormandPrinceRK78()
87 {
88  // Clear all previously allocated memory for stepper and DistChord
89 
90  delete [] ak2;
91  delete [] ak3;
92  delete [] ak4;
93  delete [] ak5;
94  delete [] ak6;
95  delete [] ak7;
96  delete [] ak8;
97  delete [] ak9;
98  delete [] ak10;
99  delete [] ak11;
100  delete [] ak12;
101  delete [] ak13;
102  delete [] yTemp;
103  delete [] yIn;
104 
105  delete [] fLastInitialVector;
106  delete [] fLastFinalVector;
107  delete [] fLastDyDx;
108  delete [] fMidVector;
109  delete [] fMidError;
110 
111  delete fAuxStepper;
112 }
113 
114 
115 // The following scheme and the set of coefficients have been obtained from
116 // Table2. RK8(7)13M (Rational approximations) from:
117 // P. J. Prince and J. R. Dormand, "High order embedded Runge-Kutta formulae"
118 // Journal of Computational and Applied Math., vol.7, no.1, pp.67-75, 1980.
119 //
120 // Stepper :
121 //
122 // Passing in the value of yInput[],the first time dydx[] and Step length
123 // Giving back yOut and yErr arrays for output and error respectively
124 //
125 void G4DormandPrinceRK78::Stepper(const G4double yInput[],
126  const G4double dydx[],
127  G4double Step,
128  G4double yOut[],
129  G4double yErr[])
130 {
131  G4int i;
132 
133  // The various constants defined on the basis of butcher tableu
134  //
135  const G4double b21 = 1.0/18,
136  b31 = 1.0/48.0 ,
137  b32 = 1.0/16.0 ,
138 
139  b41 = 1.0/32.0 ,
140  b42 = 0.0 ,
141  b43 = 3.0/32.0 ,
142 
143  b51 = 5.0/16.0 ,
144  b52 = 0.0 ,
145  b53 = -75.0/64.0 ,
146  b54 = 75.0/64.0 ,
147 
148  b61 = 3.0/80.0 ,
149  b62 = 0.0 ,
150  b63 = 0.0 ,
151  b64 = 3.0/16.0 ,
152  b65 = 3.0/20.0 ,
153 
154  b71 = 29443841.0/614563906.0 ,
155  b72 = 0.0 ,
156  b73 = 0.0 ,
157  b74 = 77736538.0/692538347.0 ,
158  b75 = -28693883.0/1125000000.0 ,
159  b76 = 23124283.0/1800000000.0 ,
160 
161  b81 = 16016141.0/946692911.0 ,
162  b82 = 0.0 ,
163  b83 = 0.0 ,
164  b84 = 61564180.0/158732637.0 ,
165  b85 = 22789713.0/633445777.0 ,
166  b86 = 545815736.0/2771057229.0 ,
167  b87 = -180193667.0/1043307555.0 ,
168 
169  b91 = 39632708.0/573591083.0 ,
170  b92 = 0.0 ,
171  b93 = 0.0 ,
172  b94 = -433636366.0/683701615.0 ,
173  b95 = -421739975.0/2616292301.0 ,
174  b96 = 100302831.0/723423059.0 ,
175  b97 = 790204164.0/839813087.0 ,
176  b98 = 800635310.0/3783071287.0 ,
177 
178  b101 = 246121993.0/1340847787.0 ,
179  b102 = 0.0 ,
180  b103 = 0.0 ,
181  b104 = -37695042795.0/15268766246.0 ,
182  b105 = -309121744.0/1061227803.0 ,
183  b106 = -12992083.0/490766935.0 ,
184  b107 = 6005943493.0/2108947869.0 ,
185  b108 = 393006217.0/1396673457.0 ,
186  b109 = 123872331.0/1001029789.0 ,
187 
188  b111 = -1028468189.0/846180014.0 ,
189  b112 = 0.0 ,
190  b113 = 0.0 ,
191  b114 = 8478235783.0/508512852.0 ,
192  b115 = 1311729495.0/1432422823.0 ,
193  b116 = -10304129995.0/1701304382.0 ,
194  b117 = -48777925059.0/3047939560.0 ,
195  b118 = 15336726248.0/1032824649.0 ,
196  b119 = -45442868181.0/3398467696.0 ,
197  b1110 = 3065993473.0/597172653.0 ,
198 
199  b121 = 185892177.0/718116043.0 ,
200  b122 = 0.0 ,
201  b123 = 0.0 ,
202  b124 = -3185094517.0/667107341.0 ,
203  b125 = -477755414.0/1098053517.0 ,
204  b126 = -703635378.0/230739211.0 ,
205  b127 = 5731566787.0/1027545527.0 ,
206  b128 = 5232866602.0/850066563.0 ,
207  b129 = -4093664535.0/808688257.0 ,
208  b1210 = 3962137247.0/1805957418.0 ,
209  b1211 = 65686358.0/487910083.0 ,
210 
211  b131 = 403863854.0/491063109.0 ,
212  b132 = 0.0 ,
213  b133 = 0.0 ,
214  b134 = -5068492393.0/434740067.0 ,
215  b135 = -411421997.0/543043805.0 ,
216  b136 = 652783627.0/914296604.0 ,
217  b137 = 11173962825.0/925320556.0 ,
218  b138 = -13158990841.0/6184727034.0 ,
219  b139 = 3936647629.0/1978049680.0 ,
220  b1310 = -160528059.0/685178525.0 ,
221  b1311 = 248638103.0/1413531060.0 ,
222  b1312 = 0.0 ,
223 
224  c1 = 14005451.0/335480064.0 ,
225  // c2 = 0.0 ,
226  // c3 = 0.0 ,
227  // c4 = 0.0 ,
228  // c5 = 0.0 ,
229  c6 = -59238493.0/1068277825.0 ,
230  c7 = 181606767.0/758867731.0 ,
231  c8 = 561292985.0/797845732.0 ,
232  c9 = -1041891430.0/1371343529.0 ,
233  c10 = 760417239.0/1151165299.0 ,
234  c11 = 118820643.0/751138087.0 ,
235  c12 = - 528747749.0/2220607170.0 ,
236  c13 = 1.0/4.0 ,
237 
238  c_1 = 13451932.0/455176623.0 ,
239  // c_2 = 0.0 ,
240  // c_3 = 0.0 ,
241  // c_4 = 0.0 ,
242  // c_5 = 0.0 ,
243  c_6 = -808719846.0/976000145.0 ,
244  c_7 = 1757004468.0/5645159321.0 ,
245  c_8 = 656045339.0/265891186.0 ,
246  c_9 = -3867574721.0/1518517206.0 ,
247  c_10 = 465885868.0/322736535.0 ,
248  c_11 = 53011238.0/667516719.0 ,
249  c_12 = 2.0/45.0 ,
250  c_13 = 0.0 ,
251 
252  dc1 = c_1 - c1 ,
253  // dc2 = c_2 - c2 ,
254  // dc3 = c_3 - c3 ,
255  // dc4 = c_4 - c4 ,
256  // dc5 = c_5 - c5 ,
257  dc6 = c_6 - c6 ,
258  dc7 = c_7 - c7 ,
259  dc8 = c_8 - c8 ,
260  dc9 = c_9 - c9 ,
261  dc10 = c_10 - c10 ,
262  dc11 = c_11 - c11 ,
263  dc12 = c_12 - c12 ,
264  dc13 = c_13 - c13 ;
265  //
266  // end of declaration !
267 
268  const G4int numberOfVariables = GetNumberOfVariables();
269 
270  // The number of variables to be integrated over
271  //
272  yOut[7] = yTemp[7] = yIn[7] = yInput[7];
273 
274  // Saving yInput because yInput and yOut can be aliases for same array
275  //
276  for(i=0; i<numberOfVariables; ++i)
277  {
278  yIn[i]=yInput[i];
279  }
280  // RightHandSide(yIn, dydx) ; // 1st Stage - Not doing, getting passed
281 
282  for(i=0; i<numberOfVariables; ++i)
283  {
284  yTemp[i] = yIn[i] + b21*Step*dydx[i] ;
285  }
286  RightHandSide(yTemp, ak2) ; // 2nd Stage
287 
288  for(i=0; i<numberOfVariables; ++i)
289  {
290  yTemp[i] = yIn[i] + Step*(b31*dydx[i] + b32*ak2[i]) ;
291  }
292  RightHandSide(yTemp, ak3) ; // 3rd Stage
293 
294  for(i=0; i<numberOfVariables; ++i)
295  {
296  yTemp[i] = yIn[i] + Step*(b41*dydx[i] + b42*ak2[i] + b43*ak3[i]) ;
297  }
298  RightHandSide(yTemp, ak4) ; // 4th Stage
299 
300  for(i=0; i<numberOfVariables; ++i)
301  {
302  yTemp[i] = yIn[i] + Step*(b51*dydx[i] + b52*ak2[i] + b53*ak3[i] +
303  b54*ak4[i]) ;
304  }
305  RightHandSide(yTemp, ak5) ; // 5th Stage
306 
307  for(i=0; i<numberOfVariables; ++i)
308  {
309  yTemp[i] = yIn[i] + Step*(b61*dydx[i] + b62*ak2[i] + b63*ak3[i] +
310  b64*ak4[i] + b65*ak5[i]) ;
311  }
312  RightHandSide(yTemp, ak6) ; // 6th Stage
313 
314  for(i=0; i<numberOfVariables; ++i)
315  {
316  yTemp[i] = yIn[i] + Step*(b71*dydx[i] + b72*ak2[i] + b73*ak3[i] +
317  b74*ak4[i] + b75*ak5[i] + b76*ak6[i]);
318  }
319  RightHandSide(yTemp, ak7); // 7th Stage
320 
321  for(i=0; i<numberOfVariables; ++i)
322  {
323  yTemp[i] = yIn[i] + Step*(b81*dydx[i] + b82*ak2[i] + b83*ak3[i] +
324  b84*ak4[i] + b85*ak5[i] + b86*ak6[i] +
325  b87*ak7[i]);
326  }
327  RightHandSide(yTemp, ak8); // 8th Stage
328 
329  for(i=0; i<numberOfVariables; ++i)
330  {
331  yTemp[i] = yIn[i] + Step*(b91*dydx[i] + b92*ak2[i] + b93*ak3[i] +
332  b94*ak4[i] + b95*ak5[i] + b96*ak6[i] +
333  b97*ak7[i] + b98*ak8[i] );
334  }
335  RightHandSide(yTemp, ak9); // 9th Stage
336 
337  for(i=0; i<numberOfVariables; ++i)
338  {
339  yTemp[i] = yIn[i] + Step*(b101*dydx[i] + b102*ak2[i] + b103*ak3[i] +
340  b104*ak4[i] + b105*ak5[i] + b106*ak6[i] +
341  b107*ak7[i] + b108*ak8[i] + b109*ak9[i]);
342  }
343  RightHandSide(yTemp, ak10); // 10th Stage
344 
345  for(i=0; i<numberOfVariables; ++i)
346  {
347  yTemp[i] = yIn[i] + Step*(b111*dydx[i] + b112*ak2[i] + b113*ak3[i] +
348  b114*ak4[i] + b115*ak5[i] + b116*ak6[i] +
349  b117*ak7[i] + b118*ak8[i] + b119*ak9[i] +
350  b1110*ak10[i]);
351  }
352  RightHandSide(yTemp, ak11); // 11th Stage
353 
354  for(i=0; i<numberOfVariables; ++i)
355  {
356  yTemp[i] = yIn[i] + Step*(b121*dydx[i] + b122*ak2[i] + b123*ak3[i] +
357  b124*ak4[i] + b125*ak5[i] + b126*ak6[i] +
358  b127*ak7[i] + b128*ak8[i] + b129*ak9[i] +
359  b1210*ak10[i] + b1211*ak11[i]);
360  }
361  RightHandSide(yTemp, ak12); // 12th Stage
362 
363  for(i=0; i<numberOfVariables; ++i)
364  {
365  yTemp[i] = yIn[i]+Step*(b131*dydx[i] + b132*ak2[i] + b133*ak3[i] +
366  b134*ak4[i] + b135*ak5[i] + b136*ak6[i] +
367  b137*ak7[i] + b138*ak8[i] + b139*ak9[i] +
368  b1310*ak10[i] + b1311*ak11[i] + b1312*ak12[i]);
369  }
370  RightHandSide(yTemp, ak13); // 13th and final Stage
371 
372  for(i=0; i<numberOfVariables; ++i)
373  {
374  // Accumulate increments with proper weights
375 
376  yOut[i] = yIn[i] + Step*(c1*dydx[i] + // c2 * ak2[i] + c3 * ak3[i]
377  // + c4 * ak4[i] + c5 * ak5[i]
378  + c6*ak6[i] +
379  c7*ak7[i] + c8*ak8[i] +c9*ak9[i] + c10*ak10[i]
380  + c11*ak11[i] + c12*ak12[i] + c13*ak13[i]) ;
381 
382  // Estimate error as difference between 7th and 8th order methods
383  //
384  yErr[i] = Step*(dc1*dydx[i] + // dc2*ak2[i] + dc3*ak3[i] + dc4*ak4[i]
385  // + dc5*ak5[i]
386  + dc6*ak6[i] + dc7*ak7[i] + dc8*ak8[i]
387  + dc9*ak9[i] + dc10*ak10[i] + dc11*ak11[i] + dc12*ak12[i]
388  + dc13*ak13[i] ) ;
389 
390  // Store Input and Final values, for possible use in calculating chord
391  //
392  fLastInitialVector[i] = yIn[i] ;
393  fLastFinalVector[i] = yOut[i];
394  fLastDyDx[i] = dydx[i];
395  }
396  fLastStepLength = Step;
397 
398  return ;
399 }
400 
401 // DistChord
402 //
403 G4double G4DormandPrinceRK78::DistChord() const
404 {
405  G4double distLine, distChord;
406  G4ThreeVector initialPoint, finalPoint, midPoint;
407 
408  // Store last initial and final points
409  // (they will be overwritten in self-Stepper call!)
410  //
411  initialPoint = G4ThreeVector( fLastInitialVector[0],
412  fLastInitialVector[1], fLastInitialVector[2]);
413  finalPoint = G4ThreeVector( fLastFinalVector[0],
414  fLastFinalVector[1], fLastFinalVector[2]);
415 
416  // Do half a step using StepNoErr
417 
418  fAuxStepper->Stepper( fLastInitialVector, fLastDyDx, 0.5 * fLastStepLength,
419  fMidVector, fMidError );
420 
421  midPoint = G4ThreeVector( fMidVector[0], fMidVector[1], fMidVector[2]);
422 
423  // Use stored values of Initial and Endpoint + new Midpoint to evaluate
424  // distance of Chord
425  //
426  if (initialPoint != finalPoint)
427  {
428  distLine = G4LineSection::Distline(midPoint, initialPoint, finalPoint);
429  distChord = distLine;
430  }
431  else
432  {
433  distChord = (midPoint-initialPoint).mag();
434  }
435  return distChord;
436 }