Geant4 Cross Reference |
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 58 using namespace field_utils; 59 60 // Name of this steppers 61 const G4String& G4DormandPrince745::StepperType() const 62 { 63 static G4String _stepperType("G4DormandPrince745: 5th order"); 64 return _stepperType; 65 } 66 67 // Description of this steppers - plus details of its implementation 68 const G4String& G4DormandPrince745::StepperDescription() const 69 { 70 static G4String _stepperDescription( 71 "Embedeed 5th order Runge-Kutta stepper - 7 stages, FSAL, Interpolating."); 72 return _stepperDescription; 73 } 74 75 G4DormandPrince745::G4DormandPrince745(G4EquationOfMotion* equation, 76 G4int noIntegrationVariables) 77 : G4MagIntegratorStepper(equation, noIntegrationVariables) 78 { 79 } 80 81 void G4DormandPrince745::Stepper(const G4double yInput[], 82 const G4double dydx[], 83 G4double hstep, 84 G4double yOutput[], 85 G4double yError[], 86 G4double dydxOutput[]) 87 { 88 Stepper(yInput, dydx, hstep, yOutput, yError); 89 field_utils::copy(dydxOutput, ak7); 90 } 91 92 // Stepper 93 // 94 // Passing in the value of yInput[],the first time dydx[] and Step length 95 // Giving back yOut and yErr arrays for output and error respectively 96 // 97 void G4DormandPrince745::Stepper(const G4double yInput[], 98 const G4double dydx[], 99 G4double hstep, 100 G4double yOut[], 101 G4double yErr[]) 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 } 213 214 G4double G4DormandPrince745::DistChord() const 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 } 240 241 // The lower (4th) order interpolant given by Dormand and Prince: 242 // J. R. Dormand and P. J. Prince, "Runge-Kutta triples" 243 // Computers & Mathematics with Applications, vol. 12, no. 9, 244 // pp. 1007-1017, 1986. 245 // 246 void G4DormandPrince745:: 247 Interpolate4thOrder(G4double yOut[], G4double tau) const 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 } 285 286 // Following interpolant of order 5 was given by Baker,Dormand,Gilmore, Prince : 287 // T. S. Baker, J. R. Dormand, J. P. Gilmore, and P. J. Prince, 288 // "Continuous approximation with embedded Runge-Kutta methods" 289 // Applied Numerical Mathematics, vol. 22, no. 1, pp. 51-62, 1996. 290 // 291 // Calculating the extra stages for the interpolant 292 // 293 void G4DormandPrince745::SetupInterpolation5thOrder() 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 } 339 340 // Calculating the interpolated result yOut with the coefficients 341 // 342 void G4DormandPrince745:: 343 Interpolate5thOrder(G4double yOut[], G4double tau) const 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 } 447