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
G4DormandPrince745.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// G4DormandPrince745 implementation
27//
28// DormandPrince7 - 5(4) non-FSAL
29// definition of the stepper() method that evaluates one step in
30// field propagation.
31// The coefficients and the algorithm have been adapted from
32//
33// J. R. Dormand and P. J. Prince, "A family of embedded Runge-Kutta formulae"
34// Journal of computational and applied Math., vol.6, no.1, pp.19-26, 1980.
35//
36// The Butcher table of the Dormand-Prince-7-4-5 method is as follows :
37//
38// 0 |
39// 1/5 | 1/5
40// 3/10| 3/40 9/40
41// 4/5 | 44/45 56/15 32/9
42// 8/9 | 19372/6561 25360/2187 64448/6561 212/729
43// 1 | 9017/3168 355/33 46732/5247 49/176 5103/18656
44// 1 | 35/384 0 500/1113 125/192 2187/6784 11/84
45// ------------------------------------------------------------------------
46// 35/384 0 500/1113 125/192 2187/6784 11/84 0
47// 5179/57600 0 7571/16695 393/640 92097/339200 187/2100 1/40
48//
49// Created: Somnath Banerjee, Google Summer of Code 2015, 25 May 2015
50// Supervision: John Apostolakis, CERN
51// --------------------------------------------------------------------
52
53#include "G4DormandPrince745.hh"
54#include "G4LineSection.hh"
55
56#include <cstring>
57
58using namespace field_utils;
59
60const G4String G4DormandPrince745::gStepperType =
61 G4String("G4DormandPrince745: 5th order");
62
63const G4String G4DormandPrince745::gStepperDescription= G4String(
64 "Embedeed 5th order Runge-Kutta stepper - 7 stages, FSAL, Interpolating.");
65
66
68 G4int noIntegrationVariables)
69 : G4MagIntegratorStepper(equation, noIntegrationVariables)
70{
71}
72
74 const G4double dydx[],
75 G4double hstep,
76 G4double yOutput[],
77 G4double yError[],
78 G4double dydxOutput[])
79{
80 Stepper(yInput, dydx, hstep, yOutput, yError);
81 field_utils::copy(dydxOutput, ak7);
82}
83
84// Stepper
85//
86// Passing in the value of yInput[],the first time dydx[] and Step length
87// Giving back yOut and yErr arrays for output and error respectively
88//
90 const G4double dydx[],
91 G4double hstep,
92 G4double yOut[],
93 G4double yErr[])
94{
95 // The various constants defined on the basis of butcher tableu
96 //
97 const G4double b21 = 0.2,
98 b31 = 3.0 / 40.0, b32 = 9.0 / 40.0,
99 b41 = 44.0 / 45.0, b42 = -56.0 / 15.0, b43 = 32.0/9.0,
100
101 b51 = 19372.0 / 6561.0, b52 = -25360.0 / 2187.0, b53 = 64448.0 / 6561.0,
102 b54 = -212.0 / 729.0,
103
104 b61 = 9017.0 / 3168.0 , b62 = -355.0 / 33.0,
105 b63 = 46732.0 / 5247.0, b64 = 49.0 / 176.0,
106 b65 = -5103.0 / 18656.0,
107
108 b71 = 35.0 / 384.0, b72 = 0.,
109 b73 = 500.0 / 1113.0, b74 = 125.0 / 192.0,
110 b75 = -2187.0 / 6784.0, b76 = 11.0 / 84.0,
111
112 //Sum of columns, sum(bij) = ei
113 // e1 = 0. ,
114 // e2 = 1.0/5.0 ,
115 // e3 = 3.0/10.0 ,
116 // e4 = 4.0/5.0 ,
117 // e5 = 8.0/9.0 ,
118 // e6 = 1.0 ,
119 // e7 = 1.0 ,
120
121 // Difference between the higher and the lower order method coeff. :
122 // b7j are the coefficients of higher order
123
124 dc1 = -(b71 - 5179.0 / 57600.0),
125 dc2 = -(b72 - .0),
126 dc3 = -(b73 - 7571.0 / 16695.0),
127 dc4 = -(b74 - 393.0 / 640.0),
128 dc5 = -(b75 + 92097.0 / 339200.0),
129 dc6 = -(b76 - 187.0 / 2100.0),
130 dc7 = -(- 1.0 / 40.0);
131
132 const G4int numberOfVariables = GetNumberOfVariables();
133 State yTemp = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
134
135 // The number of variables to be integrated over
136 //
137 yOut[7] = yTemp[7] = yInput[7];
138
139 // Saving yInput because yInput and yOut can be aliases for same array
140 //
141 for(G4int i = 0; i < numberOfVariables; ++i)
142 {
143 fyIn[i] = yInput[i];
144 }
145 // RightHandSide(yIn, dydx); // Not done! 1st stage
146
147 for(G4int i = 0; i < numberOfVariables; ++i)
148 {
149 yTemp[i] = fyIn[i] + b21 * hstep * dydx[i];
150 }
151 RightHandSide(yTemp, ak2); // 2nd stage
152
153 for(G4int i = 0; i < numberOfVariables; ++i)
154 {
155 yTemp[i] = fyIn[i] + hstep * (b31 * dydx[i] + b32 * ak2[i]);
156 }
157 RightHandSide(yTemp, ak3); // 3rd stage
158
159 for(G4int i = 0; i < numberOfVariables; ++i)
160 {
161 yTemp[i] = fyIn[i] + hstep * (
162 b41 * dydx[i] + b42 * ak2[i] + b43 * ak3[i]);
163 }
164 RightHandSide(yTemp, ak4); // 4th stage
165
166 for(G4int i = 0; i < numberOfVariables; ++i)
167 {
168 yTemp[i] = fyIn[i] + hstep * (
169 b51 * dydx[i] + b52 * ak2[i] + b53 * ak3[i] + b54 * ak4[i]);
170 }
171 RightHandSide(yTemp, ak5); // 5th stage
172
173 for(G4int i = 0; i < numberOfVariables; ++i)
174 {
175 yTemp[i] = fyIn[i] + hstep * (
176 b61 * dydx[i] + b62 * ak2[i] +
177 b63 * ak3[i] + b64 * ak4[i] + b65 * ak5[i]);
178 }
179 RightHandSide(yTemp, ak6); // 6th stage
180
181 for(G4int i = 0; i < numberOfVariables; ++i)
182 {
183 yOut[i] = fyIn[i] + hstep * (
184 b71 * dydx[i] + b72 * ak2[i] + b73 * ak3[i] +
185 b74 * ak4[i] + b75 * ak5[i] + b76 * ak6[i]);
186 }
187 RightHandSide(yOut, ak7); // 7th and Final stage
188
189 for(G4int i = 0; i < numberOfVariables; ++i)
190 {
191 yErr[i] = hstep * (
192 dc1 * dydx[i] + dc2 * ak2[i] +
193 dc3 * ak3[i] + dc4 * ak4[i] +
194 dc5 * ak5[i] + dc6 * ak6[i] + dc7 * ak7[i]
195 ) + 1.5e-18;
196
197 // Store Input and Final values, for possible use in calculating chord
198 //
199 fyOut[i] = yOut[i];
200 fdydxIn[i] = dydx[i];
201 }
202
203 fLastStepLength = hstep;
204}
205
207{
208 // Coefficients were taken from Some Practical Runge-Kutta Formulas
209 // by Lawrence F. Shampine, page 149, c*
210 //
211 const G4double hf1 = 6025192743.0 / 30085553152.0,
212 hf3 = 51252292925.0 / 65400821598.0,
213 hf4 = - 2691868925.0 / 45128329728.0,
214 hf5 = 187940372067.0 / 1594534317056.0,
215 hf6 = - 1776094331.0 / 19743644256.0,
216 hf7 = 11237099.0 / 235043384.0;
217
218 G4ThreeVector mid;
219
220 for(G4int i = 0; i < 3; ++i)
221 {
222 mid[i] = fyIn[i] + 0.5 * fLastStepLength * (
223 hf1 * fdydxIn[i] + hf3 * ak3[i] +
224 hf4 * ak4[i] + hf5 * ak5[i] + hf6 * ak6[i] + hf7 * ak7[i]);
225 }
226
227 const G4ThreeVector begin = makeVector(fyIn, Value3D::Position);
228 const G4ThreeVector end = makeVector(fyOut, Value3D::Position);
229
230 return G4LineSection::Distline(mid, begin, end);
231}
232
233// The lower (4th) order interpolant given by Dormand and Prince:
234// J. R. Dormand and P. J. Prince, "Runge-Kutta triples"
235// Computers & Mathematics with Applications, vol. 12, no. 9,
236// pp. 1007-1017, 1986.
237//
239Interpolate4thOrder(G4double yOut[], G4double tau) const
240{
241 const G4int numberOfVariables = GetNumberOfVariables();
242
243 const G4double tau2 = tau * tau,
244 tau3 = tau * tau2,
245 tau4 = tau2 * tau2;
246
247 const G4double bf1 = 1.0 / 11282082432.0 * (
248 157015080.0 * tau4 - 13107642775.0 * tau3 + 34969693132.0 * tau2 -
249 32272833064.0 * tau + 11282082432.0);
250
251 const G4double bf3 = - 100.0 / 32700410799.0 * tau * (
252 15701508.0 * tau3 - 914128567.0 * tau2 + 2074956840.0 * tau -
253 1323431896.0);
254
255 const G4double bf4 = 25.0 / 5641041216.0 * tau * (
256 94209048.0 * tau3 - 1518414297.0 * tau2 + 2460397220.0 * tau -
257 889289856.0);
258
259 const G4double bf5 = - 2187.0 / 199316789632.0 * tau * (
260 52338360.0 * tau3 - 451824525.0 * tau2 + 687873124.0 * tau -
261 259006536.0);
262
263 const G4double bf6 = 11.0 / 2467955532.0 * tau * (
264 106151040.0 * tau3 - 661884105.0 * tau2 +
265 946554244.0 * tau - 361440756.0);
266
267 const G4double bf7 = 1.0 / 29380423.0 * tau * (1.0 - tau) * (
268 8293050.0 * tau2 - 82437520.0 * tau + 44764047.0);
269
270 for(G4int i = 0; i < numberOfVariables; ++i)
271 {
272 yOut[i] = fyIn[i] + fLastStepLength * tau * (
273 bf1 * fdydxIn[i] + bf3 * ak3[i] + bf4 * ak4[i] +
274 bf5 * ak5[i] + bf6 * ak6[i] + bf7 * ak7[i]);
275 }
276}
277
278// Following interpolant of order 5 was given by Baker,Dormand,Gilmore, Prince :
279// T. S. Baker, J. R. Dormand, J. P. Gilmore, and P. J. Prince,
280// "Continuous approximation with embedded Runge-Kutta methods"
281// Applied Numerical Mathematics, vol. 22, no. 1, pp. 51-62, 1996.
282//
283// Calculating the extra stages for the interpolant
284//
286{
287 // Coefficients for the additional stages
288 //
289 const G4double b81 = 6245.0 / 62208.0,
290 b82 = 0.0,
291 b83 = 8875.0 / 103032.0,
292 b84 = -125.0 / 1728.0,
293 b85 = 801.0 / 13568.0,
294 b86 = -13519.0 / 368064.0,
295 b87 = 11105.0 / 368064.0,
296
297 b91 = 632855.0 / 4478976.0,
298 b92 = 0.0,
299 b93 = 4146875.0 / 6491016.0,
300 b94 = 5490625.0 /14183424.0,
301 b95 = -15975.0 / 108544.0,
302 b96 = 8295925.0 / 220286304.0,
303 b97 = -1779595.0 / 62938944.0,
304 b98 = -805.0 / 4104.0;
305
306 const G4int numberOfVariables = GetNumberOfVariables();
307 State yTemp = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
308
309 // Evaluate the extra stages
310 //
311 for(G4int i = 0; i < numberOfVariables; ++i)
312 {
313 yTemp[i] = fyIn[i] + fLastStepLength * (
314 b81 * fdydxIn[i] + b82 * ak2[i] + b83 * ak3[i] +
315 b84 * ak4[i] + b85 * ak5[i] + b86 * ak6[i] +
316 b87 * ak7[i]
317 );
318 }
319 RightHandSide(yTemp, ak8); // 8th Stage
320
321 for(G4int i = 0; i < numberOfVariables; ++i)
322 {
323 yTemp[i] = fyIn[i] + fLastStepLength * (
324 b91 * fdydxIn[i] + b92 * ak2[i] + b93 * ak3[i] +
325 b94 * ak4[i] + b95 * ak5[i] + b96 * ak6[i] +
326 b97 * ak7[i] + b98 * ak8[i]
327 );
328 }
329 RightHandSide(yTemp, ak9); // 9th Stage
330}
331
332// Calculating the interpolated result yOut with the coefficients
333//
335Interpolate5thOrder(G4double yOut[], G4double tau) const
336{
337 // Define the coefficients for the polynomials
338 //
339 G4double bi[10][5];
340
341 // COEFFICIENTS OF bi[1]
342 bi[1][0] = 1.0,
343 bi[1][1] = -38039.0 / 7040.0,
344 bi[1][2] = 125923.0 / 10560.0,
345 bi[1][3] = -19683.0 / 1760.0,
346 bi[1][4] = 3303.0 / 880.0,
347 // --------------------------------------------------------
348 //
349 // COEFFICIENTS OF bi[2]
350 bi[2][0] = 0.0,
351 bi[2][1] = 0.0,
352 bi[2][2] = 0.0,
353 bi[2][3] = 0.0,
354 bi[2][4] = 0.0,
355 // --------------------------------------------------------
356 //
357 // COEFFICIENTS OF bi[3]
358 bi[3][0] = 0.0,
359 bi[3][1] = -12500.0 / 4081.0,
360 bi[3][2] = 205000.0 / 12243.0,
361 bi[3][3] = -90000.0 / 4081.0,
362 bi[3][4] = 36000.0 / 4081.0,
363 // --------------------------------------------------------
364 //
365 // COEFFICIENTS OF bi[4]
366 bi[4][0] = 0.0,
367 bi[4][1] = -3125.0 / 704.0,
368 bi[4][2] = 25625.0 / 1056.0,
369 bi[4][3] = -5625.0 / 176.0,
370 bi[4][4] = 1125.0 / 88.0,
371 // --------------------------------------------------------
372 //
373 // COEFFICIENTS OF bi[5]
374 bi[5][0] = 0.0,
375 bi[5][1] = 164025.0 / 74624.0,
376 bi[5][2] = -448335.0 / 37312.0,
377 bi[5][3] = 295245.0 / 18656.0,
378 bi[5][4] = -59049.0 / 9328.0,
379 // --------------------------------------------------------
380 //
381 // COEFFICIENTS OF bi[6]
382 bi[6][0] = 0.0,
383 bi[6][1] = -25.0 / 28.0,
384 bi[6][2] = 205.0 / 42.0,
385 bi[6][3] = -45.0 / 7.0,
386 bi[6][4] = 18.0 / 7.0,
387 // --------------------------------------------------------
388 //
389 // COEFFICIENTS OF bi[7]
390 bi[7][0] = 0.0,
391 bi[7][1] = -2.0 / 11.0,
392 bi[7][2] = 73.0 / 55.0,
393 bi[7][3] = -171.0 / 55.0,
394 bi[7][4] = 108.0 / 55.0,
395 // --------------------------------------------------------
396 //
397 // COEFFICIENTS OF bi[8]
398 bi[8][0] = 0.0,
399 bi[8][1] = 189.0 / 22.0,
400 bi[8][2] = -1593.0 / 55.0,
401 bi[8][3] = 3537.0 / 110.0,
402 bi[8][4] = -648.0 / 55.0,
403 // --------------------------------------------------------
404 //
405 // COEFFICIENTS OF bi[9]
406 bi[9][0] = 0.0,
407 bi[9][1] = 351.0 / 110.0,
408 bi[9][2] = -999.0 / 55.0,
409 bi[9][3] = 2943.0 / 110.0,
410 bi[9][4] = -648.0 / 55.0;
411 // --------------------------------------------------------
412
413 // Calculating the polynomials
414
415 G4double b[10];
416 std::memset(b, 0.0, sizeof(b));
417
418 G4double tauPower = 1.0;
419 for(G4int j = 0; j <= 4; ++j)
420 {
421 for(G4int iStage = 1; iStage <= 9; ++iStage)
422 {
423 b[iStage] += bi[iStage][j] * tauPower;
424 }
425 tauPower *= tau;
426 }
427
428 const G4int numberOfVariables = GetNumberOfVariables();
429 const G4double stepLen = fLastStepLength * tau;
430 for(G4int i = 0; i < numberOfVariables; ++i)
431 {
432 yOut[i] = fyIn[i] + stepLen * (
433 b[1] * fdydxIn[i] + b[2] * ak2[i] + b[3] * ak3[i] +
434 b[4] * ak4[i] + b[5] * ak5[i] + b[6] * ak6[i] +
435 b[7] * ak7[i] + b[8] * ak8[i] + b[9] * ak9[i]
436 );
437 }
438}
double G4double
Definition G4Types.hh:83
int G4int
Definition G4Types.hh:85
void Stepper(const G4double yInput[], const G4double dydx[], G4double hstep, G4double yOutput[], G4double yError[]) override
G4DormandPrince745(G4EquationOfMotion *equation, G4int numberOfVariables=6)
void Interpolate4thOrder(G4double yOut[], G4double tau) const
G4double DistChord() const override
void Interpolate5thOrder(G4double yOut[], G4double tau) const
static G4double Distline(const G4ThreeVector &OtherPnt, const G4ThreeVector &LinePntA, const G4ThreeVector &LinePntB)
G4int GetNumberOfVariables() const
void RightHandSide(const G4double y[], G4double dydx[]) const
G4ThreeVector makeVector(const ArrayType &array, Value3D value)
G4double[G4FieldTrack::ncompSVEC] State
void copy(G4double dst[], const G4double src[], std::size_t size=G4FieldTrack::ncompSVEC)