Geant4 11.2.2
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
Public Member Functions | List of all members
G4FSALBogackiShampine45 Class Reference

#include <G4FSALBogackiShampine45.hh>

+ Inheritance diagram for G4FSALBogackiShampine45:

Public Member Functions

 G4FSALBogackiShampine45 (G4EquationOfMotion *EqRhs, G4int numberOfVariables=6, G4bool primary=true)
 
 ~G4FSALBogackiShampine45 () override
 
 G4FSALBogackiShampine45 (const G4FSALBogackiShampine45 &)=delete
 
G4FSALBogackiShampine45 & operator= (const G4FSALBogackiShampine45 &)=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)
 
G4double DistChord () const override
 
G4int IntegratorOrder () const override
 
- 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

Definition at line 40 of file G4FSALBogackiShampine45.hh.

Constructor & Destructor Documentation

◆ G4FSALBogackiShampine45() [1/2]

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

Definition at line 59 of file G4FSALBogackiShampine45.cc.

62 : G4VFSALIntegrationStepper(EqRhs, noIntegrationVariables)
63{
64 const G4int numberOfVariables = noIntegrationVariables;
65
66 // New Chunk of memory being created for use by the stepper
67
68 // aki - for storing intermediate RHS
69 //
70 ak2 = new G4double[numberOfVariables];
71 ak3 = new G4double[numberOfVariables];
72 ak4 = new G4double[numberOfVariables];
73 ak5 = new G4double[numberOfVariables];
74 ak6 = new G4double[numberOfVariables];
75 ak7 = new G4double[numberOfVariables];
76 ak8 = new G4double[numberOfVariables];
77
78 ak9 = new G4double[numberOfVariables];
79 ak10 = new G4double[numberOfVariables];
80 ak11 = new G4double[numberOfVariables];
81 DyDx = new G4double[numberOfVariables];
82
83 assert ( GetNumberOfStateVariables() >= 8 );
84 const G4int numStateVars = std::max(noIntegrationVariables,
85 GetNumberOfStateVariables() );
86
87 // Must ensure space extra 'state' variables exists - i.e. yIn[7]
88 //
89 yTemp = new G4double[numStateVars];
90 yIn = new G4double[numStateVars] ;
91
92 fLastInitialVector = new G4double[numStateVars] ;
93 fLastFinalVector = new G4double[numStateVars] ;
94 fLastDyDx = new G4double[numberOfVariables]; // Only derivatives
95
96 fMidVector = new G4double[numStateVars];
97 fMidError = new G4double[numStateVars];
98
99 pseudoDydx_for_DistChord = new G4double[numberOfVariables];
100
101 fMidVector = new G4double[numberOfVariables];
102 fMidError = new G4double[numberOfVariables];
103 if( primary )
104 {
105 fAuxStepper = new G4FSALBogackiShampine45(EqRhs, numberOfVariables,
106 !primary);
107 }
108 if( !fPreparedConstants )
109 {
110 PrepareConstants();
111 }
112}
G4double
double G4double
Definition G4Types.hh:83
G4int
int G4int
Definition G4Types.hh:85
G4FSALBogackiShampine45::G4FSALBogackiShampine45
G4FSALBogackiShampine45(G4EquationOfMotion *EqRhs, G4int numberOfVariables=6, G4bool primary=true)
Definition G4FSALBogackiShampine45.cc:59
G4VFSALIntegrationStepper::G4VFSALIntegrationStepper
G4VFSALIntegrationStepper(G4EquationOfMotion *Equation, G4int numIntegrationVariables, G4int numStateVariables=12)
Definition G4VFSALIntegrationStepper.cc:37
G4VFSALIntegrationStepper::GetNumberOfStateVariables
G4int GetNumberOfStateVariables() const

Referenced by G4FSALBogackiShampine45().

◆ ~G4FSALBogackiShampine45()

G4FSALBogackiShampine45::~G4FSALBogackiShampine45 ( )
override

Definition at line 116 of file G4FSALBogackiShampine45.cc.

117{
118 // Clear all previously allocated memory for stepper and DistChord
119
120 delete [] ak2;
121 delete [] ak3;
122 delete [] ak4;
123 delete [] ak5;
124 delete [] ak6;
125 delete [] ak7;
126 delete [] ak8;
127 delete [] ak9;
128 delete [] ak10;
129 delete [] ak11;
130 delete [] DyDx;
131 delete [] yTemp;
132 delete [] yIn;
133
134 delete [] fLastInitialVector;
135 delete [] fLastFinalVector;
136 delete [] fLastDyDx;
137 delete [] fMidVector;
138 delete [] fMidError;
139
140 delete fAuxStepper;
141
142 delete [] pseudoDydx_for_DistChord;
143}

◆ G4FSALBogackiShampine45() [2/2]

G4FSALBogackiShampine45::G4FSALBogackiShampine45 ( const G4FSALBogackiShampine45 & )
delete

Member Function Documentation

◆ DistChord()

G4double G4FSALBogackiShampine45::DistChord ( ) const
overridevirtual

Implements G4VFSALIntegrationStepper.

Definition at line 293 of file G4FSALBogackiShampine45.cc.

294{
295 G4double distLine, distChord;
296 G4ThreeVector initialPoint, finalPoint, midPoint;
297
298 // Store last initial and final points
299 // (they will be overwritten in self-Stepper call!)
300 //
301 initialPoint = G4ThreeVector( fLastInitialVector[0],
302 fLastInitialVector[1], fLastInitialVector[2]);
303 finalPoint = G4ThreeVector( fLastFinalVector[0],
304 fLastFinalVector[1], fLastFinalVector[2]);
305
306 // Do half a step using StepNoErr
307
308 fAuxStepper->Stepper( fLastInitialVector, fLastDyDx, 0.5 * fLastStepLength,
309 fMidVector, fMidError, pseudoDydx_for_DistChord );
310
311 midPoint = G4ThreeVector( fMidVector[0], fMidVector[1], fMidVector[2] );
312
313 // Use stored values of Initial and Endpoint + new Midpoint to evaluate
314 // distance of Chord
315 //
316 if (initialPoint != finalPoint)
317 {
318 distLine = G4LineSection::Distline(midPoint, initialPoint, finalPoint);
319 distChord = distLine;
320 }
321 else
322 {
323 distChord = (midPoint-initialPoint).mag();
324 }
325 return distChord;
326}
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition G4ThreeVector.hh:36
G4FSALBogackiShampine45::Stepper
void Stepper(const G4double y[], const G4double dydx[], G4double h, G4double yout[], G4double yerr[], G4double nextDydx[]) override
Definition G4FSALBogackiShampine45.cc:150
G4LineSection::Distline
static G4double Distline(const G4ThreeVector &OtherPnt, const G4ThreeVector &LinePntA, const G4ThreeVector &LinePntB)
Definition G4LineSection.hh:72

◆ IntegratorOrder()

G4int G4FSALBogackiShampine45::IntegratorOrder ( ) const
inlineoverridevirtual

Implements G4VFSALIntegrationStepper.

Definition at line 66 of file G4FSALBogackiShampine45.hh.

66{ return 4; }

◆ interpolate()

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

Definition at line 410 of file G4FSALBogackiShampine45.cc.

415{
416 const G4double a91 = 455.0/6144.0 ,
417 a92 = 0.0 ,
418 a93 = 10256301.0/35409920.0 ,
419 a94 = 2307361.0/17971200.0 ,
420 a95 = -387.0/102400.0 ,
421 a96 = 73.0/5130.0 ,
422 a97 = -7267.0/215040.0 ,
423 a98 = 1.0/32.0 ,
424
425 a101 = -837888343715.0/13176988637184.0 ,
426 a102 = 30409415.0/52955362.0 ,
427 a103 = -48321525963.0/759168069632.0 ,
428 a104 = 8530738453321.0/197654829557760.0 ,
429 a105 = 1361640523001.0/1626788720640.0 ,
430 a106 = -13143060689.0/38604458898.0 ,
431 a107 = 18700221969.0/379584034816.0 ,
432 a108 = -5831595.0/847285792.0 ,
433 a109 = -5183640.0/26477681.0 ,
434
435 a111 = 98719073263.0/1551965184000.0 ,
436 a112 = 1307.0/123552.0 ,
437 a113 = 4632066559387.0/70181753241600.0 ,
438 a114 = 7828594302389.0/382182512025600.0 ,
439 a115 = 40763687.0/11070259200.0 ,
440 a116 = 34872732407.0/224610586200.0 ,
441 a117 = -2561897.0/30105600.0 ,
442 a118 = 1.0/10.0 ,
443 a119 = -1.0/10.0 ,
444 a1110 = -1403317093.0/11371610250.0 ;
445
446 const G4int numberOfVariables = GetNumberOfVariables();
447
448 // Saving yInput because yInput and yOut can be aliases for same array
449 //
450 for(G4int i=0; i<numberOfVariables; ++i)
451 {
452 yIn[i]=yInput[i];
453 }
454
455 // The number of variables to be integrated over
456 //
457 yOut[7] = yTemp[7] = yIn[7];
458
459 // Calculating extra stages
460 //
461 for(G4int i=0; i<numberOfVariables; ++i)
462 {
463 yTemp[i] = yIn[i] + Step*(a91*dydx[i] + a92*ak2[i] + a93*ak3[i] +
464 a94*ak4[i] + a95*ak5[i] + a96*ak6[i] +
465 a97*ak7[i] + a98*ak8[i] );
466 }
467
468 RightHandSide(yTemp, ak9);
469
470 for(G4int i=0; i<numberOfVariables; ++i)
471 {
472 yTemp[i] = yIn[i] + Step*(a101*dydx[i] + a102*ak2[i] + a103*ak3[i] +
473 a104*ak4[i] + a105*ak5[i] + a106*ak6[i] +
474 a107*ak7[i] + a108*ak8[i] + a109*ak9[i] );
475 }
476
477 RightHandSide(yTemp, ak10);
478
479 for(G4int i=0; i<numberOfVariables; ++i)
480 {
481 yTemp[i] = yIn[i] + Step*(a111*dydx[i] + a112*ak2[i] + a113*ak3[i] +
482 a114*ak4[i] + a115*ak5[i] + a116*ak6[i] +
483 a117*ak7[i] + a118*ak8[i] + a119*ak9[i] +
484 a1110*ak10[i] );
485 }
486
487 RightHandSide(yTemp, ak11);
488
489 G4double tau0 = tau;
490
491 // Calculating the polynomials
492 //
493 for(auto i=1; i<=11; ++i) // i is NOT the coordinate no., it's stage no.
494 {
495 b[i] = 0.0;
496 tau = tau0;
497 for(auto j=1; j<=6; ++j)
498 {
499 b[i] += bi[i][j]*tau;
500 tau*=tau0;
501 }
502 }
503
504 for(G4int i=0; i<numberOfVariables; ++i)
505 {
506 yOut[i] = yIn[i] + Step*(b[1]*dydx[i] + b[2]*ak2[i] + b[3]*ak3[i] +
507 b[4]*ak4[i] + b[5]*ak5[i] + b[6]*ak6[i] +
508 b[7]*ak7[i] + b[8]*ak8[i] + b[9]*ak9[i] +
509 b[10]*ak10[i] + b[11]*ak11[i] );
510 }
511}
G4VFSALIntegrationStepper::GetNumberOfVariables
G4int GetNumberOfVariables() const
G4VFSALIntegrationStepper::RightHandSide
void RightHandSide(const double y[], double dydx[])
Definition G4VFSALIntegrationStepper.cc:51

◆ operator=()

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

◆ Stepper()

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

Implements G4VFSALIntegrationStepper.

Definition at line 150 of file G4FSALBogackiShampine45.cc.

156{
157 G4int i;
158
159 // The various constants defined on the basis of butcher tableu
160
161 const G4double b21 = 1.0/6.0 ,
162 b31 = 2.0/27.0 , b32 = 4.0/27.0,
163
164 b41 = 183.0/1372.0 , b42 = -162.0/343.0, b43 = 1053.0/1372.0,
165
166 b51 = 68.0/297.0, b52 = -4.0/11.0,
167 b53 = 42.0/143.0, b54 = 1960.0/3861.0,
168
169 b61 = 597.0/22528.0, b62 = 81.0/352.0,
170 b63 = 63099.0/585728.0, b64 = 58653.0/366080.0,
171 b65 = 4617.0/20480.0,
172
173 b71 = 174197.0/959244.0, b72 = -30942.0/79937.0,
174 b73 = 8152137.0/19744439.0, b74 = 666106.0/1039181.0,
175 b75 = -29421.0/29068.0, b76 = 482048.0/414219.0,
176
177 b81 = 587.0/8064.0, b82 = 0.0,
178 b83 = 4440339.0/15491840.0, b84 = 24353.0/124800.0,
179 b85 = 387.0/44800.0, b86 = 2152.0/5985.0,
180 b87 = 7267.0/94080.0,
181
182
183 // c1 = 2479.0/34992.0,
184 // c2 = 0.0,
185 // c3 = 123.0/416.0,
186 // c4 = 612941.0/3411720.0,
187 // c5 = 43.0/1440.0,
188 // c6 = 2272.0/6561.0,
189 // c7 = 79937.0/1113912.0,
190 // c8 = 3293.0/556956.0,
191
192 // For the embedded higher order method only the difference of values
193 // taken and is used directly later instead of defining the last row
194 // of butcher table in a separate set of variables and taking the
195 // difference there
196
197 dc1 = b81 - 2479.0/34992.0 ,
198 dc2 = 0.0,
199 dc3 = b83 - 123.0/416.0 ,
200 dc4 = b84 - 612941.0/3411720.0,
201 dc5 = b85 - 43.0/1440.0,
202 dc6 = b86 - 2272.0/6561.0,
203 dc7 = b87 - 79937.0/1113912.0,
204 dc8 = -3293.0/556956.0; // end of declaration
205
206 const G4int numberOfVariables = GetNumberOfVariables();
207
208 // The number of variables to be integrated over
209 //
210 yOut[7] = yTemp[7] = yIn[7];
211
212 // Saving yInput because yInput and yOut can be aliases for same array
213 //
214 for(i=0; i<numberOfVariables; ++i)
215 {
216 yIn[i]=yInput[i];
217 DyDx[i] = dydx[i];
218 }
219 // RightHandSide(yIn, dydx) ; // 1st Step - Not doing, getting passed
220
221 for(i=0; i<numberOfVariables; ++i)
222 {
223 yTemp[i] = yIn[i] + b21*Step*DyDx[i] ;
224 }
225 RightHandSide(yTemp, ak2) ; // 2nd Step
226
227 for(i=0; i<numberOfVariables; ++i)
228 {
229 yTemp[i] = yIn[i] + Step*(b31*DyDx[i] + b32*ak2[i]) ;
230 }
231 RightHandSide(yTemp, ak3) ; // 3rd Step
232
233 for(i=0; i<numberOfVariables; ++i)
234 {
235 yTemp[i] = yIn[i] + Step*(b41*DyDx[i] + b42*ak2[i] + b43*ak3[i]) ;
236 }
237 RightHandSide(yTemp, ak4) ; // 4th Step
238
239 for(i=0; i<numberOfVariables; ++i)
240 {
241 yTemp[i] = yIn[i] + Step*(b51*DyDx[i] + b52*ak2[i] + b53*ak3[i] +
242 b54*ak4[i]) ;
243 }
244 RightHandSide(yTemp, ak5) ; // 5th Step
245
246 for(i=0; i<numberOfVariables; ++i)
247 {
248 yTemp[i] = yIn[i] + Step*(b61*DyDx[i] + b62*ak2[i] + b63*ak3[i] +
249 b64*ak4[i] + b65*ak5[i]) ;
250 }
251 RightHandSide(yTemp, ak6) ; // 6th Step
252
253 for(i=0; i<numberOfVariables; ++i)
254 {
255 yTemp[i] = yIn[i] + Step*(b71*DyDx[i] + b72*ak2[i] + b73*ak3[i] +
256 b74*ak4[i] + b75*ak5[i] + b76*ak6[i]);
257 }
258 RightHandSide(yTemp, ak7); // 7th Step
259
260 for(i=0; i<numberOfVariables; ++i)
261 {
262 yOut[i] = yIn[i] + Step*(b81*DyDx[i] + b82*ak2[i] + b83*ak3[i] +
263 b84*ak4[i] + b85*ak5[i] + b86*ak6[i] +
264 b87*ak7[i]);
265 }
266 RightHandSide(yOut, ak8); // 8th Step - Final one Using FSAL
267
268
269 for(i=0; i<numberOfVariables; ++i)
270 {
271
272 yErr[i] = Step*(dc1*DyDx[i] + dc2*ak2[i] + dc3*ak3[i] + dc4*ak4[i] +
273 dc5*ak5[i] + dc6*ak6[i] + dc7*ak7[i] + dc8*ak8[i]) ;
274
275
276 // FSAL stepper : Must pass the last DyDx for the next step, here ak8
277 //
278 nextDydx[i] = ak8[i];
279
280 // Store Input and Final values, for possible use in calculating chord
281 //
282 fLastInitialVector[i] = yIn[i] ;
283 fLastFinalVector[i] = yOut[i];
284 fLastDyDx[i] = DyDx[i];
285 }
286 fLastStepLength = Step;
287
288 return;
289}

Referenced by DistChord().

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