Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/geometry/magneticfield/src/G4DormandPrinceRK78.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  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 }
437