Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/geometry/magneticfield/src/G4DormandPrince745.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 // 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 
 60 // Name of this steppers
 61 const G4String& G4DormandPrince745::StepperType() const
 62 {
 63   static G4String _stepperType("G4DormandPrince745: 5th order");
 64   return _stepperType;
 65 }
 66 
 67 // Description of this steppers - plus details of its implementation
 68 const G4String& G4DormandPrince745::StepperDescription() const
 69 {
 70   static G4String _stepperDescription(
 71     "Embedeed 5th order Runge-Kutta stepper - 7 stages, FSAL, Interpolating.");
 72   return _stepperDescription;
 73 }
 74 
 75 G4DormandPrince745::G4DormandPrince745(G4EquationOfMotion* equation,
 76                                        G4int noIntegrationVariables)
 77     : G4MagIntegratorStepper(equation, noIntegrationVariables)
 78 {
 79 }
 80 
 81 void G4DormandPrince745::Stepper(const G4double yInput[],
 82                                  const G4double dydx[],
 83                                        G4double hstep,
 84                                        G4double yOutput[],
 85                                        G4double yError[],
 86                                        G4double dydxOutput[])
 87 {
 88   Stepper(yInput, dydx, hstep, yOutput, yError);
 89   field_utils::copy(dydxOutput, ak7);
 90 }
 91 
 92 // Stepper
 93 //
 94 // Passing in the value of yInput[],the first time dydx[] and Step length
 95 // Giving back yOut and yErr arrays for output and error respectively
 96 //
 97 void G4DormandPrince745::Stepper(const G4double yInput[],
 98                                  const G4double dydx[],
 99                                        G4double hstep,
100                                        G4double yOut[],
101                                        G4double yErr[])
102 {
103     // The various constants defined on the basis of butcher tableu
104     //
105     const G4double b21 = 0.2,
106                    b31 = 3.0 / 40.0, b32 = 9.0 / 40.0,
107                    b41 = 44.0 / 45.0, b42 = -56.0 / 15.0, b43 = 32.0/9.0,
108 
109         b51 = 19372.0 / 6561.0, b52 = -25360.0 / 2187.0, b53 = 64448.0 / 6561.0,
110         b54 = -212.0 / 729.0,
111         
112         b61 = 9017.0 / 3168.0 , b62 =   -355.0 / 33.0,
113         b63 = 46732.0 / 5247.0, b64 = 49.0 / 176.0,
114         b65 = -5103.0 / 18656.0,
115         
116         b71 = 35.0 / 384.0, b72 = 0.,
117         b73 = 500.0 / 1113.0, b74 = 125.0 / 192.0,
118         b75 = -2187.0 / 6784.0, b76 = 11.0 / 84.0,
119     
120     //Sum of columns, sum(bij) = ei
121     //    e1 = 0. ,
122     //    e2 = 1.0/5.0 ,
123     //    e3 = 3.0/10.0 ,
124     //    e4 = 4.0/5.0 ,
125     //    e5 = 8.0/9.0 ,
126     //    e6 = 1.0 ,
127     //    e7 = 1.0 ,
128     
129     // Difference between the higher and the lower order method coeff. :
130     // b7j are the coefficients of higher order
131     
132         dc1 = -(b71 - 5179.0 / 57600.0),
133         dc2 = -(b72 - .0),
134         dc3 = -(b73 - 7571.0 / 16695.0),
135         dc4 = -(b74 - 393.0 / 640.0),
136         dc5 = -(b75 + 92097.0 / 339200.0),
137         dc6 = -(b76 - 187.0 / 2100.0),
138         dc7 = -(- 1.0 / 40.0);
139     
140     const G4int numberOfVariables = GetNumberOfVariables();
141     State yTemp = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
142     
143     // The number of variables to be integrated over
144     //
145     yOut[7] = yTemp[7]  = yInput[7];
146 
147     //  Saving yInput because yInput and yOut can be aliases for same array
148     //
149     for(G4int i = 0; i < numberOfVariables; ++i)
150     {
151         fyIn[i] = yInput[i];
152     }
153     // RightHandSide(yIn, dydx);    // Not done! 1st stage
154     
155     for(G4int i = 0; i < numberOfVariables; ++i)
156     {
157         yTemp[i] = fyIn[i] + b21 * hstep * dydx[i];
158     }
159     RightHandSide(yTemp, ak2);              // 2nd stage
160     
161     for(G4int i = 0; i < numberOfVariables; ++i)
162     {
163         yTemp[i] = fyIn[i] + hstep * (b31 * dydx[i] + b32 * ak2[i]);
164     }
165     RightHandSide(yTemp, ak3);              // 3rd stage
166     
167     for(G4int i = 0; i < numberOfVariables; ++i)
168     {
169         yTemp[i] = fyIn[i] + hstep * (
170             b41 * dydx[i] + b42 * ak2[i] + b43 * ak3[i]);
171     }
172     RightHandSide(yTemp, ak4);              // 4th stage
173     
174     for(G4int i = 0; i < numberOfVariables; ++i)
175     {
176         yTemp[i] = fyIn[i] + hstep * (
177             b51 * dydx[i] + b52 * ak2[i] + b53 * ak3[i] + b54 * ak4[i]);
178     }
179     RightHandSide(yTemp, ak5);              // 5th stage
180     
181     for(G4int i = 0; i < numberOfVariables; ++i)
182     {
183         yTemp[i] = fyIn[i] + hstep * (
184             b61 * dydx[i] + b62 * ak2[i] + 
185             b63 * ak3[i] + b64 * ak4[i] + b65 * ak5[i]);
186     }
187     RightHandSide(yTemp, ak6);              // 6th stage
188     
189     for(G4int i = 0; i < numberOfVariables; ++i)
190     {
191         yOut[i] = fyIn[i] + hstep * (
192             b71 * dydx[i] + b72 * ak2[i] + b73 * ak3[i] +
193             b74 * ak4[i] + b75 * ak5[i] + b76 * ak6[i]);
194     }
195     RightHandSide(yOut, ak7);               // 7th and Final stage
196     
197     for(G4int i = 0; i < numberOfVariables; ++i)
198     {
199         yErr[i] = hstep * (
200             dc1 * dydx[i] + dc2 * ak2[i] + 
201             dc3 * ak3[i] + dc4 * ak4[i] +
202             dc5 * ak5[i] + dc6 * ak6[i] + dc7 * ak7[i]
203         ) + 1.5e-18;
204 
205         // Store Input and Final values, for possible use in calculating chord
206         //
207         fyOut[i] = yOut[i];
208         fdydxIn[i] = dydx[i];        
209     }
210     
211     fLastStepLength = hstep;
212 }
213 
214 G4double G4DormandPrince745::DistChord() const
215 {
216     // Coefficients were taken from Some Practical Runge-Kutta Formulas
217     // by Lawrence F. Shampine, page 149, c*
218     //
219     const G4double hf1 = 6025192743.0 / 30085553152.0,
220                    hf3 = 51252292925.0 / 65400821598.0,
221                    hf4 = - 2691868925.0 / 45128329728.0,
222                    hf5 = 187940372067.0 / 1594534317056.0,
223                    hf6 = - 1776094331.0 / 19743644256.0,
224                    hf7 = 11237099.0 / 235043384.0;
225 
226     G4ThreeVector mid;
227 
228     for(G4int i = 0; i < 3; ++i) 
229     {
230         mid[i] = fyIn[i] + 0.5 * fLastStepLength * (
231             hf1 * fdydxIn[i] + hf3 * ak3[i] + 
232             hf4 * ak4[i] + hf5 * ak5[i] + hf6 * ak6[i] + hf7 * ak7[i]);
233     }
234 
235     const G4ThreeVector begin = makeVector(fyIn, Value3D::Position);
236     const G4ThreeVector end = makeVector(fyOut, Value3D::Position);
237 
238     return G4LineSection::Distline(mid, begin, end);
239 }
240 
241 // The lower (4th) order interpolant given by Dormand and Prince:
242 //        J. R. Dormand and P. J. Prince, "Runge-Kutta triples"
243 //        Computers & Mathematics with Applications, vol. 12, no. 9,
244 //        pp. 1007-1017, 1986.
245 //
246 void G4DormandPrince745::
247 Interpolate4thOrder(G4double yOut[], G4double tau) const
248 {
249     const G4int numberOfVariables = GetNumberOfVariables();
250     
251     const G4double tau2 = tau * tau,
252                    tau3 = tau * tau2,
253                    tau4 = tau2 * tau2;
254 
255     const G4double bf1 = 1.0 / 11282082432.0 * (
256         157015080.0 * tau4 - 13107642775.0 * tau3 + 34969693132.0 * tau2 - 
257         32272833064.0 * tau + 11282082432.0);
258 
259     const G4double bf3 = - 100.0 / 32700410799.0 * tau * (
260         15701508.0 * tau3 - 914128567.0 * tau2 + 2074956840.0 * tau - 
261         1323431896.0);
262     
263     const G4double bf4 = 25.0 / 5641041216.0 * tau * (
264         94209048.0 * tau3 - 1518414297.0 * tau2 + 2460397220.0 * tau - 
265         889289856.0);
266     
267     const G4double bf5 = - 2187.0 / 199316789632.0 * tau * (
268         52338360.0 * tau3 - 451824525.0 * tau2 + 687873124.0 * tau - 
269         259006536.0);
270 
271     const G4double bf6 = 11.0 / 2467955532.0 * tau * (
272         106151040.0 * tau3 - 661884105.0 * tau2 + 
273         946554244.0 * tau - 361440756.0);
274     
275     const G4double bf7 = 1.0 / 29380423.0 * tau * (1.0 - tau) * (
276         8293050.0 * tau2 - 82437520.0 * tau + 44764047.0);
277 
278     for(G4int i = 0; i < numberOfVariables; ++i)
279     {
280         yOut[i] = fyIn[i] + fLastStepLength * tau * (
281             bf1 * fdydxIn[i] + bf3 * ak3[i] + bf4 * ak4[i] + 
282             bf5 * ak5[i] + bf6 * ak6[i] + bf7 * ak7[i]);
283     }
284 }
285 
286 // Following interpolant of order 5 was given by Baker,Dormand,Gilmore, Prince :
287 //        T. S. Baker, J. R. Dormand, J. P. Gilmore, and P. J. Prince,
288 //        "Continuous approximation with embedded Runge-Kutta methods"
289 //        Applied Numerical Mathematics, vol. 22, no. 1, pp. 51-62, 1996.
290 //
291 // Calculating the extra stages for the interpolant
292 //
293 void G4DormandPrince745::SetupInterpolation5thOrder()
294 {
295     // Coefficients for the additional stages
296     //
297     const G4double b81 = 6245.0 / 62208.0,
298                    b82 = 0.0,
299                    b83 = 8875.0 / 103032.0,
300                    b84 = -125.0 / 1728.0,
301                    b85 = 801.0 / 13568.0,
302                    b86 = -13519.0 / 368064.0,
303                    b87 = 11105.0 / 368064.0,
304         
305                    b91 = 632855.0 / 4478976.0,
306                    b92 = 0.0,
307                    b93 = 4146875.0 / 6491016.0,
308                    b94 = 5490625.0 /14183424.0,
309                    b95 = -15975.0 / 108544.0,
310                    b96 = 8295925.0 / 220286304.0,
311                    b97 = -1779595.0 / 62938944.0,
312                    b98 = -805.0 / 4104.0;
313     
314     const G4int numberOfVariables = GetNumberOfVariables();
315     State yTemp = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
316 
317     // Evaluate the extra stages
318     //
319     for(G4int i = 0; i < numberOfVariables; ++i)
320     {
321         yTemp[i] = fyIn[i] + fLastStepLength * (
322             b81 * fdydxIn[i] + b82 * ak2[i] + b83 * ak3[i] +
323             b84 * ak4[i] + b85 * ak5[i] + b86 * ak6[i] +
324             b87 * ak7[i]
325         );
326     }
327     RightHandSide(yTemp, ak8);          // 8th Stage
328     
329     for(G4int i = 0; i < numberOfVariables; ++i)
330     {
331         yTemp[i] = fyIn[i] + fLastStepLength * (
332             b91 * fdydxIn[i] + b92 * ak2[i] + b93 * ak3[i] +
333             b94 * ak4[i] + b95 * ak5[i] + b96 * ak6[i] +
334             b97 * ak7[i] + b98 * ak8[i]
335         );
336     }
337     RightHandSide(yTemp, ak9);          // 9th Stage
338 }
339 
340 // Calculating the interpolated result yOut with the coefficients
341 //
342 void G4DormandPrince745::
343 Interpolate5thOrder(G4double yOut[], G4double tau) const
344 {
345     // Define the coefficients for the polynomials
346     //
347     G4double bi[10][5];
348     
349     //  COEFFICIENTS OF   bi[1]
350     bi[1][0] = 1.0,
351     bi[1][1] = -38039.0 / 7040.0,
352     bi[1][2] = 125923.0 / 10560.0,
353     bi[1][3] = -19683.0 / 1760.0,
354     bi[1][4] = 3303.0 / 880.0,
355     //  --------------------------------------------------------
356     //
357     //  COEFFICIENTS OF  bi[2]
358     bi[2][0] = 0.0,
359     bi[2][1] = 0.0,
360     bi[2][2] = 0.0,
361     bi[2][3] = 0.0,
362     bi[2][4] = 0.0,
363     //  --------------------------------------------------------
364     //
365     //  COEFFICIENTS OF  bi[3]
366     bi[3][0] = 0.0,
367     bi[3][1] = -12500.0 / 4081.0,
368     bi[3][2] = 205000.0 / 12243.0,
369     bi[3][3] = -90000.0 / 4081.0,
370     bi[3][4] = 36000.0 / 4081.0,
371     //  --------------------------------------------------------
372     //
373     //  COEFFICIENTS OF  bi[4]
374     bi[4][0] = 0.0,
375     bi[4][1] = -3125.0 / 704.0,
376     bi[4][2] = 25625.0 / 1056.0,
377     bi[4][3] = -5625.0 / 176.0,
378     bi[4][4] = 1125.0 / 88.0,
379     //  --------------------------------------------------------
380     //
381     //  COEFFICIENTS OF  bi[5]
382     bi[5][0] = 0.0,
383     bi[5][1] = 164025.0 / 74624.0,
384     bi[5][2] = -448335.0 / 37312.0,
385     bi[5][3] = 295245.0 / 18656.0,
386     bi[5][4] = -59049.0 / 9328.0,
387     //  --------------------------------------------------------
388     //
389     //  COEFFICIENTS OF  bi[6]
390     bi[6][0] = 0.0,
391     bi[6][1] = -25.0 / 28.0,
392     bi[6][2] = 205.0 / 42.0,
393     bi[6][3] = -45.0 / 7.0,
394     bi[6][4] = 18.0 / 7.0,
395     //  --------------------------------------------------------
396     //
397     //  COEFFICIENTS OF  bi[7]
398     bi[7][0] = 0.0,
399     bi[7][1] = -2.0 / 11.0,
400     bi[7][2] = 73.0 / 55.0,
401     bi[7][3] = -171.0 / 55.0,
402     bi[7][4] = 108.0 / 55.0,
403     //  --------------------------------------------------------
404     //
405     //  COEFFICIENTS OF  bi[8]
406     bi[8][0] = 0.0,
407     bi[8][1] = 189.0 / 22.0,
408     bi[8][2] = -1593.0 / 55.0,
409     bi[8][3] = 3537.0 / 110.0,
410     bi[8][4] = -648.0 / 55.0,
411     //  --------------------------------------------------------
412     //
413     //  COEFFICIENTS OF  bi[9]
414     bi[9][0] = 0.0,
415     bi[9][1] = 351.0 / 110.0,
416     bi[9][2] = -999.0 / 55.0,
417     bi[9][3] = 2943.0 / 110.0,
418     bi[9][4] = -648.0 / 55.0;
419     //  --------------------------------------------------------
420         
421     // Calculating the polynomials
422  
423     G4double b[10];
424     std::memset(b, 0.0, sizeof(b));
425 
426     G4double tauPower = 1.0;
427     for(G4int j = 0; j <= 4; ++j)
428     {       
429        for(G4int iStage = 1; iStage <= 9; ++iStage)
430        {
431           b[iStage] += bi[iStage][j] * tauPower;
432        }
433        tauPower *= tau;       
434     }
435     
436     const G4int numberOfVariables = GetNumberOfVariables();
437     const G4double stepLen = fLastStepLength * tau;
438     for(G4int i = 0; i < numberOfVariables; ++i)
439     {
440         yOut[i] = fyIn[i] + stepLen * (
441             b[1] * fdydxIn[i] + b[2] * ak2[i] + b[3] * ak3[i] +
442             b[4] * ak4[i] + b[5] * ak5[i] + b[6] * ak6[i] +
443             b[7] * ak7[i] + b[8] * ak8[i] + b[9] * ak9[i]
444         );
445     }
446 }
447