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