| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // G4FSALDormandPrince745 implementation
27 //
28 // The Butcher table of the FDormand-Prince-7-
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 2
35 // 1 | 9017/3168 355/33 46732/5247 4
36 // 1 | 35/384 0 500/1113 1
37 // ------------------------------------------
38 // 35/384 0 500/1113 1
39 // 5179/57600 0 7571/16695 3
40 //
41 // Created: Somnath Banerjee, Google Summer of
42 // Supervision: John Apostolakis, CERN
43 // -------------------------------------------
44
45 #include "G4FSALDormandPrince745.hh"
46 #include "G4LineSection.hh"
47 #include <cmath>
48
49 // Constructor
50 //
51 G4FSALDormandPrince745::G4FSALDormandPrince745
52
53
54 : G4VFSALIntegrationStepper(EqRhs, noIntegra
55 {
56 const G4int numberOfVariables = noIntegratio
57
58 // New Chunk of memory being created for use
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 interpola
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[numb
78
79 fInitialDyDx = new G4double[numberOfVariable
80 fLastInitialVector = new G4double[numberOfVa
81 fLastFinalVector = new G4double[numberOfVari
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(E
90 }
91 }
92
93 // Destructor
94 //
95 G4FSALDormandPrince745::~G4FSALDormandPrince74
96 {
97 // Clear all previously allocated memory f
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; pseud
112 delete [] fInitialDyDx; fInit
113
114 delete [] fLastInitialVector; fLastInit
115 delete [] fLastFinalVector; fLastFina
116 delete [] fLastDyDx; fLastDyDx
117 delete [] fMidVector; fMidVecto
118 delete [] fMidError; fMidError
119
120 delete fAuxStepper; fAuxStepper = nullptr;
121 }
122
123 // Stepper
124 //
125 // Passing in the value of yInput[],the first
126 // Giving back yOut and yErr arrays for output
127 //
128 void G4FSALDormandPrince745::Stepper(const G4d
129 const G4d
130 G4d
131 G4d
132 G4d
133 G4d
134 {
135 G4int i;
136
137 // The various constants defined on the ba
138
139 const G4double b21 = 0.2 ,
140 b31 = 3.0/40.0, b32 = 9.0/4
141
142 b41 = 44.0/45.0, b42 = -56.
143
144 b51 = 19372.0/6561.0, b52 =
145 b53 = 64448.0/6561.0, b54 =
146
147 b61 = 9017.0/3168.0 , b62 =
148 b63 = 46732.0/5247.0 ,
149 b65 = -5103.0/18656.0 ,
150
151 b71 = 35.0/384.0, b72 = 0.,
152 b73 = 500.0/1113.0, b74 = 1
153 b75 = -2187.0/6784.0, b76 =
154
155 // c1 = 35.0/384.0, c2 = .0,
156 // c3 = 500.0/1113.0, c4 = 125
157 // c5 = -2187.0/6784.0, c6 = 1
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.
165 dc6 = b76 - 187.0/2100.0,
166 dc7 = - 1.0/40.0 ; //end o
167
168 const G4int numberOfVariables = GetNumberO
169 // The number of variables to be integra
170
171 // Saving yInput because yInput and yOut c
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 ca
179 //
180 yOut[7] = yTemp[7] = yInput[7];
181 // RightHandSide(yIn, DyDx) ; // 1st Ste
182
183 for(i=0; i<numberOfVariables; ++i)
184 {
185 yTemp[i] = yIn[i] + b21*Step*fInitialD
186 }
187 RightHandSide(yTemp, ak2) ; /
188
189 for(i=0; i<numberOfVariables; ++i)
190 {
191 yTemp[i] = yIn[i] + Step*(b31*fInitial
192 }
193 RightHandSide(yTemp, ak3) ; /
194
195 for(i=0; i<numberOfVariables; ++i)
196 {
197 yTemp[i] = yIn[i] + Step*(b41*fInitial
198 + b42*ak2[i] +
199 }
200 RightHandSide(yTemp, ak4) ; /
201
202 for(i=0; i<numberOfVariables; ++i)
203 {
204 yTemp[i] = yIn[i] + Step*(b51*fInitial
205 + b52*ak2[i] +
206 }
207 RightHandSide(yTemp, ak5) ; /
208
209 for(i=0; i<numberOfVariables; ++i)
210 {
211 yTemp[i] = yIn[i] + Step*(b61*fInitial
212 + b63*ak3[i] +
213 }
214 RightHandSide(yTemp, ak6) ; /
215
216 for(i=0; i<numberOfVariables; ++i)
217 {
218 yOut[i] = yIn[i] + Step*(b71*fInitialD
219 + b74*ak4[i] +
220 }
221 RightHandSide(yOut, ak7); //
222
223 for(i=0; i<numberOfVariables; ++i)
224 {
225
226 yErr[i] = Step*(dc1*fInitialDyDx[i] +
227 + dc4*ak4[i] + dc5*ak5[i
228
229 // Store Input and Final values, for p
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 //
243 G4double G4FSALDormandPrince745::DistChord() c
244 {
245 G4double distLine, distChord;
246 G4ThreeVector initialPoint, finalPoint, mi
247
248 // Store last initial and final points
249 // (they will be overwritten in self-Stepp
250 //
251 initialPoint = G4ThreeVector(fLastInitialV
252 fLastInitialV
253 finalPoint = G4ThreeVector(fLastFinalVec
254 fLastFinalVec
255
256 // Do half a step using StepNoErr
257
258 fAuxStepper->Stepper( fLastInitialVector,
259 fMidVector, fMidErro
260
261 midPoint = G4ThreeVector( fMidVector[0], f
262
263 // Use stored values of Initial and Endpoi
264 // distance of Chord
265 //
266 if (initialPoint != finalPoint)
267 {
268 distLine = G4LineSection::Distline( m
269 distChord = distLine;
270 }
271 else
272 {
273 distChord = (midPoint-initialPoint).ma
274 }
275 return distChord;
276 }
277
278 // interpolate
279 //
280 void G4FSALDormandPrince745::interpolate( cons
281 cons
282
283
284
285 {
286 G4double bf1, bf2, bf3, bf4, bf5, bf6, bf7
287
288 const G4int numberOfVariables = GetNumberO
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*t
302 + 34969693132.0*tau_2- 32272833064.0*
303 / 11282082432.0;
304 bf2 = 0.0;
305 bf3 = - 100.0*tau*(15701508.0*tau_3 - 9141
306 + 2074956840.0*tau - 1323
307 bf4 = 25.0*tau*(94209048.0*tau_3- 15184142
308 + 2460397220.0*tau - 8892898
309 / 5641041216.0;
310 bf5 = -2187.0*tau*(52338360.0*tau_3 - 4518
311 + 687873124.0*tau - 25900
312 / 199316789632.0;
313 bf6 = 11.0*tau*(106151040.0*tau_3- 661884
314 + 946554244.0*tau - 3614407
315 / 2467955532.0;
316 bf7 = tau*(1.0 - tau)*(8293050.0*tau_2 - 8
317 / 29380423.0;
318
319 for(G4int i=0; i<numberOfVariables; ++i)
320 {
321 yOut[i] = yIn[i] + Step*tau*(bf1*dydx[i]
322 + bf4*ak4[i]
323 + bf7*ak7[i]
324 }
325 }
326
327 // SetupInterpolate
328 //
329 void G4FSALDormandPrince745::SetupInterpolate(
330
331
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 = GetNumberO
353
354 // Saving yInput because yInput and yOut c
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]
368 b84*ak4[i]
369 b87*ak7[i]
370 }
371 RightHandSide( yTemp, ak8 );
372
373 for(G4int i=0; i<numberOfVariables; ++i)
374 {
375 yTemp[i] = yIn[i] + Step * ( b91*dydx[
376 b94*ak4[i]
377 b97*ak7[i]
378 }
379 RightHandSide( yTemp, ak9 );
380 }
381
382 // Interpolate
383 //
384 void G4FSALDormandPrince745::Interpolate( cons
385 cons
386 cons
387
388
389 {
390 // Define the coefficients for the polynom
391
392 G4double bi[10][5], b[10];
393 G4int numberOfVariables = GetNumberOfVaria
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
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) /
488 {
489 yOut[i] = yIn[i] + Step*tau0*(b[1]*dyd
490 b[4]*ak4
491 b[7]*ak7
492 }
493 }
494