ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4FSALBogackiShampine45.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4FSALBogackiShampine45.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 // G4FSALBogackiShampine45 implementation
27 //
28 // The Butcher table of the Bogacki-Shampine-8-4-5 method is as follows:
29 //
30 // 0 |
31 // 1/6 | 1/6
32 // 2/9 | 2/27 4/27
33 // 3/7 | 183/1372 -162/343 1053/1372
34 // 2/3 | 68/297 -4/11 42/143 1960/3861
35 // 3/4 | 597/22528 81/352 63099/585728 58653/366080 4617/20480
36 // 1 | 174197/959244 -30942/79937 8152137/19744439 666106/1039181 -29421/29068 482048/414219
37 // 1 | 587/8064 0 4440339/15491840 24353/124800 387/44800 2152/5985 7267/94080
38 // -------------------------------------------------------------------------------------------------------------------
39 // 587/8064 0 4440339/15491840 24353/124800 387/44800 2152/5985 7267/94080 0
40 // 2479/34992 0 123/416 612941/3411720 43/1440 2272/6561 79937/1113912 3293/556956
41 //
42 // Created: Somnath Banerjee, Google Summer of Code 2015, 26 May 2015
43 // Supervision: John Apostolakis, CERN
44 // --------------------------------------------------------------------
45 
46 // Plan is that this source file / class will be merged with the updated
47 // BogackiShampine45 class, which contains improvements (May 2016)
48 
49 #include <cassert>
50 
51 #include "G4FSALBogackiShampine45.hh"
52 #include "G4LineSection.hh"
53 
54 G4bool G4FSALBogackiShampine45::fPreparedConstants = false;
55 G4double G4FSALBogackiShampine45::bi[12][7];
56 
57 // Constructor
58 //
59 G4FSALBogackiShampine45::G4FSALBogackiShampine45(G4EquationOfMotion* EqRhs,
60  G4int noIntegrationVariables,
61  G4bool primary)
62  : G4VFSALIntegrationStepper(EqRhs, noIntegrationVariables)
63 {
64  const G4int numberOfVariables = noIntegrationVariables;
65 
66  // New Chunk of memory being created for use by the stepper
67 
68  // aki - for storing intermediate RHS
69  //
70  ak2 = new G4double[numberOfVariables];
71  ak3 = new G4double[numberOfVariables];
72  ak4 = new G4double[numberOfVariables];
73  ak5 = new G4double[numberOfVariables];
74  ak6 = new G4double[numberOfVariables];
75  ak7 = new G4double[numberOfVariables];
76  ak8 = new G4double[numberOfVariables];
77 
78  ak9 = new G4double[numberOfVariables];
79  ak10 = new G4double[numberOfVariables];
80  ak11 = new G4double[numberOfVariables];
81  DyDx = new G4double[numberOfVariables];
82 
83  assert ( GetNumberOfStateVariables() >= 8 );
84  const G4int numStateVars = std::max(noIntegrationVariables,
85  GetNumberOfStateVariables() );
86 
87  // Must ensure space extra 'state' variables exists - i.e. yIn[7]
88  //
89  yTemp = new G4double[numStateVars];
90  yIn = new G4double[numStateVars] ;
91 
92  fLastInitialVector = new G4double[numStateVars] ;
93  fLastFinalVector = new G4double[numStateVars] ;
94  fLastDyDx = new G4double[numberOfVariables]; // Only derivatives
95 
96  fMidVector = new G4double[numStateVars];
97  fMidError = new G4double[numStateVars];
98 
99  pseudoDydx_for_DistChord = new G4double[numberOfVariables];
100 
101  fMidVector = new G4double[numberOfVariables];
102  fMidError = new G4double[numberOfVariables];
103  if( primary )
104  {
105  fAuxStepper = new G4FSALBogackiShampine45(EqRhs, numberOfVariables,
106  !primary);
107  }
108  if( !fPreparedConstants )
109  {
110  PrepareConstants();
111  }
112 }
113 
114 // Destructor
115 //
116 G4FSALBogackiShampine45::~G4FSALBogackiShampine45()
117 {
118  // Clear all previously allocated memory for stepper and DistChord
119 
120  delete [] ak2;
121  delete [] ak3;
122  delete [] ak4;
123  delete [] ak5;
124  delete [] ak6;
125  delete [] ak7;
126  delete [] ak8;
127  delete [] ak9;
128  delete [] ak10;
129  delete [] ak11;
130  delete [] DyDx;
131  delete [] yTemp;
132  delete [] yIn;
133 
134  delete [] fLastInitialVector;
135  delete [] fLastFinalVector;
136  delete [] fLastDyDx;
137  delete [] fMidVector;
138  delete [] fMidError;
139 
140  delete fAuxStepper;
141 
142  delete [] pseudoDydx_for_DistChord;
143 }
144 
145 // Stepper
146 //
147 // Passing in the value of yInput[],the first time dydx[] and Step length
148 // Giving back yOut and yErr arrays for output and error respectively
149 //
150 void G4FSALBogackiShampine45::Stepper(const G4double yInput[],
151  const G4double dydx[],
152  G4double Step,
153  G4double yOut[],
154  G4double yErr[],
155  G4double nextDydx[])
156 {
157  G4int i;
158 
159  // The various constants defined on the basis of butcher tableu
160 
161  const G4double b21 = 1.0/6.0 ,
162  b31 = 2.0/27.0 , b32 = 4.0/27.0,
163 
164  b41 = 183.0/1372.0 , b42 = -162.0/343.0, b43 = 1053.0/1372.0,
165 
166  b51 = 68.0/297.0, b52 = -4.0/11.0,
167  b53 = 42.0/143.0, b54 = 1960.0/3861.0,
168 
169  b61 = 597.0/22528.0, b62 = 81.0/352.0,
170  b63 = 63099.0/585728.0, b64 = 58653.0/366080.0,
171  b65 = 4617.0/20480.0,
172 
173  b71 = 174197.0/959244.0, b72 = -30942.0/79937.0,
174  b73 = 8152137.0/19744439.0, b74 = 666106.0/1039181.0,
175  b75 = -29421.0/29068.0, b76 = 482048.0/414219.0,
176 
177  b81 = 587.0/8064.0, b82 = 0.0,
178  b83 = 4440339.0/15491840.0, b84 = 24353.0/124800.0,
179  b85 = 387.0/44800.0, b86 = 2152.0/5985.0,
180  b87 = 7267.0/94080.0,
181 
182 
183  // c1 = 2479.0/34992.0,
184  // c2 = 0.0,
185  // c3 = 123.0/416.0,
186  // c4 = 612941.0/3411720.0,
187  // c5 = 43.0/1440.0,
188  // c6 = 2272.0/6561.0,
189  // c7 = 79937.0/1113912.0,
190  // c8 = 3293.0/556956.0,
191 
192  // For the embedded higher order method only the difference of values
193  // taken and is used directly later instead of defining the last row
194  // of butcher table in a separate set of variables and taking the
195  // difference there
196 
197  dc1 = b81 - 2479.0/34992.0 ,
198  dc2 = 0.0,
199  dc3 = b83 - 123.0/416.0 ,
200  dc4 = b84 - 612941.0/3411720.0,
201  dc5 = b85 - 43.0/1440.0,
202  dc6 = b86 - 2272.0/6561.0,
203  dc7 = b87 - 79937.0/1113912.0,
204  dc8 = -3293.0/556956.0; // end of declaration
205 
206  const G4int numberOfVariables = GetNumberOfVariables();
207 
208  // The number of variables to be integrated over
209  //
210  yOut[7] = yTemp[7] = yIn[7];
211 
212  // Saving yInput because yInput and yOut can be aliases for same array
213  //
214  for(i=0; i<numberOfVariables; ++i)
215  {
216  yIn[i]=yInput[i];
217  DyDx[i] = dydx[i];
218  }
219  // RightHandSide(yIn, dydx) ; // 1st Step - Not doing, getting passed
220 
221  for(i=0; i<numberOfVariables; ++i)
222  {
223  yTemp[i] = yIn[i] + b21*Step*DyDx[i] ;
224  }
225  RightHandSide(yTemp, ak2) ; // 2nd Step
226 
227  for(i=0; i<numberOfVariables; ++i)
228  {
229  yTemp[i] = yIn[i] + Step*(b31*DyDx[i] + b32*ak2[i]) ;
230  }
231  RightHandSide(yTemp, ak3) ; // 3rd Step
232 
233  for(i=0; i<numberOfVariables; ++i)
234  {
235  yTemp[i] = yIn[i] + Step*(b41*DyDx[i] + b42*ak2[i] + b43*ak3[i]) ;
236  }
237  RightHandSide(yTemp, ak4) ; // 4th Step
238 
239  for(i=0; i<numberOfVariables; ++i)
240  {
241  yTemp[i] = yIn[i] + Step*(b51*DyDx[i] + b52*ak2[i] + b53*ak3[i] +
242  b54*ak4[i]) ;
243  }
244  RightHandSide(yTemp, ak5) ; // 5th Step
245 
246  for(i=0; i<numberOfVariables; ++i)
247  {
248  yTemp[i] = yIn[i] + Step*(b61*DyDx[i] + b62*ak2[i] + b63*ak3[i] +
249  b64*ak4[i] + b65*ak5[i]) ;
250  }
251  RightHandSide(yTemp, ak6) ; // 6th Step
252 
253  for(i=0; i<numberOfVariables; ++i)
254  {
255  yTemp[i] = yIn[i] + Step*(b71*DyDx[i] + b72*ak2[i] + b73*ak3[i] +
256  b74*ak4[i] + b75*ak5[i] + b76*ak6[i]);
257  }
258  RightHandSide(yTemp, ak7); // 7th Step
259 
260  for(i=0; i<numberOfVariables; ++i)
261  {
262  yOut[i] = yIn[i] + Step*(b81*DyDx[i] + b82*ak2[i] + b83*ak3[i] +
263  b84*ak4[i] + b85*ak5[i] + b86*ak6[i] +
264  b87*ak7[i]);
265  }
266  RightHandSide(yOut, ak8); // 8th Step - Final one Using FSAL
267 
268 
269  for(i=0; i<numberOfVariables; ++i)
270  {
271 
272  yErr[i] = Step*(dc1*DyDx[i] + dc2*ak2[i] + dc3*ak3[i] + dc4*ak4[i] +
273  dc5*ak5[i] + dc6*ak6[i] + dc7*ak7[i] + dc8*ak8[i]) ;
274 
275 
276  // FSAL stepper : Must pass the last DyDx for the next step, here ak8
277  //
278  nextDydx[i] = ak8[i];
279 
280  // Store Input and Final values, for possible use in calculating chord
281  //
282  fLastInitialVector[i] = yIn[i] ;
283  fLastFinalVector[i] = yOut[i];
284  fLastDyDx[i] = DyDx[i];
285  }
286  fLastStepLength = Step;
287 
288  return;
289 }
290 
291 // DistChord
292 //
293 G4double G4FSALBogackiShampine45::DistChord() const
294 {
295  G4double distLine, distChord;
296  G4ThreeVector initialPoint, finalPoint, midPoint;
297 
298  // Store last initial and final points
299  // (they will be overwritten in self-Stepper call!)
300  //
301  initialPoint = G4ThreeVector( fLastInitialVector[0],
302  fLastInitialVector[1], fLastInitialVector[2]);
303  finalPoint = G4ThreeVector( fLastFinalVector[0],
304  fLastFinalVector[1], fLastFinalVector[2]);
305 
306  // Do half a step using StepNoErr
307 
308  fAuxStepper->Stepper( fLastInitialVector, fLastDyDx, 0.5 * fLastStepLength,
309  fMidVector, fMidError, pseudoDydx_for_DistChord );
310 
311  midPoint = G4ThreeVector( fMidVector[0], fMidVector[1], fMidVector[2] );
312 
313  // Use stored values of Initial and Endpoint + new Midpoint to evaluate
314  // distance of Chord
315  //
316  if (initialPoint != finalPoint)
317  {
318  distLine = G4LineSection::Distline(midPoint, initialPoint, finalPoint);
319  distChord = distLine;
320  }
321  else
322  {
323  distChord = (midPoint-initialPoint).mag();
324  }
325  return distChord;
326 }
327 
328 // PrepareConstants
329 //
330 void G4FSALBogackiShampine45::PrepareConstants()
331 {
332  // --------------------------------------------------------
333  // COEFFICIENTS FOR INTERPOLANT bi WITH 11 STAGES
334  // --------------------------------------------------------
335 
336  // Initialise all values of G4double bi[12][7]
337  //
338  for(auto i=1; i<12; ++i)
339  {
340  for(auto j=1; j<7; ++j)
341  {
342  bi[i][j] = 0.0 ;
343  }
344  }
345 
346  bi[1][6] = -12134338393.0/1050809760.0 ,
347  bi[1][5] = -1620741229.0/50038560.0 ,
348  bi[1][4] = -2048058893.0/59875200.0 ,
349  bi[1][3] = -87098480009.0/5254048800.0 ,
350  bi[1][2] = -11513270273.0/3502699200.0 ,
351  //
352  bi[3][6] = -33197340367.0/1218433216.0 ,
353  bi[3][5] = -539868024987.0/6092166080.0 ,
354  bi[3][4] = -39991188681.0/374902528.0 ,
355  bi[3][3] = -69509738227.0/1218433216.0 ,
356  bi[3][2] = -29327744613.0/2436866432.0 ,
357  //
358  bi[4][6] = -284800997201.0/19905339168.0 ,
359  bi[4][5] = -7896875450471.0/165877826400.0 ,
360  bi[4][4] = -333945812879.0/5671036800.0 ,
361  bi[4][3] = -16209923456237.0/497633479200.0 ,
362  bi[4][2] = -2382590741699.0/331755652800.0 ,
363  //
364  bi[5][6] = -540919.0/741312.0 ,
365  bi[5][5] = -103626067.0/43243200.0 ,
366  bi[5][4] = -633779.0/211200.0 ,
367  bi[5][3] = -32406787.0/18532800.0 ,
368  bi[5][2] = -36591193.0/86486400.0 ,
369  //
370  bi[6][6] = 7157998304.0/374350977.0 ,
371  bi[6][5] = 30405842464.0/623918295.0 ,
372  bi[6][4] = 183022264.0/5332635.0 ,
373  bi[6][3] = -3357024032.0/1871754885.0 ,
374  bi[6][2] = -611586736.0/89131185.0 ,
375  //
376  bi[7][6] = -138073.0/9408.0 ,
377  bi[7][5] = -719433.0/15680.0 ,
378  bi[7][4] = -1620541.0/31360.0 ,
379  bi[7][3] = -385151.0/15680.0 ,
380  bi[7][2] = -65403.0/15680.0 ,
381  //
382  bi[8][6] = 1245.0/64.0 ,
383  bi[8][5] = 3991.0/64.0 ,
384  bi[8][4] = 4715.0/64.0 ,
385  bi[8][3] = 2501.0/64.0 ,
386  bi[8][2] = 149.0/16.0 ,
387  bi[8][1] = 1.0 ,
388  //
389  bi[9][6] = 55.0/3.0 ,
390  bi[9][5] = 71.0 ,
391  bi[9][4] = 103.0 ,
392  bi[9][3] = 199.0/3.0 ,
393  bi[9][2] = 16.0 ,
394  //
395  bi[10][6] = -1774004627.0/75810735.0 ,
396  bi[10][5] = -1774004627.0/25270245.0 ,
397  bi[10][4] = -26477681.0/359975.0 ,
398  bi[10][3] = -11411880511.0/379053675.0 ,
399  bi[10][2] = -423642896.0/126351225.0 ,
400  //
401  bi[11][6] = 35.0 ,
402  bi[11][5] = 105.0 ,
403  bi[11][4] = 117.0 ,
404  bi[11][3] = 59.0 ,
405  bi[11][2] = 12.0 ;
406 }
407 
408 // ---------------------------------------------------------------------------------------
409 
410 void G4FSALBogackiShampine45::interpolate( const G4double yInput[],
411  const G4double dydx[],
412  G4double yOut[],
413  G4double Step,
414  G4double tau )
415 {
416  const G4double a91 = 455.0/6144.0 ,
417  a92 = 0.0 ,
418  a93 = 10256301.0/35409920.0 ,
419  a94 = 2307361.0/17971200.0 ,
420  a95 = -387.0/102400.0 ,
421  a96 = 73.0/5130.0 ,
422  a97 = -7267.0/215040.0 ,
423  a98 = 1.0/32.0 ,
424 
425  a101 = -837888343715.0/13176988637184.0 ,
426  a102 = 30409415.0/52955362.0 ,
427  a103 = -48321525963.0/759168069632.0 ,
428  a104 = 8530738453321.0/197654829557760.0 ,
429  a105 = 1361640523001.0/1626788720640.0 ,
430  a106 = -13143060689.0/38604458898.0 ,
431  a107 = 18700221969.0/379584034816.0 ,
432  a108 = -5831595.0/847285792.0 ,
433  a109 = -5183640.0/26477681.0 ,
434 
435  a111 = 98719073263.0/1551965184000.0 ,
436  a112 = 1307.0/123552.0 ,
437  a113 = 4632066559387.0/70181753241600.0 ,
438  a114 = 7828594302389.0/382182512025600.0 ,
439  a115 = 40763687.0/11070259200.0 ,
440  a116 = 34872732407.0/224610586200.0 ,
441  a117 = -2561897.0/30105600.0 ,
442  a118 = 1.0/10.0 ,
443  a119 = -1.0/10.0 ,
444  a1110 = -1403317093.0/11371610250.0 ;
445 
446  const G4int numberOfVariables = GetNumberOfVariables();
447 
448  // Saving yInput because yInput and yOut can be aliases for same array
449  //
450  for(G4int i=0; i<numberOfVariables; ++i)
451  {
452  yIn[i]=yInput[i];
453  }
454 
455  // The number of variables to be integrated over
456  //
457  yOut[7] = yTemp[7] = yIn[7];
458 
459  // Calculating extra stages
460  //
461  for(G4int i=0; i<numberOfVariables; ++i)
462  {
463  yTemp[i] = yIn[i] + Step*(a91*dydx[i] + a92*ak2[i] + a93*ak3[i] +
464  a94*ak4[i] + a95*ak5[i] + a96*ak6[i] +
465  a97*ak7[i] + a98*ak8[i] );
466  }
467 
468  RightHandSide(yTemp, ak9);
469 
470  for(G4int i=0; i<numberOfVariables; ++i)
471  {
472  yTemp[i] = yIn[i] + Step*(a101*dydx[i] + a102*ak2[i] + a103*ak3[i] +
473  a104*ak4[i] + a105*ak5[i] + a106*ak6[i] +
474  a107*ak7[i] + a108*ak8[i] + a109*ak9[i] );
475  }
476 
477  RightHandSide(yTemp, ak10);
478 
479  for(G4int i=0; i<numberOfVariables; ++i)
480  {
481  yTemp[i] = yIn[i] + Step*(a111*dydx[i] + a112*ak2[i] + a113*ak3[i] +
482  a114*ak4[i] + a115*ak5[i] + a116*ak6[i] +
483  a117*ak7[i] + a118*ak8[i] + a119*ak9[i] +
484  a1110*ak10[i] );
485  }
486 
487  RightHandSide(yTemp, ak11);
488 
489  G4double tau0 = tau;
490 
491  // Calculating the polynomials
492  //
493  for(auto i=1; i<=11; ++i) // i is NOT the coordinate no., it's stage no.
494  {
495  b[i] = 0.0;
496  tau = tau0;
497  for(auto j=1; j<=6; ++j)
498  {
499  b[i] += bi[i][j]*tau;
500  tau*=tau0;
501  }
502  }
503 
504  for(G4int i=0; i<numberOfVariables; ++i)
505  {
506  yOut[i] = yIn[i] + Step*(b[1]*dydx[i] + b[2]*ak2[i] + b[3]*ak3[i] +
507  b[4]*ak4[i] + b[5]*ak5[i] + b[6]*ak6[i] +
508  b[7]*ak7[i] + b[8]*ak8[i] + b[9]*ak9[i] +
509  b[10]*ak10[i] + b[11]*ak11[i] );
510  }
511 }