Geant4 11.4.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4FSALDormandPrince745 Class Reference

G4FSALDormandPrince745 is an integrator of particle's equation of motion based on the DormandPrince7 - 5(4) FSAL implementation. More...

#include <G4FSALDormandPrince745.hh>

Inheritance diagram for G4FSALDormandPrince745:

Public Member Functions

 G4FSALDormandPrince745 (G4EquationOfMotion *EqRhs, G4int numberOfVariables=6, G4bool primary=true)
 ~G4FSALDormandPrince745 () override
 G4FSALDormandPrince745 (const G4FSALDormandPrince745 &)=delete
G4FSALDormandPrince745 & operator= (const G4FSALDormandPrince745 &)=delete
void Stepper (const G4double y[], const G4double dydx[], G4double h, G4double yout[], G4double yerr[], G4double nextDydx[]) override
void interpolate (const G4double yInput[], const G4double dydx[], G4double yOut[], G4double Step, G4double tau)
void Interpolate (const G4double yInput[], const G4double dydx[], const G4double Step, G4double yOut[], G4double tau)
void SetupInterpolate (const G4double yInput[], const G4double dydx[], const G4double Step)
G4double DistChord () const override
G4int IntegratorOrder () const override
G4bool isFSAL () const
Public Member Functions inherited from G4VFSALIntegrationStepper
 G4VFSALIntegrationStepper (G4EquationOfMotion *Equation, G4int numIntegrationVariables, G4int numStateVariables=12)
virtual ~G4VFSALIntegrationStepper ()=default
 G4VFSALIntegrationStepper (const G4VFSALIntegrationStepper &)=delete
G4VFSALIntegrationStepper & operator= (const G4VFSALIntegrationStepper &)=delete
void NormaliseTangentVector (G4double vec[6])
void NormalisePolarizationVector (G4double vec[12])
void RightHandSide (const double y[], double dydx[])
G4int GetNumberOfVariables () const
G4int GetNumberOfStateVariables () const
G4EquationOfMotion * GetEquationOfMotion ()
void SetEquationOfMotion (G4EquationOfMotion *newEquation)
G4int GetfNoRHSCalls ()
void increasefNORHSCalls ()
void ResetfNORHSCalls ()

Detailed Description

G4FSALDormandPrince745 is an integrator of particle's equation of motion based on the DormandPrince7 - 5(4) FSAL implementation.

Definition at line 45 of file G4FSALDormandPrince745.hh.

Constructor & Destructor Documentation

◆ G4FSALDormandPrince745() [1/2]

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

Constructor for G4FSALDormandPrince745.

Parameters
[in]EqRhsPointer to the provided equation of motion.
[in]numberOfVariablesThe number of integration variables.
[in]primaryFlag for initialisation of the auxiliary stepper.

Definition at line 51 of file G4FSALDormandPrince745.cc.

54 : G4VFSALIntegrationStepper(EqRhs, noIntegrationVariables)
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 //
62 ak2 = new G4double[numberOfVariables];
63 ak3 = new G4double[numberOfVariables];
64 ak4 = new G4double[numberOfVariables];
65 ak5 = new G4double[numberOfVariables];
66 ak6 = new G4double[numberOfVariables];
67 ak7 = new G4double[numberOfVariables];
68
69 // Also always allocate arrays for interpolation stages
70 //
71 ak8 = new G4double[numberOfVariables];
72 ak9 = new G4double[numberOfVariables];
73
74 yTemp = new G4double[numberOfVariables] ;
75 yIn = new G4double[numberOfVariables] ;
76
77 pseudoDydx_for_DistChord = new G4double[numberOfVariables];
78
79 fInitialDyDx = new G4double[numberOfVariables];
80 fLastInitialVector = new G4double[numberOfVariables] ;
81 fLastFinalVector = new G4double[numberOfVariables] ;
82 fLastDyDx = new G4double[numberOfVariables];
83
84 fMidVector = new G4double[numberOfVariables];
85 fMidError = new G4double[numberOfVariables];
86
87 if( primary )
88 {
89 fAuxStepper = new G4FSALDormandPrince745(EqRhs,numberOfVariables,!primary);
90 }
91}
G4double
double G4double
Definition G4Types.hh:83
G4int
int G4int
Definition G4Types.hh:85
G4FSALDormandPrince745::G4FSALDormandPrince745
G4FSALDormandPrince745(G4EquationOfMotion *EqRhs, G4int numberOfVariables=6, G4bool primary=true)
Definition G4FSALDormandPrince745.cc:51
G4VFSALIntegrationStepper::G4VFSALIntegrationStepper
G4VFSALIntegrationStepper(G4EquationOfMotion *Equation, G4int numIntegrationVariables, G4int numStateVariables=12)
Definition G4VFSALIntegrationStepper.cc:37

Referenced by G4FSALDormandPrince745(), G4FSALDormandPrince745(), and operator=().

◆ ~G4FSALDormandPrince745()

G4FSALDormandPrince745::~G4FSALDormandPrince745 ( )
override

Destructor.

Definition at line 95 of file G4FSALDormandPrince745.cc.

96{
97 // Clear all previously allocated memory for stepper and DistChord
98
99 delete [] ak2; ak2 = nullptr;
100 delete [] ak3; ak3 = nullptr;
101 delete [] ak4; ak4 = nullptr;
102 delete [] ak5; ak5 = nullptr;
103 delete [] ak6; ak6 = nullptr;
104 delete [] ak7; ak7 = nullptr;
105 delete [] ak8; ak8 = nullptr;
106 delete [] ak9; ak9 = nullptr;
107
108 delete [] yTemp; yTemp = nullptr;
109 delete [] yIn; yIn = nullptr;
110
111 delete [] pseudoDydx_for_DistChord; pseudoDydx_for_DistChord = nullptr;
112 delete [] fInitialDyDx; fInitialDyDx = nullptr;
113
114 delete [] fLastInitialVector; fLastInitialVector = nullptr;
115 delete [] fLastFinalVector; fLastFinalVector = nullptr;
116 delete [] fLastDyDx; fLastDyDx = nullptr;
117 delete [] fMidVector; fMidVector = nullptr;
118 delete [] fMidError; fMidError = nullptr;
119
120 delete fAuxStepper; fAuxStepper = nullptr;
121}

◆ G4FSALDormandPrince745() [2/2]

G4FSALDormandPrince745::G4FSALDormandPrince745 ( const G4FSALDormandPrince745 & )
delete

Copy constructor and assignment operator not allowed.

Member Function Documentation

◆ DistChord()

G4double G4FSALDormandPrince745::DistChord ( ) const
overridevirtual

Returns the distance from chord line.

Implements G4VFSALIntegrationStepper.

Definition at line 243 of file G4FSALDormandPrince745.cc.

244{
245 G4double distLine, distChord;
246 G4ThreeVector initialPoint, finalPoint, midPoint;
247
248 // Store last initial and final points
249 // (they will be overwritten in self-Stepper call!)
250 //
251 initialPoint = G4ThreeVector(fLastInitialVector[0],
252 fLastInitialVector[1], fLastInitialVector[2]);
253 finalPoint = G4ThreeVector(fLastFinalVector[0],
254 fLastFinalVector[1], fLastFinalVector[2]);
255
256 // Do half a step using StepNoErr
257
258 fAuxStepper->Stepper( fLastInitialVector, fLastDyDx, 0.5 * fLastStepLength,
259 fMidVector, fMidError, pseudoDydx_for_DistChord );
260
261 midPoint = G4ThreeVector( fMidVector[0], fMidVector[1], fMidVector[2] );
262
263 // Use stored values of Initial and Endpoint + new Midpoint to evaluate
264 // distance of Chord
265 //
266 if (initialPoint != finalPoint)
267 {
268 distLine = G4LineSection::Distline( midPoint,initialPoint,finalPoint );
269 distChord = distLine;
270 }
271 else
272 {
273 distChord = (midPoint-initialPoint).mag();
274 }
275 return distChord;
276}
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition G4ThreeVector.hh:36
G4LineSection::Distline
static G4double Distline(const G4ThreeVector &OtherPnt, const G4ThreeVector &LinePntA, const G4ThreeVector &LinePntB)
Definition G4LineSection.hh:96

◆ IntegratorOrder()

G4int G4FSALDormandPrince745::IntegratorOrder ( ) const
inlineoverridevirtual

Returns the order, 4, of integration.

Implements G4VFSALIntegrationStepper.

Definition at line 131 of file G4FSALDormandPrince745.hh.

131{ return 4; }

◆ Interpolate()

void G4FSALDormandPrince745::Interpolate ( const G4double yInput[],
const G4double dydx[],
const G4double Step,
G4double yOut[],
G4double tau )

Calculates the output at the tau fraction of step. Same as above for higher order interpolant.

Definition at line 384 of file G4FSALDormandPrince745.cc.

389{
390 // Define the coefficients for the polynomials
391
392 G4double bi[10][5], b[10];
393 G4int numberOfVariables = GetNumberOfVariables();
394
395 // COEFFICIENTS OF bi[1]
396 bi[1][0] = 1.0 ,
397 bi[1][1] = -38039.0/7040.0 ,
398 bi[1][2] = 125923.0/10560.0 ,
399 bi[1][3] = -19683.0/1760.0 ,
400 bi[1][4] = 3303.0/880.0 ,
401 // --------------------------------------------------------
402 //
403 // COEFFICIENTS OF bi[2]
404 bi[2][0] = 0.0 ,
405 bi[2][1] = 0.0 ,
406 bi[2][2] = 0.0 ,
407 bi[2][3] = 0.0 ,
408 bi[2][4] = 0.0 ,
409 // --------------------------------------------------------
410 //
411 // COEFFICIENTS OF bi[3]
412 bi[3][0] = 0.0 ,
413 bi[3][1] = -12500.0/4081.0 ,
414 bi[3][2] = 205000.0/12243.0 ,
415 bi[3][3] = -90000.0/4081.0 ,
416 bi[3][4] = 36000.0/4081.0 ,
417 // --------------------------------------------------------
418 //
419 // COEFFICIENTS OF bi[4]
420 bi[4][0] = 0.0 ,
421 bi[4][1] = -3125.0/704.0 ,
422 bi[4][2] = 25625.0/1056.0 ,
423 bi[4][3] = -5625.0/176.0 ,
424 bi[4][4] = 1125.0/88.0 ,
425 // --------------------------------------------------------
426 //
427 // COEFFICIENTS OF bi[5]
428 bi[5][0] = 0.0 ,
429 bi[5][1] = 164025.0/74624.0 ,
430 bi[5][2] = -448335.0/37312.0 ,
431 bi[5][3] = 295245.0/18656.0 ,
432 bi[5][4] = -59049.0/9328.0 ,
433 // --------------------------------------------------------
434 //
435 // COEFFICIENTS OF bi[6]
436 bi[6][0] = 0.0 ,
437 bi[6][1] = -25.0/28.0 ,
438 bi[6][2] = 205.0/42.0 ,
439 bi[6][3] = -45.0/7.0 ,
440 bi[6][4] = 18.0/7.0 ,
441 // --------------------------------------------------------
442 //
443 // COEFFICIENTS OF bi[7]
444 bi[7][0] = 0.0 ,
445 bi[7][1] = -2.0/11.0 ,
446 bi[7][2] = 73.0/55.0 ,
447 bi[7][3] = -171.0/55.0 ,
448 bi[7][4] = 108.0/55.0 ,
449 // --------------------------------------------------------
450 //
451 // COEFFICIENTS OF bi[8]
452 bi[8][0] = 0.0 ,
453 bi[8][1] = 189.0/22.0 ,
454 bi[8][2] = -1593.0/55.0 ,
455 bi[8][3] = 3537.0/110.0 ,
456 bi[8][4] = -648.0/55.0 ,
457 // --------------------------------------------------------
458 //
459 // COEFFICIENTS OF bi[9]
460 bi[9][0] = 0.0 ,
461 bi[9][1] = 351.0/110.0 ,
462 bi[9][2] = -999.0/55.0 ,
463 bi[9][3] = 2943.0/110.0 ,
464 bi[9][4] = -648.0/55.0 ;
465 // --------------------------------------------------------
466
467 for(G4int i = 0; i< numberOfVariables; ++i)
468 {
469 yIn[i] = yInput[i];
470 }
471
472 G4double tau0 = tau;
473
474 // Calculating the polynomials
475 //
476 for(auto i=1; i<=9; ++i) // i is NOT the coordinate no., it's stage no.
477 {
478 b[i] = 0;
479 tau = 1.0;
480 for(auto j=0; j<=4; ++j)
481 {
482 b[i] += bi[i][j]*tau;
483 tau*=tau0;
484 }
485 }
486
487 for(G4int i=0; i<numberOfVariables; ++i) // Here i IS the coordinate no.
488 {
489 yOut[i] = yIn[i] + Step*tau0*(b[1]*dydx[i] + b[2]*ak2[i] + b[3]*ak3[i] +
490 b[4]*ak4[i] + b[5]*ak5[i] + b[6]*ak6[i] +
491 b[7]*ak7[i] + b[8]*ak8[i] + b[9]*ak9[i] );
492 }
493}
G4VFSALIntegrationStepper::GetNumberOfVariables
G4int GetNumberOfVariables() const

◆ interpolate()

void G4FSALDormandPrince745::interpolate ( const G4double yInput[],
const G4double dydx[],
G4double yOut[],
G4double Step,
G4double tau )

Calculates the output at the tau fraction of step.

Parameters
[in]yInputStarting values array of integration variables.
[in]dydxDerivatives array.
[out]yOutInterpolation output.
[in]StepThe given step size.
[in]tauThe tau fraction of the step.

Definition at line 280 of file G4FSALDormandPrince745.cc.

285{
286 G4double bf1, bf2, bf3, bf4, bf5, bf6, bf7;
287
288 const G4int numberOfVariables = GetNumberOfVariables();
289
290 G4double tau0 = tau;
291
292 for(G4int i=0;i<numberOfVariables; ++i)
293 {
294 yIn[i]=yInput[i];
295 }
296
297 G4double tau_2 = tau0*tau0 ,
298 tau_3 = tau0*tau_2,
299 tau_4 = tau_2*tau_2;
300
301 bf1 = (157015080.0*tau_4 - 13107642775.0*tau_3
302 + 34969693132.0*tau_2- 32272833064.0*tau + 11282082432.0)
303 / 11282082432.0;
304 bf2 = 0.0;
305 bf3 = - 100.0*tau*(15701508.0*tau_3 - 914128567.0*tau_2
306 + 2074956840.0*tau - 1323431896.0) / 32700410799.0;
307 bf4 = 25.0*tau*(94209048.0*tau_3- 1518414297.0*tau_2
308 + 2460397220.0*tau - 889289856.0)
309 / 5641041216.0;
310 bf5 = -2187.0*tau*(52338360.0*tau_3 - 451824525.0*tau_2
311 + 687873124.0*tau - 259006536.0)
312 / 199316789632.0;
313 bf6 = 11.0*tau*(106151040.0*tau_3- 661884105.0*tau_2
314 + 946554244.0*tau - 361440756.0)
315 / 2467955532.0;
316 bf7 = tau*(1.0 - tau)*(8293050.0*tau_2 - 82437520.0*tau + 44764047.0)
317 / 29380423.0;
318
319 for(G4int i=0; i<numberOfVariables; ++i)
320 {
321 yOut[i] = yIn[i] + Step*tau*(bf1*dydx[i] + bf2*ak2[i] + bf3*ak3[i]
322 + bf4*ak4[i] + bf5*ak5[i] + bf6*ak6[i]
323 + bf7*ak7[i] );
324 }
325}

◆ isFSAL()

G4bool G4FSALDormandPrince745::isFSAL ( ) const
inline

Returns true as this is a FSAL integrator.

Definition at line 136 of file G4FSALDormandPrince745.hh.

136{ return true; }

◆ operator=()

G4FSALDormandPrince745 & G4FSALDormandPrince745::operator= ( const G4FSALDormandPrince745 & )
delete

◆ SetupInterpolate()

void G4FSALDormandPrince745::SetupInterpolate ( const G4double yInput[],
const G4double dydx[],
const G4double Step )

Setup method for higher order interpolant.

Parameters
[in]yInputStarting values array of integration variables.
[in]dydxDerivatives array.
[in]StepThe given step size.

Definition at line 329 of file G4FSALDormandPrince745.cc.

332{
333 // Coefficients for the additional stages
334 //
335 G4double b81 = 6245.0/62208.0 ,
336 b82 = 0.0 ,
337 b83 = 8875.0/103032.0 ,
338 b84 = -125.0/1728.0 ,
339 b85 = 801.0/13568.0 ,
340 b86 = -13519.0/368064.0 ,
341 b87 = 11105.0/368064.0 ,
342
343 b91 = 632855.0/4478976.0 ,
344 b92 = 0.0 ,
345 b93 = 4146875.0/6491016.0 ,
346 b94 = 5490625.0/14183424.0 ,
347 b95 = -15975.0/108544.0 ,
348 b96 = 8295925.0/220286304.0 ,
349 b97 = -1779595.0/62938944.0 ,
350 b98 = -805.0/4104.0 ;
351
352 const G4int numberOfVariables = GetNumberOfVariables();
353
354 // Saving yInput because yInput and yOut can be aliases for same array
355 //
356 for(G4int i=0; i<numberOfVariables; ++i)
357 {
358 yIn[i] = yInput[i];
359 }
360
361 yTemp[7] = yIn[7];
362
363 // Evaluate the extra stages
364 //
365 for(G4int i=0; i<numberOfVariables; ++i)
366 {
367 yTemp[i] = yIn[i] + Step*( b81*dydx[i] + b82*ak2[i] + b83*ak3[i] +
368 b84*ak4[i] + b85*ak5[i] + b86*ak6[i] +
369 b87*ak7[i] );
370 }
371 RightHandSide( yTemp, ak8 ); // 8th Stage
372
373 for(G4int i=0; i<numberOfVariables; ++i)
374 {
375 yTemp[i] = yIn[i] + Step * ( b91*dydx[i] + b92*ak2[i] + b93*ak3[i] +
376 b94*ak4[i] + b95*ak5[i] + b96*ak6[i] +
377 b97*ak7[i] + b98*ak8[i] );
378 }
379 RightHandSide( yTemp, ak9 ); // 9th Stage
380}
G4VFSALIntegrationStepper::RightHandSide
void RightHandSide(const double y[], double dydx[])
Definition G4VFSALIntegrationStepper.cc:51

◆ Stepper()

void G4FSALDormandPrince745::Stepper ( const G4double y[],
const G4double dydx[],
G4double h,
G4double yout[],
G4double yerr[],
G4double nextDydx[] )
overridevirtual

The stepper for the Runge Kutta integration. The stepsize is fixed, with the step size given by 'h'. Integrates ODE starting values y[0 to 6]. Outputs yout[] and its estimated error yerr[].

Parameters
[in]yStarting values array of integration variables.
[in]dydxDerivatives array.
[in]hThe given step size.
[out]youtIntegration output.
[out]yerrThe estimated error.
[out]nextDydxLast derivatives array for the next step.

Implements G4VFSALIntegrationStepper.

Definition at line 128 of file G4FSALDormandPrince745.cc.

134{
135 G4int i;
136
137 // The various constants defined on the basis of butcher tableu
138
139 const G4double b21 = 0.2 ,
140 b31 = 3.0/40.0, b32 = 9.0/40.0 ,
141
142 b41 = 44.0/45.0, b42 = -56.0/15.0, b43 = 32.0/9.0,
143
144 b51 = 19372.0/6561.0, b52 = -25360.0/2187.0,
145 b53 = 64448.0/6561.0, b54 = -212.0/729.0 ,
146
147 b61 = 9017.0/3168.0 , b62 = -355.0/33.0,
148 b63 = 46732.0/5247.0 , b64 = 49.0/176.0 ,
149 b65 = -5103.0/18656.0 ,
150
151 b71 = 35.0/384.0, b72 = 0.,
152 b73 = 500.0/1113.0, b74 = 125.0/192.0,
153 b75 = -2187.0/6784.0, b76 = 11.0/84.0,
154
155 // c1 = 35.0/384.0, c2 = .0,
156 // c3 = 500.0/1113.0, c4 = 125.0/192.0,
157 // c5 = -2187.0/6784.0, c6 = 11.0/84.0,
158 // c7 = 0,
159
160 dc1 = b71 - 5179.0/57600.0,
161 dc2 = b72 - .0,
162 dc3 = b73 - 7571.0/16695.0,
163 dc4 = b74 - 393.0/640.0,
164 dc5 = b75 + 92097.0/339200.0,
165 dc6 = b76 - 187.0/2100.0,
166 dc7 = - 1.0/40.0 ; //end of declaration
167
168 const G4int numberOfVariables = GetNumberOfVariables();
169 // The number of variables to be integrated over
170
171 // Saving yInput because yInput and yOut can be aliases for same array
172 //
173 for(i=0; i<numberOfVariables; ++i)
174 {
175 yIn[i] = yInput[i];
176 fInitialDyDx[i] = dydx[i];
177 }
178 // Ensure that time is initialised - in case it is not integrated
179 //
180 yOut[7] = yTemp[7] = yInput[7];
181 // RightHandSide(yIn, DyDx) ; // 1st Step - Not doing, getting passed
182
183 for(i=0; i<numberOfVariables; ++i)
184 {
185 yTemp[i] = yIn[i] + b21*Step*fInitialDyDx[i] ;
186 }
187 RightHandSide(yTemp, ak2) ; // 2nd Step
188
189 for(i=0; i<numberOfVariables; ++i)
190 {
191 yTemp[i] = yIn[i] + Step*(b31*fInitialDyDx[i] + b32*ak2[i]) ;
192 }
193 RightHandSide(yTemp, ak3) ; // 3rd Step
194
195 for(i=0; i<numberOfVariables; ++i)
196 {
197 yTemp[i] = yIn[i] + Step*(b41*fInitialDyDx[i]
198 + b42*ak2[i] + b43*ak3[i]) ;
199 }
200 RightHandSide(yTemp, ak4) ; // 4th Step
201
202 for(i=0; i<numberOfVariables; ++i)
203 {
204 yTemp[i] = yIn[i] + Step*(b51*fInitialDyDx[i]
205 + b52*ak2[i] + b53*ak3[i] + b54*ak4[i]) ;
206 }
207 RightHandSide(yTemp, ak5) ; // 5th Step
208
209 for(i=0; i<numberOfVariables; ++i)
210 {
211 yTemp[i] = yIn[i] + Step*(b61*fInitialDyDx[i] + b62*ak2[i]
212 + b63*ak3[i] + b64*ak4[i] + b65*ak5[i]) ;
213 }
214 RightHandSide(yTemp, ak6) ; // 6th Step
215
216 for(i=0; i<numberOfVariables; ++i)
217 {
218 yOut[i] = yIn[i] + Step*(b71*fInitialDyDx[i] + b72*ak2[i] + b73*ak3[i]
219 + b74*ak4[i] + b75*ak5[i] + b76*ak6[i]);
220 }
221 RightHandSide(yOut, ak7); //7th and Final step
222
223 for(i=0; i<numberOfVariables; ++i)
224 {
225
226 yErr[i] = Step*(dc1*fInitialDyDx[i] + dc2*ak2[i] + dc3*ak3[i]
227 + dc4*ak4[i] + dc5*ak5[i] + dc6*ak6[i] + dc7*ak7[i] ) ;
228
229 // Store Input and Final values, for possible use in calculating chord
230 //
231 fLastInitialVector[i] = yIn[i] ;
232 fLastFinalVector[i] = yOut[i];
233 fLastDyDx[i] = fInitialDyDx[i];
234 nextDydx[i] = ak7[i];
235 }
236 fLastStepLength = Step;
237
238 return ;
239}
The documentation for this class was generated from the following files: