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

G4TsitourasRK45 is an implementation of the 5(4) Runge-Kutta stepper (non-FSAL version). More...

#include <G4TsitourasRK45.hh>

Inheritance diagram for G4TsitourasRK45:

Public Member Functions

 G4TsitourasRK45 (G4EquationOfMotion *EqRhs, G4int numberOfVariables=6, G4bool primary=true)
 ~G4TsitourasRK45 () override
 G4TsitourasRK45 (const G4TsitourasRK45 &)=delete
G4TsitourasRK45 & operator= (const G4TsitourasRK45 &)=delete
void Stepper (const G4double y[], const G4double dydx[], G4double h, G4double yout[], G4double yerr[]) override
void SetupInterpolation ()
void Interpolate (const G4double yInput[], const G4double dydx[], const G4double Step, G4double yOut[], G4double tau)
void interpolate (const G4double yInput[], const G4double dydx[], G4double yOut[], G4double Step, G4double tau)
G4double DistChord () const override
G4int IntegratorOrder () const override
G4StepperType StepperType () const override
Public Member Functions inherited from G4MagIntegratorStepper
 G4MagIntegratorStepper (G4EquationOfMotion *Equation, G4int numIntegrationVariables, G4int numStateVariables=12, G4bool isFSAL=false)
virtual ~G4MagIntegratorStepper ()=default
 G4MagIntegratorStepper (const G4MagIntegratorStepper &)=delete
G4MagIntegratorStepper & operator= (const G4MagIntegratorStepper &)=delete
void NormaliseTangentVector (G4double vec[6])
void NormalisePolarizationVector (G4double vec[12])
void RightHandSide (const G4double y[], G4double dydx[]) const
void RightHandSide (const G4double y[], G4double dydx[], G4double field[]) const
G4int GetNumberOfVariables () const
G4int GetNumberOfStateVariables () const
G4int IntegrationOrder ()
G4EquationOfMotion * GetEquationOfMotion ()
const G4EquationOfMotion * GetEquationOfMotion () const
void SetEquationOfMotion (G4EquationOfMotion *newEquation)
unsigned long GetfNoRHSCalls ()
void ResetfNORHSCalls ()
G4bool IsFSAL () const
G4bool isQSS () const
void SetIsQSS (G4bool val)

Additional Inherited Members

Protected Member Functions inherited from G4MagIntegratorStepper
void SetIntegrationOrder (G4int order)
void SetFSAL (G4bool flag=true)

Detailed Description

G4TsitourasRK45 is an implementation of the 5(4) Runge-Kutta stepper (non-FSAL version).

Definition at line 50 of file G4TsitourasRK45.hh.

Constructor & Destructor Documentation

◆ G4TsitourasRK45() [1/2]

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

Constructor for G4TsitourasRK45.

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 39 of file G4TsitourasRK45.cc.

42 : G4MagIntegratorStepper(EqRhs, noIntegrationVariables)
43{
44 const G4int numberOfVariables = noIntegrationVariables;
45
46 ak2 = new G4double[numberOfVariables] ;
47 ak3 = new G4double[numberOfVariables] ;
48 ak4 = new G4double[numberOfVariables] ;
49 ak5 = new G4double[numberOfVariables] ;
50 ak6 = new G4double[numberOfVariables] ;
51 ak7 = new G4double[numberOfVariables] ;
52 ak8 = new G4double[numberOfVariables] ;
53
54 // Must ensure space extra 'state' variables exists - i.e. yIn[7]
55 //
56 const G4int numStateMax = std::max(GetNumberOfStateVariables(), 8);
57 const G4int numStateVars = std::max(noIntegrationVariables,
58 numStateMax );
59
60 yTemp = new G4double[numStateVars] ;
61 yIn = new G4double[numStateVars] ;
62
63 fLastInitialVector = new G4double[numberOfVariables] ;
64 fLastFinalVector = new G4double[numberOfVariables] ;
65
66 fLastDyDx = new G4double[numberOfVariables];
67
68 fMidVector = new G4double[numberOfVariables];
69 fMidError = new G4double[numberOfVariables];
70
71 if( primary )
72 {
73 fAuxStepper = new G4TsitourasRK45(EqRhs, numberOfVariables, !primary);
74 }
75}
G4double
double G4double
Definition G4Types.hh:83
G4int
int G4int
Definition G4Types.hh:85
G4MagIntegratorStepper::G4MagIntegratorStepper
G4MagIntegratorStepper(G4EquationOfMotion *Equation, G4int numIntegrationVariables, G4int numStateVariables=12, G4bool isFSAL=false)
Definition G4MagIntegratorStepper.cc:37
G4MagIntegratorStepper::GetNumberOfStateVariables
G4int GetNumberOfStateVariables() const
G4TsitourasRK45::G4TsitourasRK45
G4TsitourasRK45(G4EquationOfMotion *EqRhs, G4int numberOfVariables=6, G4bool primary=true)
Definition G4TsitourasRK45.cc:39

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

◆ ~G4TsitourasRK45()

G4TsitourasRK45::~G4TsitourasRK45 ( )
override

Destructor.

Definition at line 81 of file G4TsitourasRK45.cc.

82{
83 delete [] ak2;
84 delete [] ak3;
85 delete [] ak4;
86 delete [] ak5;
87 delete [] ak6;
88 delete [] ak7;
89 delete [] ak8;
90
91 delete [] yTemp;
92 delete [] yIn;
93
94 delete [] fLastInitialVector;
95 delete [] fLastFinalVector;
96 delete [] fLastDyDx;
97 delete [] fMidVector;
98 delete [] fMidError;
99
100 delete fAuxStepper;
101}

◆ G4TsitourasRK45() [2/2]

G4TsitourasRK45::G4TsitourasRK45 ( const G4TsitourasRK45 & )
delete

Copy constructor and assignment operator not allowed.

Member Function Documentation

◆ DistChord()

G4double G4TsitourasRK45::DistChord ( ) const
overridevirtual

Returns the distance from chord line.

Implements G4MagIntegratorStepper.

Definition at line 296 of file G4TsitourasRK45.cc.

297{
298 G4double distLine, distChord;
299 G4ThreeVector initialPoint, finalPoint, midPoint;
300
301 // Store last initial and final points (they will be
302 // overwritten in self-Stepper call!)
303 //
304 initialPoint = G4ThreeVector( fLastInitialVector[0],
305 fLastInitialVector[1], fLastInitialVector[2]);
306 finalPoint = G4ThreeVector( fLastFinalVector[0],
307 fLastFinalVector[1], fLastFinalVector[2]);
308
309 // Do half a step using StepNoErr
310 //
311 fAuxStepper->Stepper( fLastInitialVector, fLastDyDx, 0.5 * fLastStepLength,
312 fMidVector, fMidError );
313
314 midPoint = G4ThreeVector( fMidVector[0], fMidVector[1], fMidVector[2]);
315
316 // Use stored values of Initial and Endpoint + new Midpoint to evaluate
317 // distance of Chord
318 //
319 if (initialPoint != finalPoint)
320 {
321 distLine = G4LineSection::Distline( midPoint, initialPoint, finalPoint );
322 distChord = distLine;
323 }
324 else
325 {
326 distChord = (midPoint-initialPoint).mag();
327 }
328 return distChord;
329}
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 G4TsitourasRK45::IntegratorOrder ( ) const
inlineoverridevirtual

Returns the order, 4, of integration.

Implements G4MagIntegratorStepper.

Definition at line 125 of file G4TsitourasRK45.hh.

125{ return 4; }

◆ Interpolate()

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

Calculates the output at the tau fraction of step.

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

Definition at line 249 of file G4TsitourasRK45.cc.

254{
255 G4double bf1, bf2, bf3, bf4, bf5, bf6, bf7;
256 // Coefficients for all the seven stages.
257
258 const G4int numberOfVariables = GetNumberOfVariables();
259
260 G4double tau0 = tau;
261
262 for(G4int i=0; i<numberOfVariables; ++i)
263 {
264 yIn[i] = yInput[i];
265 }
266
267 G4double tau_2 = tau0*tau0 ;
268 // tau_3 = tau0*tau_2,
269 // tau_4 = tau_2*tau_2;
270
271 bf1 = -1.0530884977290216*tau*(tau - 1.3299890189751412)*(tau_2 -
272 1.4364028541716351*tau + 0.7139816917074209);
273 bf2 = 0.1017*tau_2*(tau_2 - 2.1966568338249754*tau +
274 1.2949852507374631);
275 bf3 = 2.490627285651252793*tau_2*(tau_2 - 2.38535645472061657*tau
276 + 1.57803468208092486);
277 bf4 = -16.54810288924490272*(tau - 1.21712927295533244)*
278 (tau - 0.61620406037800089)*tau_2;
279 bf5 = 47.37952196281928122*(tau - 1.203071208372362603)*
280 (tau - 0.658047292653547382)*tau_2;
281 bf6 = -34.87065786149660974*(tau - 1.2)*(tau -
282 0.666666666666666667)*tau_2;
283 bf7 = 2.5*(tau - 1.0)*(tau - 0.6)*tau_2;
284
285 // Putting together the coefficients calculated as the respective
286 // stage coefficients
287 //
288 for(G4int i=0; i<numberOfVariables; ++i)
289 {
290 yOut[i] = yIn[i] + Step*( bf1*dydx[i] + bf2*ak2[i] + bf3*ak3[i]
291 + bf4*ak4[i] + bf5*ak5[i] + bf6*ak6[i]
292 + bf7*ak7[i] ) ;
293 }
294}
G4MagIntegratorStepper::GetNumberOfVariables
G4int GetNumberOfVariables() const

◆ interpolate()

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

◆ operator=()

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

◆ SetupInterpolation()

void G4TsitourasRK45::SetupInterpolation ( )

Setup all coefficients for interpolation.

Definition at line 243 of file G4TsitourasRK45.cc.

245{
246 // Nothing to be done
247}

◆ Stepper()

void G4TsitourasRK45::Stepper ( const G4double y[],
const G4double dydx[],
G4double h,
G4double yout[],
G4double yerr[] )
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.

Implements G4MagIntegratorStepper.

Definition at line 116 of file G4TsitourasRK45.cc.

121{
122 const G4double b21 = 0.161 ,
123 b31 = -0.00848065549235698854 ,
124 b32 = 0.335480655492356989 ,
125
126 b41 = 2.89715305710549343 ,
127 b42 = -6.35944848997507484 ,
128 b43 = 4.36229543286958141 ,
129
130 b51 = 5.325864828439257,
131 b52 = -11.748883564062828,
132 b53 = 7.49553934288983621 ,
133 b54 = -0.09249506636175525,
134
135 b61 = 5.8614554429464200,
136 b62 = -12.9209693178471093 ,
137 b63 = 8.1593678985761586 ,
138 b64 = -0.071584973281400997,
139 b65 = -0.0282690503940683829,
140
141 b71 = 0.0964607668180652295 ,
142 b72 = 0.01,
143 b73 = 0.479889650414499575,
144 b74 = 1.37900857410374189,
145 b75 = -3.2900695154360807,
146 b76 = 2.32471052409977398,
147
148 // c1 = 0.001780011052226 ,
149 // c2 = 0.000816434459657 ,
150 // c3 = -0.007880878010262 ,
151 // c4 = 0.144711007173263 ,
152 // c5 = -0.582357165452555 ,
153 // c6 = 0.458082105929187 ,
154 // c7 = 1.0/66.0 ;
155
156 dc1 = 0.0935237485818927066 - b71 , // - 0.001780011052226,
157 dc2 = 0.00865288314156636761 - b72, // - 0.000816434459657,
158 dc3 = 0.492893099131431868 - b73 , // + 0.007880878010262,
159 dc4 = 1.14023541226785810 - b74 , // 0.144711007173263,
160 dc5 = - 2.3291801924393646 - b75, // + 0.582357165452555,
161 dc6 = 1.56887504931661552 - b76 , // - 0.458082105929187,
162 dc7 = 0.025; //- 1.0/66.0 ;
163
164 // dc1 = -3.0/1280.0,
165 // dc2 = 0.0,
166 // dc3 = 6561.0/632320.0,
167 // dc4 = -343.0/20800.0,
168 // dc5 = 243.0/12800.0,
169 // dc6 = -1.0/95.0,
170 // dc7 = 0.0 ;
171
172 const G4int numberOfVariables = GetNumberOfVariables();
173
174 // The number of variables to be integrated over
175 //
176 yOut[7] = yTemp[7] = yIn[7] = yInput[7];
177
178 // Saving yInput because yInput and yOut can be aliases for same array
179 //
180 for(G4int i=0; i<numberOfVariables; ++i)
181 {
182 yIn[i]=yInput[i];
183 }
184 // RightHandSide(yIn, dydx) ; // 1st Step - Not doing, getting passed
185
186 for(G4int i=0; i<numberOfVariables; ++i)
187 {
188 yTemp[i] = yIn[i] + b21*Step*dydx[i] ;
189 }
190 RightHandSide(yTemp, ak2) ; // 2nd Stage
191
192 for(G4int i=0; i<numberOfVariables; ++i)
193 {
194 yTemp[i] = yIn[i] + Step*(b31*dydx[i] + b32*ak2[i]) ;
195 }
196 RightHandSide(yTemp, ak3) ; // 3rd Stage
197
198 for(G4int i=0; i<numberOfVariables; ++i)
199 {
200 yTemp[i] = yIn[i] + Step*(b41*dydx[i] + b42*ak2[i] + b43*ak3[i]) ;
201 }
202 RightHandSide(yTemp, ak4) ; // 4th Stage
203
204 for(G4int i=0; i<numberOfVariables; ++i)
205 {
206 yTemp[i] = yIn[i] + Step*(b51*dydx[i] + b52*ak2[i] + b53*ak3[i] +
207 b54*ak4[i]) ;
208 }
209 RightHandSide(yTemp, ak5) ; // 5th Stage
210
211 for(G4int i=0; i<numberOfVariables; ++i)
212 {
213 yTemp[i] = yIn[i] + Step*(b61*dydx[i] + b62*ak2[i] + b63*ak3[i] +
214 b64*ak4[i] + b65*ak5[i]) ;
215 }
216 RightHandSide(yTemp, ak6) ; // 6th Stage
217
218 for(G4int i=0; i<numberOfVariables; ++i)
219 {
220 yOut[i] = yIn[i] + Step*(b71*dydx[i] + b72*ak2[i] + b73*ak3[i] +
221 b74*ak4[i] + b75*ak5[i] + b76*ak6[i]);
222 }
223 RightHandSide(yOut, ak7); // 7th Stage
224
225 //Calculate the error in the step:
226 for(G4int i=0; i<numberOfVariables; ++i)
227 {
228 yErr[i] = Step*(dc1*dydx[i] + dc2*ak2[i] + dc3*ak3[i] + dc4*ak4[i] +
229 dc5*ak5[i] + dc6*ak6[i] + dc7*ak7[i] ) ;
230
231 // Store Input and Final values, for possible use in calculating chord
232 //
233 fLastInitialVector[i] = yIn[i] ;
234 fLastFinalVector[i] = yOut[i];
235 fLastDyDx[i] = dydx[i];
236 }
237
238 fLastStepLength = Step;
239
240 return ;
241}
G4MagIntegratorStepper::RightHandSide
void RightHandSide(const G4double y[], G4double dydx[]) const

◆ StepperType()

G4StepperType G4TsitourasRK45::StepperType ( ) const
inlineoverridevirtual

Returns the stepper type-ID, "kTsitourasRK45".

Reimplemented from G4MagIntegratorStepper.

Definition at line 130 of file G4TsitourasRK45.hh.

130{ return kTsitourasRK45; }
kTsitourasRK45
@ kTsitourasRK45
G4TsitourasRK45.
Definition G4FieldParameters.hh:112
The documentation for this class was generated from the following files: