Geant4 11.2.2
Toolkit for the simulation of the passage of particles through matter
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
G4FSALDormandPrince745.cc
Go to the documentation of this file.
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26// G4FSALDormandPrince745 implementation
27//
28// The Butcher table of the FDormand-Prince-7-4-5 method is as follows:
29//
30// 0 |
31// 1/5 | 1/5
32// 3/10| 3/40 9/40
33// 4/5 | 44/45 56/15 32/9
34// 8/9 | 19372/6561 25360/2187 64448/6561 212/729
35// 1 | 9017/3168 355/33 46732/5247 49/176 5103/18656
36// 1 | 35/384 0 500/1113 125/192 2187/6784 11/84
37// ---------------------------------------------------------------------------
38// 35/384 0 500/1113 125/192 2187/6784 11/84 0
39// 5179/57600 0 7571/16695 393/640 92097/339200 187/2100 1/40
40//
41// Created: Somnath Banerjee, Google Summer of Code 2015, 25 May 2015
42// Supervision: John Apostolakis, CERN
43// --------------------------------------------------------------------
44
45#include "G4FSALDormandPrince745.hh"
46#include "G4LineSection.hh"
47#include <cmath>
48
49// Constructor
50//
51G4FSALDormandPrince745::G4FSALDormandPrince745(G4EquationOfMotion* EqRhs,
52 G4int noIntegrationVariables,
53 G4bool primary)
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}
92
93// Destructor
94//
95G4FSALDormandPrince745::~G4FSALDormandPrince745()
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}
122
123// Stepper
124//
125// Passing in the value of yInput[],the first time dydx[] and Step length
126// Giving back yOut and yErr arrays for output and error respectively
127//
128void G4FSALDormandPrince745::Stepper(const G4double yInput[],
129 const G4double dydx[],
130 G4double Step,
131 G4double yOut[],
132 G4double yErr[],
133 G4double nextDydx[] )
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}
240
241// DistChord
242//
243G4double G4FSALDormandPrince745::DistChord() const
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}
277
278// interpolate
279//
280void G4FSALDormandPrince745::interpolate( const G4double yInput[],
281 const G4double dydx[],
282 G4double yOut[],
283 G4double Step,
284 G4double tau)
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}
326
327// SetupInterpolate
328//
329void G4FSALDormandPrince745::SetupInterpolate(const G4double yInput[],
330 const G4double dydx[],
331 const G4double Step )
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}
381
382// Interpolate
383//
384void G4FSALDormandPrince745::Interpolate( const G4double yInput[],
385 const G4double dydx[],
386 const G4double Step,
387 G4double yOut[],
388 G4double tau )
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}
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition G4ThreeVector.hh:36
G4double
double G4double
Definition G4Types.hh:83
G4bool
bool G4bool
Definition G4Types.hh:86
G4int
int G4int
Definition G4Types.hh:85
G4FSALDormandPrince745::G4FSALDormandPrince745
G4FSALDormandPrince745(G4EquationOfMotion *EqRhs, G4int numberOfVariables=6, G4bool primary=true)
Definition G4FSALDormandPrince745.cc:51
G4FSALDormandPrince745::Interpolate
void Interpolate(const G4double yInput[], const G4double dydx[], const G4double Step, G4double yOut[], G4double tau)
Definition G4FSALDormandPrince745.cc:384
G4FSALDormandPrince745::Stepper
void Stepper(const G4double y[], const G4double dydx[], G4double h, G4double yout[], G4double yerr[], G4double nextDydx[]) override
Definition G4FSALDormandPrince745.cc:128
G4FSALDormandPrince745::interpolate
void interpolate(const G4double yInput[], const G4double dydx[], G4double yOut[], G4double Step, G4double tau)
Definition G4FSALDormandPrince745.cc:280
G4FSALDormandPrince745::DistChord
G4double DistChord() const override
Definition G4FSALDormandPrince745.cc:243
G4FSALDormandPrince745::SetupInterpolate
void SetupInterpolate(const G4double yInput[], const G4double dydx[], const G4double Step)
Definition G4FSALDormandPrince745.cc:329
G4LineSection::Distline
static G4double Distline(const G4ThreeVector &OtherPnt, const G4ThreeVector &LinePntA, const G4ThreeVector &LinePntB)
Definition G4LineSection.hh:72
G4VFSALIntegrationStepper::GetNumberOfVariables
G4int GetNumberOfVariables() const
G4VFSALIntegrationStepper::RightHandSide
void RightHandSide(const double y[], double dydx[])
Definition G4VFSALIntegrationStepper.cc:51