ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4FSALDormandPrince745.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4FSALDormandPrince745.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 // G4FSALDormandPrince745 implementation
27 //
28 // The Butcher table of the FDormand-Prince-7-4-5 method is as follows:
29 //
30 // 0 |
31 // 1/5 | 1/5
32 // 3/10| 3/40 9/40
33 // 4/5 | 44/45 56/15 32/9
34 // 8/9 | 19372/6561 25360/2187 64448/6561 212/729
35 // 1 | 9017/3168 355/33 46732/5247 49/176 5103/18656
36 // 1 | 35/384 0 500/1113 125/192 2187/6784 11/84
37 // ---------------------------------------------------------------------------
38 // 35/384 0 500/1113 125/192 2187/6784 11/84 0
39 // 5179/57600 0 7571/16695 393/640 92097/339200 187/2100 1/40
40 //
41 // Created: Somnath Banerjee, Google Summer of Code 2015, 25 May 2015
42 // Supervision: John Apostolakis, CERN
43 // --------------------------------------------------------------------
44 
45 #include "G4FSALDormandPrince745.hh"
46 #include "G4LineSection.hh"
47 #include <cmath>
48 
49 // Constructor
50 //
51 G4FSALDormandPrince745::G4FSALDormandPrince745(G4EquationOfMotion* EqRhs,
52  G4int noIntegrationVariables,
53  G4bool primary)
54  : G4VFSALIntegrationStepper(EqRhs, noIntegrationVariables)
55 {
56  const G4int numberOfVariables = noIntegrationVariables;
57 
58  // New Chunk of memory being created for use by the stepper
59 
60  // aki - for storing intermediate RHS
61  //
62  ak2 = new G4double[numberOfVariables];
63  ak3 = new G4double[numberOfVariables];
64  ak4 = new G4double[numberOfVariables];
65  ak5 = new G4double[numberOfVariables];
66  ak6 = new G4double[numberOfVariables];
67  ak7 = new G4double[numberOfVariables];
68 
69  // Also always allocate arrays for interpolation stages
70  //
71  ak8 = new G4double[numberOfVariables];
72  ak9 = new G4double[numberOfVariables];
73 
74  yTemp = new G4double[numberOfVariables] ;
75  yIn = new G4double[numberOfVariables] ;
76 
77  pseudoDydx_for_DistChord = new G4double[numberOfVariables];
78 
79  fInitialDyDx = new G4double[numberOfVariables];
80  fLastInitialVector = new G4double[numberOfVariables] ;
81  fLastFinalVector = new G4double[numberOfVariables] ;
82  fLastDyDx = new G4double[numberOfVariables];
83 
84  fMidVector = new G4double[numberOfVariables];
85  fMidError = new G4double[numberOfVariables];
86 
87  if( primary )
88  {
89  fAuxStepper = new G4FSALDormandPrince745(EqRhs,numberOfVariables,!primary);
90  }
91 }
92 
93 // Destructor
94 //
95 G4FSALDormandPrince745::~G4FSALDormandPrince745()
96 {
97  // Clear all previously allocated memory for stepper and DistChord
98 
99  delete [] ak2; ak2 = nullptr;
100  delete [] ak3; ak3 = nullptr;
101  delete [] ak4; ak4 = nullptr;
102  delete [] ak5; ak5 = nullptr;
103  delete [] ak6; ak6 = nullptr;
104  delete [] ak7; ak7 = nullptr;
105  delete [] ak8; ak8 = nullptr;
106  delete [] ak9; ak9 = nullptr;
107 
108  delete [] yTemp; yTemp = nullptr;
109  delete [] yIn; yIn = nullptr;
110 
111  delete [] pseudoDydx_for_DistChord; pseudoDydx_for_DistChord = nullptr;
112  delete [] fInitialDyDx; fInitialDyDx = nullptr;
113 
114  delete [] fLastInitialVector; fLastInitialVector = nullptr;
115  delete [] fLastFinalVector; fLastFinalVector = nullptr;
116  delete [] fLastDyDx; fLastDyDx = nullptr;
117  delete [] fMidVector; fMidVector = nullptr;
118  delete [] fMidError; fMidError = nullptr;
119 
120  delete fAuxStepper; fAuxStepper = nullptr;
121 }
122 
123 // Stepper
124 //
125 // Passing in the value of yInput[],the first time dydx[] and Step length
126 // Giving back yOut and yErr arrays for output and error respectively
127 //
128 void G4FSALDormandPrince745::Stepper(const G4double yInput[],
129  const G4double dydx[],
130  G4double Step,
131  G4double yOut[],
132  G4double yErr[],
133  G4double nextDydx[] )
134 {
135  G4int i;
136 
137  // The various constants defined on the basis of butcher tableu
138 
139  const G4double b21 = 0.2 ,
140  b31 = 3.0/40.0, b32 = 9.0/40.0 ,
141 
142  b41 = 44.0/45.0, b42 = -56.0/15.0, b43 = 32.0/9.0,
143 
144  b51 = 19372.0/6561.0, b52 = -25360.0/2187.0,
145  b53 = 64448.0/6561.0, b54 = -212.0/729.0 ,
146 
147  b61 = 9017.0/3168.0 , b62 = -355.0/33.0,
148  b63 = 46732.0/5247.0 , b64 = 49.0/176.0 ,
149  b65 = -5103.0/18656.0 ,
150 
151  b71 = 35.0/384.0, b72 = 0.,
152  b73 = 500.0/1113.0, b74 = 125.0/192.0,
153  b75 = -2187.0/6784.0, b76 = 11.0/84.0,
154 
155  // c1 = 35.0/384.0, c2 = .0,
156  // c3 = 500.0/1113.0, c4 = 125.0/192.0,
157  // c5 = -2187.0/6784.0, c6 = 11.0/84.0,
158  // c7 = 0,
159 
160  dc1 = b71 - 5179.0/57600.0,
161  dc2 = b72 - .0,
162  dc3 = b73 - 7571.0/16695.0,
163  dc4 = b74 - 393.0/640.0,
164  dc5 = b75 + 92097.0/339200.0,
165  dc6 = b76 - 187.0/2100.0,
166  dc7 = - 1.0/40.0 ; //end of declaration
167 
168  const G4int numberOfVariables = GetNumberOfVariables();
169  // The number of variables to be integrated over
170 
171  // Saving yInput because yInput and yOut can be aliases for same array
172  //
173  for(i=0; i<numberOfVariables; ++i)
174  {
175  yIn[i] = yInput[i];
176  fInitialDyDx[i] = dydx[i];
177  }
178  // Ensure that time is initialised - in case it is not integrated
179  //
180  yOut[7] = yTemp[7] = yInput[7];
181  // RightHandSide(yIn, DyDx) ; // 1st Step - Not doing, getting passed
182 
183  for(i=0; i<numberOfVariables; ++i)
184  {
185  yTemp[i] = yIn[i] + b21*Step*fInitialDyDx[i] ;
186  }
187  RightHandSide(yTemp, ak2) ; // 2nd Step
188 
189  for(i=0; i<numberOfVariables; ++i)
190  {
191  yTemp[i] = yIn[i] + Step*(b31*fInitialDyDx[i] + b32*ak2[i]) ;
192  }
193  RightHandSide(yTemp, ak3) ; // 3rd Step
194 
195  for(i=0; i<numberOfVariables; ++i)
196  {
197  yTemp[i] = yIn[i] + Step*(b41*fInitialDyDx[i]
198  + b42*ak2[i] + b43*ak3[i]) ;
199  }
200  RightHandSide(yTemp, ak4) ; // 4th Step
201 
202  for(i=0; i<numberOfVariables; ++i)
203  {
204  yTemp[i] = yIn[i] + Step*(b51*fInitialDyDx[i]
205  + b52*ak2[i] + b53*ak3[i] + b54*ak4[i]) ;
206  }
207  RightHandSide(yTemp, ak5) ; // 5th Step
208 
209  for(i=0; i<numberOfVariables; ++i)
210  {
211  yTemp[i] = yIn[i] + Step*(b61*fInitialDyDx[i] + b62*ak2[i]
212  + b63*ak3[i] + b64*ak4[i] + b65*ak5[i]) ;
213  }
214  RightHandSide(yTemp, ak6) ; // 6th Step
215 
216  for(i=0; i<numberOfVariables; ++i)
217  {
218  yOut[i] = yIn[i] + Step*(b71*fInitialDyDx[i] + b72*ak2[i] + b73*ak3[i]
219  + b74*ak4[i] + b75*ak5[i] + b76*ak6[i]);
220  }
221  RightHandSide(yOut, ak7); //7th and Final step
222 
223  for(i=0; i<numberOfVariables; ++i)
224  {
225 
226  yErr[i] = Step*(dc1*fInitialDyDx[i] + dc2*ak2[i] + dc3*ak3[i]
227  + dc4*ak4[i] + dc5*ak5[i] + dc6*ak6[i] + dc7*ak7[i] ) ;
228 
229  // Store Input and Final values, for possible use in calculating chord
230  //
231  fLastInitialVector[i] = yIn[i] ;
232  fLastFinalVector[i] = yOut[i];
233  fLastDyDx[i] = fInitialDyDx[i];
234  nextDydx[i] = ak7[i];
235  }
236  fLastStepLength = Step;
237 
238  return ;
239 }
240 
241 // DistChord
242 //
243 G4double G4FSALDormandPrince745::DistChord() const
244 {
245  G4double distLine, distChord;
246  G4ThreeVector initialPoint, finalPoint, midPoint;
247 
248  // Store last initial and final points
249  // (they will be overwritten in self-Stepper call!)
250  //
251  initialPoint = G4ThreeVector(fLastInitialVector[0],
252  fLastInitialVector[1], fLastInitialVector[2]);
253  finalPoint = G4ThreeVector(fLastFinalVector[0],
254  fLastFinalVector[1], fLastFinalVector[2]);
255 
256  // Do half a step using StepNoErr
257 
258  fAuxStepper->Stepper( fLastInitialVector, fLastDyDx, 0.5 * fLastStepLength,
259  fMidVector, fMidError, pseudoDydx_for_DistChord );
260 
261  midPoint = G4ThreeVector( fMidVector[0], fMidVector[1], fMidVector[2] );
262 
263  // Use stored values of Initial and Endpoint + new Midpoint to evaluate
264  // distance of Chord
265  //
266  if (initialPoint != finalPoint)
267  {
268  distLine = G4LineSection::Distline( midPoint,initialPoint,finalPoint );
269  distChord = distLine;
270  }
271  else
272  {
273  distChord = (midPoint-initialPoint).mag();
274  }
275  return distChord;
276 }
277 
278 // interpolate
279 //
280 void G4FSALDormandPrince745::interpolate( const G4double yInput[],
281  const G4double dydx[],
282  G4double yOut[],
283  G4double Step,
284  G4double tau)
285 {
286  G4double bf1, bf2, bf3, bf4, bf5, bf6, bf7;
287 
288  const G4int numberOfVariables = GetNumberOfVariables();
289 
290  G4double tau0 = tau;
291 
292  for(G4int i=0;i<numberOfVariables; ++i)
293  {
294  yIn[i]=yInput[i];
295  }
296 
297  G4double tau_2 = tau0*tau0 ,
298  tau_3 = tau0*tau_2,
299  tau_4 = tau_2*tau_2;
300 
301  bf1 = (157015080.0*tau_4 - 13107642775.0*tau_3
302  + 34969693132.0*tau_2- 32272833064.0*tau + 11282082432.0)
303  / 11282082432.0;
304  bf2 = 0.0;
305  bf3 = - 100.0*tau*(15701508.0*tau_3 - 914128567.0*tau_2
306  + 2074956840.0*tau - 1323431896.0) / 32700410799.0;
307  bf4 = 25.0*tau*(94209048.0*tau_3- 1518414297.0*tau_2
308  + 2460397220.0*tau - 889289856.0)
309  / 5641041216.0;
310  bf5 = -2187.0*tau*(52338360.0*tau_3 - 451824525.0*tau_2
311  + 687873124.0*tau - 259006536.0)
312  / 199316789632.0;
313  bf6 = 11.0*tau*(106151040.0*tau_3- 661884105.0*tau_2
314  + 946554244.0*tau - 361440756.0)
315  / 2467955532.0;
316  bf7 = tau*(1.0 - tau)*(8293050.0*tau_2 - 82437520.0*tau + 44764047.0)
317  / 29380423.0;
318 
319  for(G4int i=0; i<numberOfVariables; ++i)
320  {
321  yOut[i] = yIn[i] + Step*tau*(bf1*dydx[i] + bf2*ak2[i] + bf3*ak3[i]
322  + bf4*ak4[i] + bf5*ak5[i] + bf6*ak6[i]
323  + bf7*ak7[i] );
324  }
325 }
326 
327 // SetupInterpolate
328 //
329 void G4FSALDormandPrince745::SetupInterpolate(const G4double yInput[],
330  const G4double dydx[],
331  const G4double Step )
332 {
333  // Coefficients for the additional stages
334  //
335  G4double b81 = 6245.0/62208.0 ,
336  b82 = 0.0 ,
337  b83 = 8875.0/103032.0 ,
338  b84 = -125.0/1728.0 ,
339  b85 = 801.0/13568.0 ,
340  b86 = -13519.0/368064.0 ,
341  b87 = 11105.0/368064.0 ,
342 
343  b91 = 632855.0/4478976.0 ,
344  b92 = 0.0 ,
345  b93 = 4146875.0/6491016.0 ,
346  b94 = 5490625.0/14183424.0 ,
347  b95 = -15975.0/108544.0 ,
348  b96 = 8295925.0/220286304.0 ,
349  b97 = -1779595.0/62938944.0 ,
350  b98 = -805.0/4104.0 ;
351 
352  const G4int numberOfVariables = GetNumberOfVariables();
353 
354  // Saving yInput because yInput and yOut can be aliases for same array
355  //
356  for(G4int i=0; i<numberOfVariables; ++i)
357  {
358  yIn[i] = yInput[i];
359  }
360 
361  yTemp[7] = yIn[7];
362 
363  // Evaluate the extra stages
364  //
365  for(G4int i=0; i<numberOfVariables; ++i)
366  {
367  yTemp[i] = yIn[i] + Step*( b81*dydx[i] + b82*ak2[i] + b83*ak3[i] +
368  b84*ak4[i] + b85*ak5[i] + b86*ak6[i] +
369  b87*ak7[i] );
370  }
371  RightHandSide( yTemp, ak8 ); // 8th Stage
372 
373  for(G4int i=0; i<numberOfVariables; ++i)
374  {
375  yTemp[i] = yIn[i] + Step * ( b91*dydx[i] + b92*ak2[i] + b93*ak3[i] +
376  b94*ak4[i] + b95*ak5[i] + b96*ak6[i] +
377  b97*ak7[i] + b98*ak8[i] );
378  }
379  RightHandSide( yTemp, ak9 ); // 9th Stage
380 }
381 
382 // Interpolate
383 //
384 void G4FSALDormandPrince745::Interpolate( const G4double yInput[],
385  const G4double dydx[],
386  const G4double Step,
387  G4double yOut[],
388  G4double tau )
389 {
390  // Define the coefficients for the polynomials
391 
392  G4double bi[10][5], b[10];
393  G4int numberOfVariables = GetNumberOfVariables();
394 
395  // COEFFICIENTS OF bi[1]
396  bi[1][0] = 1.0 ,
397  bi[1][1] = -38039.0/7040.0 ,
398  bi[1][2] = 125923.0/10560.0 ,
399  bi[1][3] = -19683.0/1760.0 ,
400  bi[1][4] = 3303.0/880.0 ,
401  // --------------------------------------------------------
402  //
403  // COEFFICIENTS OF bi[2]
404  bi[2][0] = 0.0 ,
405  bi[2][1] = 0.0 ,
406  bi[2][2] = 0.0 ,
407  bi[2][3] = 0.0 ,
408  bi[2][4] = 0.0 ,
409  // --------------------------------------------------------
410  //
411  // COEFFICIENTS OF bi[3]
412  bi[3][0] = 0.0 ,
413  bi[3][1] = -12500.0/4081.0 ,
414  bi[3][2] = 205000.0/12243.0 ,
415  bi[3][3] = -90000.0/4081.0 ,
416  bi[3][4] = 36000.0/4081.0 ,
417  // --------------------------------------------------------
418  //
419  // COEFFICIENTS OF bi[4]
420  bi[4][0] = 0.0 ,
421  bi[4][1] = -3125.0/704.0 ,
422  bi[4][2] = 25625.0/1056.0 ,
423  bi[4][3] = -5625.0/176.0 ,
424  bi[4][4] = 1125.0/88.0 ,
425  // --------------------------------------------------------
426  //
427  // COEFFICIENTS OF bi[5]
428  bi[5][0] = 0.0 ,
429  bi[5][1] = 164025.0/74624.0 ,
430  bi[5][2] = -448335.0/37312.0 ,
431  bi[5][3] = 295245.0/18656.0 ,
432  bi[5][4] = -59049.0/9328.0 ,
433  // --------------------------------------------------------
434  //
435  // COEFFICIENTS OF bi[6]
436  bi[6][0] = 0.0 ,
437  bi[6][1] = -25.0/28.0 ,
438  bi[6][2] = 205.0/42.0 ,
439  bi[6][3] = -45.0/7.0 ,
440  bi[6][4] = 18.0/7.0 ,
441  // --------------------------------------------------------
442  //
443  // COEFFICIENTS OF bi[7]
444  bi[7][0] = 0.0 ,
445  bi[7][1] = -2.0/11.0 ,
446  bi[7][2] = 73.0/55.0 ,
447  bi[7][3] = -171.0/55.0 ,
448  bi[7][4] = 108.0/55.0 ,
449  // --------------------------------------------------------
450  //
451  // COEFFICIENTS OF bi[8]
452  bi[8][0] = 0.0 ,
453  bi[8][1] = 189.0/22.0 ,
454  bi[8][2] = -1593.0/55.0 ,
455  bi[8][3] = 3537.0/110.0 ,
456  bi[8][4] = -648.0/55.0 ,
457  // --------------------------------------------------------
458  //
459  // COEFFICIENTS OF bi[9]
460  bi[9][0] = 0.0 ,
461  bi[9][1] = 351.0/110.0 ,
462  bi[9][2] = -999.0/55.0 ,
463  bi[9][3] = 2943.0/110.0 ,
464  bi[9][4] = -648.0/55.0 ;
465  // --------------------------------------------------------
466 
467  for(G4int i = 0; i< numberOfVariables; ++i)
468  {
469  yIn[i] = yInput[i];
470  }
471 
472  G4double tau0 = tau;
473 
474  // Calculating the polynomials
475  //
476  for(auto i=1; i<=9; ++i) // i is NOT the coordinate no., it's stage no.
477  {
478  b[i] = 0;
479  tau = 1.0;
480  for(auto j=0; j<=4; ++j)
481  {
482  b[i] += bi[i][j]*tau;
483  tau*=tau0;
484  }
485  }
486 
487  for(G4int i=0; i<numberOfVariables; ++i) // Here i IS the coordinate no.
488  {
489  yOut[i] = yIn[i] + Step*tau0*(b[1]*dydx[i] + b[2]*ak2[i] + b[3]*ak3[i] +
490  b[4]*ak4[i] + b[5]*ak5[i] + b[6]*ak6[i] +
491  b[7]*ak7[i] + b[8]*ak8[i] + b[9]*ak9[i] );
492  }
493 }