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

#include <G4DormandPrince745.hh>

+ Inheritance diagram for G4DormandPrince745:

Public Member Functions

 G4DormandPrince745 (G4EquationOfMotion *equation, G4int numberOfVariables=6)
 
void Stepper (const G4double yInput[], const G4double dydx[], G4double hstep, G4double yOutput[], G4double yError[]) override
 
void Stepper (const G4double yInput[], const G4double dydx[], G4double hstep, G4double yOutput[], G4double yError[], G4double dydxOutput[])
 
void SetupInterpolation ()
 
void Interpolate (G4double tau, G4double yOut[]) const
 
G4double DistChord () const override
 
G4int IntegratorOrder () const override
 
const G4String & StepperType () const
 
const G4String & StepperDescription () const
 
const field_utils::State & GetYOut () const
 
void Interpolate4thOrder (G4double yOut[], G4double tau) const
 
void SetupInterpolation5thOrder ()
 
void Interpolate5thOrder (G4double yOut[], G4double tau) const
 
G4EquationOfMotion * GetSpecificEquation ()
 
- 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

Definition at line 46 of file G4DormandPrince745.hh.

Constructor & Destructor Documentation

◆ G4DormandPrince745()

G4DormandPrince745::G4DormandPrince745 ( G4EquationOfMotion * equation,
G4int numberOfVariables = 6 )

Definition at line 75 of file G4DormandPrince745.cc.

77 : G4MagIntegratorStepper(equation, noIntegrationVariables)
78{
79}
G4MagIntegratorStepper::G4MagIntegratorStepper
G4MagIntegratorStepper(G4EquationOfMotion *Equation, G4int numIntegrationVariables, G4int numStateVariables=12, G4bool isFSAL=false)
Definition G4MagIntegratorStepper.cc:38

Member Function Documentation

◆ DistChord()

G4double G4DormandPrince745::DistChord ( ) const
overridevirtual

Implements G4MagIntegratorStepper.

Definition at line 214 of file G4DormandPrince745.cc.

215{
216 // Coefficients were taken from Some Practical Runge-Kutta Formulas
217 // by Lawrence F. Shampine, page 149, c*
218 //
219 const G4double hf1 = 6025192743.0 / 30085553152.0,
220 hf3 = 51252292925.0 / 65400821598.0,
221 hf4 = - 2691868925.0 / 45128329728.0,
222 hf5 = 187940372067.0 / 1594534317056.0,
223 hf6 = - 1776094331.0 / 19743644256.0,
224 hf7 = 11237099.0 / 235043384.0;
225
226 G4ThreeVector mid;
227
228 for(G4int i = 0; i < 3; ++i)
229 {
230 mid[i] = fyIn[i] + 0.5 * fLastStepLength * (
231 hf1 * fdydxIn[i] + hf3 * ak3[i] +
232 hf4 * ak4[i] + hf5 * ak5[i] + hf6 * ak6[i] + hf7 * ak7[i]);
233 }
234
235 const G4ThreeVector begin = makeVector(fyIn, Value3D::Position);
236 const G4ThreeVector end = makeVector(fyOut, Value3D::Position);
237
238 return G4LineSection::Distline(mid, begin, end);
239}
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition G4ThreeVector.hh:36
G4double
double G4double
Definition G4Types.hh:83
G4int
int G4int
Definition G4Types.hh:85
G4LineSection::Distline
static G4double Distline(const G4ThreeVector &OtherPnt, const G4ThreeVector &LinePntA, const G4ThreeVector &LinePntB)
Definition G4LineSection.hh:72
field_utils::makeVector
G4ThreeVector makeVector(const ArrayType &array, Value3D value)

◆ GetSpecificEquation()

G4EquationOfMotion * G4DormandPrince745::GetSpecificEquation ( )
inline

Definition at line 88 of file G4DormandPrince745.hh.

88{ return GetEquationOfMotion(); }
G4MagIntegratorStepper::GetEquationOfMotion
G4EquationOfMotion * GetEquationOfMotion()

◆ GetYOut()

const field_utils::State & G4DormandPrince745::GetYOut ( ) const
inline

Definition at line 81 of file G4DormandPrince745.hh.

81{ return fyOut; }

◆ IntegratorOrder()

G4int G4DormandPrince745::IntegratorOrder ( ) const
inlineoverridevirtual

Implements G4MagIntegratorStepper.

Definition at line 76 of file G4DormandPrince745.hh.

76{ return 4; }

◆ Interpolate()

void G4DormandPrince745::Interpolate ( G4double tau,
G4double yOut[] ) const
inline

Definition at line 68 of file G4DormandPrince745.hh.

69 {
70 Interpolate4thOrder(yOut, tau);
71 }
G4DormandPrince745::Interpolate4thOrder
void Interpolate4thOrder(G4double yOut[], G4double tau) const
Definition G4DormandPrince745.cc:247

◆ Interpolate4thOrder()

void G4DormandPrince745::Interpolate4thOrder ( G4double yOut[],
G4double tau ) const

Definition at line 246 of file G4DormandPrince745.cc.

248{
249 const G4int numberOfVariables = GetNumberOfVariables();
250
251 const G4double tau2 = tau * tau,
252 tau3 = tau * tau2,
253 tau4 = tau2 * tau2;
254
255 const G4double bf1 = 1.0 / 11282082432.0 * (
256 157015080.0 * tau4 - 13107642775.0 * tau3 + 34969693132.0 * tau2 -
257 32272833064.0 * tau + 11282082432.0);
258
259 const G4double bf3 = - 100.0 / 32700410799.0 * tau * (
260 15701508.0 * tau3 - 914128567.0 * tau2 + 2074956840.0 * tau -
261 1323431896.0);
262
263 const G4double bf4 = 25.0 / 5641041216.0 * tau * (
264 94209048.0 * tau3 - 1518414297.0 * tau2 + 2460397220.0 * tau -
265 889289856.0);
266
267 const G4double bf5 = - 2187.0 / 199316789632.0 * tau * (
268 52338360.0 * tau3 - 451824525.0 * tau2 + 687873124.0 * tau -
269 259006536.0);
270
271 const G4double bf6 = 11.0 / 2467955532.0 * tau * (
272 106151040.0 * tau3 - 661884105.0 * tau2 +
273 946554244.0 * tau - 361440756.0);
274
275 const G4double bf7 = 1.0 / 29380423.0 * tau * (1.0 - tau) * (
276 8293050.0 * tau2 - 82437520.0 * tau + 44764047.0);
277
278 for(G4int i = 0; i < numberOfVariables; ++i)
279 {
280 yOut[i] = fyIn[i] + fLastStepLength * tau * (
281 bf1 * fdydxIn[i] + bf3 * ak3[i] + bf4 * ak4[i] +
282 bf5 * ak5[i] + bf6 * ak6[i] + bf7 * ak7[i]);
283 }
284}
G4MagIntegratorStepper::GetNumberOfVariables
G4int GetNumberOfVariables() const

Referenced by Interpolate().

◆ Interpolate5thOrder()

void G4DormandPrince745::Interpolate5thOrder ( G4double yOut[],
G4double tau ) const

Definition at line 342 of file G4DormandPrince745.cc.

344{
345 // Define the coefficients for the polynomials
346 //
347 G4double bi[10][5];
348
349 // COEFFICIENTS OF bi[1]
350 bi[1][0] = 1.0,
351 bi[1][1] = -38039.0 / 7040.0,
352 bi[1][2] = 125923.0 / 10560.0,
353 bi[1][3] = -19683.0 / 1760.0,
354 bi[1][4] = 3303.0 / 880.0,
355 // --------------------------------------------------------
356 //
357 // COEFFICIENTS OF bi[2]
358 bi[2][0] = 0.0,
359 bi[2][1] = 0.0,
360 bi[2][2] = 0.0,
361 bi[2][3] = 0.0,
362 bi[2][4] = 0.0,
363 // --------------------------------------------------------
364 //
365 // COEFFICIENTS OF bi[3]
366 bi[3][0] = 0.0,
367 bi[3][1] = -12500.0 / 4081.0,
368 bi[3][2] = 205000.0 / 12243.0,
369 bi[3][3] = -90000.0 / 4081.0,
370 bi[3][4] = 36000.0 / 4081.0,
371 // --------------------------------------------------------
372 //
373 // COEFFICIENTS OF bi[4]
374 bi[4][0] = 0.0,
375 bi[4][1] = -3125.0 / 704.0,
376 bi[4][2] = 25625.0 / 1056.0,
377 bi[4][3] = -5625.0 / 176.0,
378 bi[4][4] = 1125.0 / 88.0,
379 // --------------------------------------------------------
380 //
381 // COEFFICIENTS OF bi[5]
382 bi[5][0] = 0.0,
383 bi[5][1] = 164025.0 / 74624.0,
384 bi[5][2] = -448335.0 / 37312.0,
385 bi[5][3] = 295245.0 / 18656.0,
386 bi[5][4] = -59049.0 / 9328.0,
387 // --------------------------------------------------------
388 //
389 // COEFFICIENTS OF bi[6]
390 bi[6][0] = 0.0,
391 bi[6][1] = -25.0 / 28.0,
392 bi[6][2] = 205.0 / 42.0,
393 bi[6][3] = -45.0 / 7.0,
394 bi[6][4] = 18.0 / 7.0,
395 // --------------------------------------------------------
396 //
397 // COEFFICIENTS OF bi[7]
398 bi[7][0] = 0.0,
399 bi[7][1] = -2.0 / 11.0,
400 bi[7][2] = 73.0 / 55.0,
401 bi[7][3] = -171.0 / 55.0,
402 bi[7][4] = 108.0 / 55.0,
403 // --------------------------------------------------------
404 //
405 // COEFFICIENTS OF bi[8]
406 bi[8][0] = 0.0,
407 bi[8][1] = 189.0 / 22.0,
408 bi[8][2] = -1593.0 / 55.0,
409 bi[8][3] = 3537.0 / 110.0,
410 bi[8][4] = -648.0 / 55.0,
411 // --------------------------------------------------------
412 //
413 // COEFFICIENTS OF bi[9]
414 bi[9][0] = 0.0,
415 bi[9][1] = 351.0 / 110.0,
416 bi[9][2] = -999.0 / 55.0,
417 bi[9][3] = 2943.0 / 110.0,
418 bi[9][4] = -648.0 / 55.0;
419 // --------------------------------------------------------
420
421 // Calculating the polynomials
422
423 G4double b[10];
424 std::memset(b, 0.0, sizeof(b));
425
426 G4double tauPower = 1.0;
427 for(G4int j = 0; j <= 4; ++j)
428 {
429 for(G4int iStage = 1; iStage <= 9; ++iStage)
430 {
431 b[iStage] += bi[iStage][j] * tauPower;
432 }
433 tauPower *= tau;
434 }
435
436 const G4int numberOfVariables = GetNumberOfVariables();
437 const G4double stepLen = fLastStepLength * tau;
438 for(G4int i = 0; i < numberOfVariables; ++i)
439 {
440 yOut[i] = fyIn[i] + stepLen * (
441 b[1] * fdydxIn[i] + b[2] * ak2[i] + b[3] * ak3[i] +
442 b[4] * ak4[i] + b[5] * ak5[i] + b[6] * ak6[i] +
443 b[7] * ak7[i] + b[8] * ak8[i] + b[9] * ak9[i]
444 );
445 }
446}

◆ SetupInterpolation()

void G4DormandPrince745::SetupInterpolation ( )
inline

Definition at line 66 of file G4DormandPrince745.hh.

66{}

◆ SetupInterpolation5thOrder()

void G4DormandPrince745::SetupInterpolation5thOrder ( )

Definition at line 293 of file G4DormandPrince745.cc.

294{
295 // Coefficients for the additional stages
296 //
297 const G4double b81 = 6245.0 / 62208.0,
298 b82 = 0.0,
299 b83 = 8875.0 / 103032.0,
300 b84 = -125.0 / 1728.0,
301 b85 = 801.0 / 13568.0,
302 b86 = -13519.0 / 368064.0,
303 b87 = 11105.0 / 368064.0,
304
305 b91 = 632855.0 / 4478976.0,
306 b92 = 0.0,
307 b93 = 4146875.0 / 6491016.0,
308 b94 = 5490625.0 /14183424.0,
309 b95 = -15975.0 / 108544.0,
310 b96 = 8295925.0 / 220286304.0,
311 b97 = -1779595.0 / 62938944.0,
312 b98 = -805.0 / 4104.0;
313
314 const G4int numberOfVariables = GetNumberOfVariables();
315 State yTemp = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
316
317 // Evaluate the extra stages
318 //
319 for(G4int i = 0; i < numberOfVariables; ++i)
320 {
321 yTemp[i] = fyIn[i] + fLastStepLength * (
322 b81 * fdydxIn[i] + b82 * ak2[i] + b83 * ak3[i] +
323 b84 * ak4[i] + b85 * ak5[i] + b86 * ak6[i] +
324 b87 * ak7[i]
325 );
326 }
327 RightHandSide(yTemp, ak8); // 8th Stage
328
329 for(G4int i = 0; i < numberOfVariables; ++i)
330 {
331 yTemp[i] = fyIn[i] + fLastStepLength * (
332 b91 * fdydxIn[i] + b92 * ak2[i] + b93 * ak3[i] +
333 b94 * ak4[i] + b95 * ak5[i] + b96 * ak6[i] +
334 b97 * ak7[i] + b98 * ak8[i]
335 );
336 }
337 RightHandSide(yTemp, ak9); // 9th Stage
338}
State
#define State(X)
Definition G4ITMultiNavigator.cc:45
G4MagIntegratorStepper::RightHandSide
void RightHandSide(const G4double y[], G4double dydx[]) const

◆ Stepper() [1/2]

void G4DormandPrince745::Stepper ( const G4double yInput[],
const G4double dydx[],
G4double hstep,
G4double yOutput[],
G4double yError[] )
overridevirtual

Implements G4MagIntegratorStepper.

Definition at line 97 of file G4DormandPrince745.cc.

102{
103 // The various constants defined on the basis of butcher tableu
104 //
105 const G4double b21 = 0.2,
106 b31 = 3.0 / 40.0, b32 = 9.0 / 40.0,
107 b41 = 44.0 / 45.0, b42 = -56.0 / 15.0, b43 = 32.0/9.0,
108
109 b51 = 19372.0 / 6561.0, b52 = -25360.0 / 2187.0, b53 = 64448.0 / 6561.0,
110 b54 = -212.0 / 729.0,
111
112 b61 = 9017.0 / 3168.0 , b62 = -355.0 / 33.0,
113 b63 = 46732.0 / 5247.0, b64 = 49.0 / 176.0,
114 b65 = -5103.0 / 18656.0,
115
116 b71 = 35.0 / 384.0, b72 = 0.,
117 b73 = 500.0 / 1113.0, b74 = 125.0 / 192.0,
118 b75 = -2187.0 / 6784.0, b76 = 11.0 / 84.0,
119
120 //Sum of columns, sum(bij) = ei
121 // e1 = 0. ,
122 // e2 = 1.0/5.0 ,
123 // e3 = 3.0/10.0 ,
124 // e4 = 4.0/5.0 ,
125 // e5 = 8.0/9.0 ,
126 // e6 = 1.0 ,
127 // e7 = 1.0 ,
128
129 // Difference between the higher and the lower order method coeff. :
130 // b7j are the coefficients of higher order
131
132 dc1 = -(b71 - 5179.0 / 57600.0),
133 dc2 = -(b72 - .0),
134 dc3 = -(b73 - 7571.0 / 16695.0),
135 dc4 = -(b74 - 393.0 / 640.0),
136 dc5 = -(b75 + 92097.0 / 339200.0),
137 dc6 = -(b76 - 187.0 / 2100.0),
138 dc7 = -(- 1.0 / 40.0);
139
140 const G4int numberOfVariables = GetNumberOfVariables();
141 State yTemp = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
142
143 // The number of variables to be integrated over
144 //
145 yOut[7] = yTemp[7] = yInput[7];
146
147 // Saving yInput because yInput and yOut can be aliases for same array
148 //
149 for(G4int i = 0; i < numberOfVariables; ++i)
150 {
151 fyIn[i] = yInput[i];
152 }
153 // RightHandSide(yIn, dydx); // Not done! 1st stage
154
155 for(G4int i = 0; i < numberOfVariables; ++i)
156 {
157 yTemp[i] = fyIn[i] + b21 * hstep * dydx[i];
158 }
159 RightHandSide(yTemp, ak2); // 2nd stage
160
161 for(G4int i = 0; i < numberOfVariables; ++i)
162 {
163 yTemp[i] = fyIn[i] + hstep * (b31 * dydx[i] + b32 * ak2[i]);
164 }
165 RightHandSide(yTemp, ak3); // 3rd stage
166
167 for(G4int i = 0; i < numberOfVariables; ++i)
168 {
169 yTemp[i] = fyIn[i] + hstep * (
170 b41 * dydx[i] + b42 * ak2[i] + b43 * ak3[i]);
171 }
172 RightHandSide(yTemp, ak4); // 4th stage
173
174 for(G4int i = 0; i < numberOfVariables; ++i)
175 {
176 yTemp[i] = fyIn[i] + hstep * (
177 b51 * dydx[i] + b52 * ak2[i] + b53 * ak3[i] + b54 * ak4[i]);
178 }
179 RightHandSide(yTemp, ak5); // 5th stage
180
181 for(G4int i = 0; i < numberOfVariables; ++i)
182 {
183 yTemp[i] = fyIn[i] + hstep * (
184 b61 * dydx[i] + b62 * ak2[i] +
185 b63 * ak3[i] + b64 * ak4[i] + b65 * ak5[i]);
186 }
187 RightHandSide(yTemp, ak6); // 6th stage
188
189 for(G4int i = 0; i < numberOfVariables; ++i)
190 {
191 yOut[i] = fyIn[i] + hstep * (
192 b71 * dydx[i] + b72 * ak2[i] + b73 * ak3[i] +
193 b74 * ak4[i] + b75 * ak5[i] + b76 * ak6[i]);
194 }
195 RightHandSide(yOut, ak7); // 7th and Final stage
196
197 for(G4int i = 0; i < numberOfVariables; ++i)
198 {
199 yErr[i] = hstep * (
200 dc1 * dydx[i] + dc2 * ak2[i] +
201 dc3 * ak3[i] + dc4 * ak4[i] +
202 dc5 * ak5[i] + dc6 * ak6[i] + dc7 * ak7[i]
203 ) + 1.5e-18;
204
205 // Store Input and Final values, for possible use in calculating chord
206 //
207 fyOut[i] = yOut[i];
208 fdydxIn[i] = dydx[i];
209 }
210
211 fLastStepLength = hstep;
212}

Referenced by Stepper().

◆ Stepper() [2/2]

void G4DormandPrince745::Stepper ( const G4double yInput[],
const G4double dydx[],
G4double hstep,
G4double yOutput[],
G4double yError[],
G4double dydxOutput[] )

Definition at line 81 of file G4DormandPrince745.cc.

87{
88 Stepper(yInput, dydx, hstep, yOutput, yError);
89 field_utils::copy(dydxOutput, ak7);
90}
G4DormandPrince745::Stepper
void Stepper(const G4double yInput[], const G4double dydx[], G4double hstep, G4double yOutput[], G4double yError[]) override
Definition G4DormandPrince745.cc:97
field_utils::copy
void copy(G4double dst[], const G4double src[], std::size_t size=G4FieldTrack::ncompSVEC)
Definition G4FieldUtils.cc:98

◆ StepperDescription()

const G4String & G4DormandPrince745::StepperDescription ( ) const

Definition at line 68 of file G4DormandPrince745.cc.

69{
70 static G4String _stepperDescription(
71 "Embedeed 5th order Runge-Kutta stepper - 7 stages, FSAL, Interpolating.");
72 return _stepperDescription;
73}

◆ StepperType()

const G4String & G4DormandPrince745::StepperType ( ) const

Definition at line 61 of file G4DormandPrince745.cc.

62{
63 static G4String _stepperType("G4DormandPrince745: 5th order");
64 return _stepperType;
65}
The documentation for this class was generated from the following files: