Geant4  10.03.p01
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Public Member Functions | List of all members
G4DormandPrinceRK78 Class Reference

#include <G4DormandPrinceRK78.hh>

Inheritance diagram for G4DormandPrinceRK78:
Inheritance graph
[legend]
Collaboration diagram for G4DormandPrinceRK78:
Collaboration graph
[legend]

Public Member Functions

 G4DormandPrinceRK78 (G4EquationOfMotion *EqRhs, G4int numberOfVariables=6, G4bool primary=true)
 
 ~G4DormandPrinceRK78 ()
 
void Stepper (const G4double y[], const G4double dydx[], G4double h, G4double yout[], G4double yerr[])
 
G4double DistChord () const
 
G4int IntegratorOrder () const
 
- Public Member Functions inherited from G4MagIntegratorStepper
 G4MagIntegratorStepper (G4EquationOfMotion *Equation, G4int numIntegrationVariables, G4int numStateVariables=12, bool isFSAL=false)
 
virtual ~G4MagIntegratorStepper ()
 
virtual void ComputeRightHandSide (const G4double y[], G4double dydx[])
 
void NormaliseTangentVector (G4double vec[6])
 
void NormalisePolarizationVector (G4double vec[12])
 
void RightHandSide (const double y[], double dydx[])
 
G4int GetNumberOfVariables () const
 
G4int GetNumberOfStateVariables () const
 
G4int IntegrationOrder ()
 
G4EquationOfMotion * GetEquationOfMotion ()
 
void SetEquationOfMotion (G4EquationOfMotion *newEquation)
 
unsigned long GetfNoRHSCalls ()
 
void ResetfNORHSCalls ()
 
bool IsFSAL ()
 

Additional Inherited Members

- Protected Member Functions inherited from G4MagIntegratorStepper
void SetIntegrationOrder (int order)
 
void SetFSAL (bool flag=true)
 

Detailed Description

Definition at line 35 of file G4DormandPrinceRK78.hh.

Constructor & Destructor Documentation

G4DormandPrinceRK78::G4DormandPrinceRK78 ( G4EquationOfMotion *  EqRhs,
G4int  numberOfVariables = 6,
G4bool  primary = true 
)

Definition at line 49 of file G4DormandPrinceRK78.cc.

52 : G4MagIntegratorStepper(EqRhs, noIntegrationVariables)
53 {
54  const G4int numberOfVariables = noIntegrationVariables;
55 
56  //New Chunk of memory being created for use by the stepper
57 
58  //aki - for storing intermediate RHS
59  ak2 = new G4double[numberOfVariables];
60  ak3 = new G4double[numberOfVariables];
61  ak4 = new G4double[numberOfVariables];
62  ak5 = new G4double[numberOfVariables];
63  ak6 = new G4double[numberOfVariables];
64  ak7 = new G4double[numberOfVariables];
65  ak8 = new G4double[numberOfVariables];
66  ak9 = new G4double[numberOfVariables];
67  ak10 = new G4double[numberOfVariables];
68  ak11 = new G4double[numberOfVariables];
69  ak12 = new G4double[numberOfVariables];
70  ak13 = new G4double[numberOfVariables];
71 
72  const G4int numStateVars = std::max(noIntegrationVariables, 8);
73  yTemp = new G4double[numStateVars];
74  yIn = new G4double[numStateVars] ;
75 
76  fLastInitialVector = new G4double[numStateVars] ;
77  fLastFinalVector = new G4double[numStateVars] ;
78 
79  fLastDyDx = new G4double[numStateVars];
80 
81  fMidVector = new G4double[numStateVars];
82  fMidError = new G4double[numStateVars];
83 
84  if( primary )
85  {
86  fAuxStepper = new G4DormandPrinceRK78(EqRhs, numberOfVariables,
87  !primary);
88  }
89 }
G4DormandPrinceRK78::G4DormandPrinceRK78
G4DormandPrinceRK78(G4EquationOfMotion *EqRhs, G4int numberOfVariables=6, G4bool primary=true)
Definition: G4DormandPrinceRK78.cc:49
G4int
int G4int
Definition: G4Types.hh:78
G4MagIntegratorStepper::G4MagIntegratorStepper
G4MagIntegratorStepper(G4EquationOfMotion *Equation, G4int numIntegrationVariables, G4int numStateVariables=12, bool isFSAL=false)
Definition: G4MagIntegratorStepper.cc:36
G4INCL::Math::max
T max(const T t1, const T t2)
brief Return the largest of the two arguments
Definition: G4INCLGlobals.hh:112
G4double
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

G4DormandPrinceRK78::~G4DormandPrinceRK78 ( )

Definition at line 92 of file G4DormandPrinceRK78.cc.

92  {
93  //clear all previously allocated memory for stepper and DistChord
94  delete[] ak2;
95  delete[] ak3;
96  delete[] ak4;
97  delete[] ak5;
98  delete[] ak6;
99  delete[] ak7;
100  delete[] ak8;
101  delete[] ak9;
102  delete[] ak10;
103  delete[] ak11;
104  delete[] ak12;
105  delete[] ak13;
106  delete[] yTemp;
107  delete[] yIn;
108 
109  delete[] fLastInitialVector;
110  delete[] fLastFinalVector;
111  delete[] fLastDyDx;
112  delete[] fMidVector;
113  delete[] fMidError;
114 
115  delete fAuxStepper;
116 
117 }

Member Function Documentation

G4double G4DormandPrinceRK78::DistChord ( ) const
virtual

Implements G4MagIntegratorStepper.

Definition at line 416 of file G4DormandPrinceRK78.cc.

417 {
418  G4double distLine, distChord;
419  G4ThreeVector initialPoint, finalPoint, midPoint;
420 
421  // Store last initial and final points (they will be overwritten in self-Stepper call!)
422  initialPoint = G4ThreeVector( fLastInitialVector[0],
423  fLastInitialVector[1], fLastInitialVector[2]);
424  finalPoint = G4ThreeVector( fLastFinalVector[0],
425  fLastFinalVector[1], fLastFinalVector[2]);
426 
427  // Do half a step using StepNoErr
428 
429  fAuxStepper->Stepper( fLastInitialVector, fLastDyDx, 0.5 * fLastStepLength,
430  fMidVector, fMidError );
431 
432  midPoint = G4ThreeVector( fMidVector[0], fMidVector[1], fMidVector[2]);
433 
434  // Use stored values of Initial and Endpoint + new Midpoint to evaluate
435  // distance of Chord
436 
437 
438  if (initialPoint != finalPoint)
439  {
440  distLine = G4LineSection::Distline( midPoint, initialPoint, finalPoint );
441  distChord = distLine;
442  }
443  else
444  {
445  distChord = (midPoint-initialPoint).mag();
446  }
447  return distChord;
448 }
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:42
G4LineSection::Distline
static G4double Distline(const G4ThreeVector &OtherPnt, const G4ThreeVector &LinePntA, const G4ThreeVector &LinePntB)
Definition: G4LineSection.hh:84
G4DormandPrinceRK78::Stepper
void Stepper(const G4double y[], const G4double dydx[], G4double h, G4double yout[], G4double yerr[])
Definition: G4DormandPrinceRK78.cc:132
G4double
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

G4int G4DormandPrinceRK78::IntegratorOrder ( ) const
inlinevirtual

Implements G4MagIntegratorStepper.

Definition at line 50 of file G4DormandPrinceRK78.hh.

50 {return 7; }
void G4DormandPrinceRK78::Stepper ( const G4double  y[],
const G4double  dydx[],
G4double  h,
G4double  yout[],
G4double  yerr[] 
)
virtual

Implements G4MagIntegratorStepper.

Definition at line 132 of file G4DormandPrinceRK78.cc.

137 {
138  G4int i;
139 
140  //The various constants defined on the basis of butcher tableu
141  //G4double - only once
142  const G4double
143 
144  b21 = 1.0/18,
145 
146  b31 = 1.0/48.0 ,
147  b32 = 1.0/16.0 ,
148 
149  b41 = 1.0/32.0 ,
150  b42 = 0.0 ,
151  b43 = 3.0/32.0 ,
152 
153  b51 = 5.0/16.0 ,
154  b52 = 0.0 ,
155  b53 = -75.0/64.0 ,
156  b54 = 75.0/64.0 ,
157 
158  b61 = 3.0/80.0 ,
159  b62 = 0.0 ,
160  b63 = 0.0 ,
161  b64 = 3.0/16.0 ,
162  b65 = 3.0/20.0 ,
163 
164  b71 = 29443841.0/614563906.0 ,
165  b72 = 0.0 ,
166  b73 = 0.0 ,
167  b74 = 77736538.0/692538347.0 ,
168  b75 = -28693883.0/1125000000.0 ,
169  b76 = 23124283.0/1800000000.0 ,
170 
171  b81 = 16016141.0/946692911.0 ,
172  b82 = 0.0 ,
173  b83 = 0.0 ,
174  b84 = 61564180.0/158732637.0 ,
175  b85 = 22789713.0/633445777.0 ,
176  b86 = 545815736.0/2771057229.0 ,
177  b87 = -180193667.0/1043307555.0 ,
178 
179  b91 = 39632708.0/573591083.0 ,
180  b92 = 0.0 ,
181  b93 = 0.0 ,
182  b94 = -433636366.0/683701615.0 ,
183  b95 = -421739975.0/2616292301.0 ,
184  b96 = 100302831.0/723423059.0 ,
185  b97 = 790204164.0/839813087.0 ,
186  b98 = 800635310.0/3783071287.0 ,
187 
188  b101 = 246121993.0/1340847787.0 ,
189  b102 = 0.0 ,
190  b103 = 0.0 ,
191  b104 = -37695042795.0/15268766246.0 ,
192  b105 = -309121744.0/1061227803.0 ,
193  b106 = -12992083.0/490766935.0 ,
194  b107 = 6005943493.0/2108947869.0 ,
195  b108 = 393006217.0/1396673457.0 ,
196  b109 = 123872331.0/1001029789.0 ,
197 
198  b111 = -1028468189.0/846180014.0 ,
199  b112 = 0.0 ,
200  b113 = 0.0 ,
201  b114 = 8478235783.0/508512852.0 ,
202  b115 = 1311729495.0/1432422823.0 ,
203  b116 = -10304129995.0/1701304382.0 ,
204  b117 = -48777925059.0/3047939560.0 ,
205  b118 = 15336726248.0/1032824649.0 ,
206  b119 = -45442868181.0/3398467696.0 ,
207  b1110 = 3065993473.0/597172653.0 ,
208 
209  b121 = 185892177.0/718116043.0 ,
210  b122 = 0.0 ,
211  b123 = 0.0 ,
212  b124 = -3185094517.0/667107341.0 ,
213  b125 = -477755414.0/1098053517.0 ,
214  b126 = -703635378.0/230739211.0 ,
215  b127 = 5731566787.0/1027545527.0 ,
216  b128 = 5232866602.0/850066563.0 ,
217  b129 = -4093664535.0/808688257.0 ,
218  b1210 = 3962137247.0/1805957418.0 ,
219  b1211 = 65686358.0/487910083.0 ,
220 
221  b131 = 403863854.0/491063109.0 ,
222  b132 = 0.0 ,
223  b133 = 0.0 ,
224  b134 = -5068492393.0/434740067.0 ,
225  b135 = -411421997.0/543043805.0 ,
226  b136 = 652783627.0/914296604.0 ,
227  b137 = 11173962825.0/925320556.0 ,
228  b138 = -13158990841.0/6184727034.0 ,
229  b139 = 3936647629.0/1978049680.0 ,
230  b1310 = -160528059.0/685178525.0 ,
231  b1311 = 248638103.0/1413531060.0 ,
232  b1312 = 0.0 ,
233 
234  c1 = 14005451.0/335480064.0 ,
235  // c2 = 0.0 ,
236  // c3 = 0.0 ,
237  // c4 = 0.0 ,
238  // c5 = 0.0 ,
239  c6 = -59238493.0/1068277825.0 ,
240  c7 = 181606767.0/758867731.0 ,
241  c8 = 561292985.0/797845732.0 ,
242  c9 = -1041891430.0/1371343529.0 ,
243  c10 = 760417239.0/1151165299.0 ,
244  c11 = 118820643.0/751138087.0 ,
245  c12 = - 528747749.0/2220607170.0 ,
246  c13 = 1.0/4.0 ,
247 
248  c_1 = 13451932.0/455176623.0 ,
249  // c_2 = 0.0 ,
250  // c_3 = 0.0 ,
251  // c_4 = 0.0 ,
252  // c_5 = 0.0 ,
253  c_6 = -808719846.0/976000145.0 ,
254  c_7 = 1757004468.0/5645159321.0 ,
255  c_8 = 656045339.0/265891186.0 ,
256  c_9 = -3867574721.0/1518517206.0 ,
257  c_10 = 465885868.0/322736535.0 ,
258  c_11 = 53011238.0/667516719.0 ,
259  c_12 = 2.0/45.0 ,
260  c_13 = 0.0 ,
261 
262  dc1 = c_1 - c1 ,
263  // dc2 = c_2 - c2 ,
264  // dc3 = c_3 - c3 ,
265  // dc4 = c_4 - c4 ,
266  // dc5 = c_5 - c5 ,
267  dc6 = c_6 - c6 ,
268  dc7 = c_7 - c7 ,
269  dc8 = c_8 - c8 ,
270  dc9 = c_9 - c9 ,
271  dc10 = c_10 - c10 ,
272  dc11 = c_11 - c11 ,
273  dc12 = c_12 - c12 ,
274  dc13 = c_13 - c13 ;
275 
276  //end of declaration !
277 
278  const G4int numberOfVariables= this->GetNumberOfVariables();
279 
280  // The number of variables to be integrated over
281  yOut[7] = yTemp[7] = yIn[7];
282  // Saving yInput because yInput and yOut can be aliases for same array
283 
284  for(i=0;i<numberOfVariables;i++)
285  {
286  yIn[i]=yInput[i];
287  }
288 
289  // RightHandSide(yIn, dydx) ;
290  // 1st Stage - Not doing, getting passed
291 
292  for(i=0;i<numberOfVariables;i++)
293  {
294  yTemp[i] = yIn[i] + b21*Step*dydx[i] ;
295  }
296  RightHandSide(yTemp, ak2) ; // 2nd Stage
297 
298  for(i=0;i<numberOfVariables;i++)
299  {
300  yTemp[i] = yIn[i] + Step*(b31*dydx[i] + b32*ak2[i]) ;
301  }
302  RightHandSide(yTemp, ak3) ; // 3rd Stage
303 
304  for(i=0;i<numberOfVariables;i++)
305  {
306  yTemp[i] = yIn[i] + Step*(b41*dydx[i] + b42*ak2[i] + b43*ak3[i]) ;
307  }
308  RightHandSide(yTemp, ak4) ; // 4th Stage
309 
310  for(i=0;i<numberOfVariables;i++)
311  {
312  yTemp[i] = yIn[i] + Step*(b51*dydx[i] + b52*ak2[i] + b53*ak3[i] +
313  b54*ak4[i]) ;
314  }
315  RightHandSide(yTemp, ak5) ; // 5th Stage
316 
317  for(i=0;i<numberOfVariables;i++)
318  {
319  yTemp[i] = yIn[i] + Step*(b61*dydx[i] + b62*ak2[i] + b63*ak3[i] +
320  b64*ak4[i] + b65*ak5[i]) ;
321  }
322  RightHandSide(yTemp, ak6) ; // 6th Stage
323 
324  for(i=0;i<numberOfVariables;i++)
325  {
326  yTemp[i] = yIn[i] + Step*(b71*dydx[i] + b72*ak2[i] + b73*ak3[i] +
327  b74*ak4[i] + b75*ak5[i] + b76*ak6[i]);
328  }
329  RightHandSide(yTemp, ak7); //7th Stage
330 
331  for(i=0;i<numberOfVariables;i++)
332  {
333  yTemp[i] = yIn[i] + Step*(b81*dydx[i] + b82*ak2[i] + b83*ak3[i] +
334  b84*ak4[i] + b85*ak5[i] + b86*ak6[i] +
335  b87*ak7[i]);
336  }
337  RightHandSide(yTemp, ak8); //8th Stage
338 
339  for(i=0;i<numberOfVariables;i++)
340  {
341  yTemp[i] = yIn[i] + Step*(b91*dydx[i] + b92*ak2[i] + b93*ak3[i] +
342  b94*ak4[i] + b95*ak5[i] + b96*ak6[i] +
343  b97*ak7[i] + b98*ak8[i] );
344  }
345  RightHandSide(yTemp, ak9); //9th Stage
346 
347  for(i=0;i<numberOfVariables;i++)
348  {
349  yTemp[i] = yIn[i] + Step*(b101*dydx[i] + b102*ak2[i] + b103*ak3[i] +
350  b104*ak4[i] + b105*ak5[i] + b106*ak6[i] +
351  b107*ak7[i] + b108*ak8[i] + b109*ak9[i]);
352  }
353  RightHandSide(yTemp, ak10); //10th Stage
354 
355  for(i=0;i<numberOfVariables;i++)
356  {
357  yTemp[i] = yIn[i] + Step*(b111*dydx[i] + b112*ak2[i] + b113*ak3[i] +
358  b114*ak4[i] + b115*ak5[i] + b116*ak6[i] +
359  b117*ak7[i] + b118*ak8[i] + b119*ak9[i] +
360  b1110*ak10[i]);
361  }
362  RightHandSide(yTemp, ak11); //11th Stage
363 
364  for(i=0;i<numberOfVariables;i++)
365  {
366  yTemp[i] = yIn[i] + Step*(b121*dydx[i] + b122*ak2[i] + b123*ak3[i] +
367  b124*ak4[i] + b125*ak5[i] + b126*ak6[i] +
368  b127*ak7[i] + b128*ak8[i] + b129*ak9[i] +
369  b1210*ak10[i] + b1211*ak11[i]);
370  }
371  RightHandSide(yTemp, ak12); //12th Stage
372 
373  for(i=0;i<numberOfVariables;i++)
374  {
375  yTemp[i] = yIn[i] + Step*(b131*dydx[i] + b132*ak2[i] + b133*ak3[i] +
376  b134*ak4[i] + b135*ak5[i] + b136*ak6[i] +
377  b137*ak7[i] + b138*ak8[i] + b139*ak9[i] +
378  b1310*ak10[i] + b1311*ak11[i] + b1312*ak12[i]);
379  }
380  RightHandSide(yTemp, ak13); //13th and final Stage
381 
382  for(i=0;i<numberOfVariables;i++)
383  {
384  // Accumulate increments with proper weights
385 
386  yOut[i] = yIn[i] + Step*(c1*dydx[i] + // c2 * ak2[i] + c3 * ak3[i]
387  // + c4 * ak4[i] + c5 * ak5[i]
388  + c6*ak6[i] +
389  c7*ak7[i] + c8*ak8[i] +c9*ak9[i] + c10*ak10[i]
390  + c11*ak11[i] + c12*ak12[i] + c13*ak13[i]) ;
391 
392  // Estimate error as difference between 7th and
393  // 8th order methods
394 
395  yErr[i] = Step*(dc1*dydx[i] + // dc2*ak2[i] + dc3*ak3[i] + dc4*ak4[i]
396  // + dc5*ak5[i]
397  + dc6*ak6[i] + dc7*ak7[i] + dc8*ak8[i]
398  + dc9*ak9[i] + dc10*ak10[i] + dc11*ak11[i] + dc12*ak12[i]
399  + dc13*ak13[i] ) ;
400  // Store Input and Final values, for possible use in calculating chord
401  fLastInitialVector[i] = yIn[i] ;
402  fLastFinalVector[i] = yOut[i];
403  fLastDyDx[i] = dydx[i];
404 
405 
406  }
407 
408  fLastStepLength = Step;
409 
410  return ;
411 }
G4int
int G4int
Definition: G4Types.hh:78
G4MagIntegratorStepper::GetNumberOfVariables
G4int GetNumberOfVariables() const
G4MagIntegratorStepper::RightHandSide
void RightHandSide(const double y[], double dydx[])
G4double
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

Here is the caller graph for this function:

The documentation for this class was generated from the following files: