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 // G4FSALBogackiShampine45 implementation 27 // 28 // The Butcher table of the Bogacki-Shampine-8 29 // 30 // 0 | 31 // 1/6 | 1/6 32 // 2/9 | 2/27 4/27 33 // 3/7 | 183/1372 -162/343 1053/1372 34 // 2/3 | 68/297 -4/11 42/143 196 35 // 3/4 | 597/22528 81/352 63099/5857 36 // 1 | 174197/959244 -30942/79937 8152137/19 37 // 1 | 587/8064 0 4440339/15 38 // ------------------------------------------- 39 // 587/8064 0 4440339/15 40 // 2479/34992 0 123/416 41 // 42 // Created: Somnath Banerjee, Google Summer of 43 // Supervision: John Apostolakis, CERN 44 // ------------------------------------------- 45 46 // Plan is that this source file / class will 47 // BogackiShampine45 class, which contains imp 48 49 #include <cassert> 50 51 #include "G4FSALBogackiShampine45.hh" 52 #include "G4LineSection.hh" 53 54 G4bool G4FSALBogackiShampine45::fPreparedCon 55 G4double G4FSALBogackiShampine45::bi[12][7]; 56 57 // Constructor 58 // 59 G4FSALBogackiShampine45::G4FSALBogackiShampine 60 61 62 : G4VFSALIntegrationStepper(EqRhs, noIntegr 63 { 64 const G4int numberOfVariables = noIntegrat 65 66 // New Chunk of memory being created for u 67 68 // aki - for storing intermediate RHS 69 // 70 ak2 = new G4double[numberOfVariables]; 71 ak3 = new G4double[numberOfVariables]; 72 ak4 = new G4double[numberOfVariables]; 73 ak5 = new G4double[numberOfVariables]; 74 ak6 = new G4double[numberOfVariables]; 75 ak7 = new G4double[numberOfVariables]; 76 ak8 = new G4double[numberOfVariables]; 77 78 ak9 = new G4double[numberOfVariables]; 79 ak10 = new G4double[numberOfVariables]; 80 ak11 = new G4double[numberOfVariables]; 81 DyDx = new G4double[numberOfVariables]; 82 83 assert ( GetNumberOfStateVariables() >= 8 84 const G4int numStateVars = std::max(noInte 85 GetNum 86 87 // Must ensure space extra 'state' variabl 88 // 89 yTemp = new G4double[numStateVars]; 90 yIn = new G4double[numStateVars] ; 91 92 fLastInitialVector = new G4double[numState 93 fLastFinalVector = new G4double[numStateVa 94 fLastDyDx = new G4double[numberOfVariables 95 96 fMidVector = new G4double[numStateVars]; 97 fMidError = new G4double[numStateVars]; 98 99 pseudoDydx_for_DistChord = new G4double[nu 100 101 fMidVector = new G4double[numberOfVariable 102 fMidError = new G4double[numberOfVariable 103 if( primary ) 104 { 105 fAuxStepper = new G4FSALBogackiShampine4 106 107 } 108 if( !fPreparedConstants ) 109 { 110 PrepareConstants(); 111 } 112 } 113 114 // Destructor 115 // 116 G4FSALBogackiShampine45::~G4FSALBogackiShampin 117 { 118 // Clear all previously allocated memory f 119 120 delete [] ak2; 121 delete [] ak3; 122 delete [] ak4; 123 delete [] ak5; 124 delete [] ak6; 125 delete [] ak7; 126 delete [] ak8; 127 delete [] ak9; 128 delete [] ak10; 129 delete [] ak11; 130 delete [] DyDx; 131 delete [] yTemp; 132 delete [] yIn; 133 134 delete [] fLastInitialVector; 135 delete [] fLastFinalVector; 136 delete [] fLastDyDx; 137 delete [] fMidVector; 138 delete [] fMidError; 139 140 delete fAuxStepper; 141 142 delete [] pseudoDydx_for_DistChord; 143 } 144 145 // Stepper 146 // 147 // Passing in the value of yInput[],the first 148 // Giving back yOut and yErr arrays for output 149 // 150 void G4FSALBogackiShampine45::Stepper(const G4 151 const G4 152 G4 153 G4 154 G4 155 G4 156 { 157 G4int i; 158 159 // The various constants defined on the ba 160 161 const G4double b21 = 1.0/6.0 , 162 b31 = 2.0/27.0 , b32 = 4.0/ 163 164 b41 = 183.0/1372.0 , b42 = 165 166 b51 = 68.0/297.0, b52 = -4. 167 b53 = 42.0/143.0, b54 = 196 168 169 b61 = 597.0/22528.0, b62 = 170 b63 = 63099.0/585728.0, b64 171 b65 = 4617.0/20480.0, 172 173 b71 = 174197.0/959244.0, b7 174 b73 = 8152137.0/19744439.0, 175 b75 = -29421.0/29068.0, b7 176 177 b81 = 587.0/8064.0, b82 = 178 b83 = 4440339.0/15491840.0, 179 b85 = 387.0/44800.0, b86 = 180 b87 = 7267.0/94080.0, 181 182 183 // c1 = 2479.0/34992.0, 184 // c2 = 0.0, 185 // c3 = 123.0/416.0, 186 // c4 = 612941.0/3411720.0, 187 // c5 = 43.0/1440.0, 188 // c6 = 2272.0/6561.0, 189 // c7 = 79937.0/1113912.0, 190 // c8 = 3293.0/556956.0, 191 192 // For the embedded higher order method on 193 // taken and is used directly later instea 194 // of butcher table in a separate set of v 195 // difference there 196 197 dc1 = b81 - 2479.0/34992.0 198 dc2 = 0.0, 199 dc3 = b83 - 123.0/416.0 , 200 dc4 = b84 - 612941.0/341172 201 dc5 = b85 - 43.0/1440.0, 202 dc6 = b86 - 2272.0/6561.0, 203 dc7 = b87 - 79937.0/1113912 204 dc8 = -3293.0/556956.0; / 205 206 const G4int numberOfVariables = GetNumberO 207 208 // The number of variables to be integrate 209 // 210 yOut[7] = yTemp[7] = yIn[7]; 211 212 // Saving yInput because yInput and yOut 213 // 214 for(i=0; i<numberOfVariables; ++i) 215 { 216 yIn[i]=yInput[i]; 217 DyDx[i] = dydx[i]; 218 } 219 // RightHandSide(yIn, dydx) ; // 1st Ste 220 221 for(i=0; i<numberOfVariables; ++i) 222 { 223 yTemp[i] = yIn[i] + b21*Step*DyDx[i] ; 224 } 225 RightHandSide(yTemp, ak2) ; / 226 227 for(i=0; i<numberOfVariables; ++i) 228 { 229 yTemp[i] = yIn[i] + Step*(b31*DyDx[i] 230 } 231 RightHandSide(yTemp, ak3) ; / 232 233 for(i=0; i<numberOfVariables; ++i) 234 { 235 yTemp[i] = yIn[i] + Step*(b41*DyDx[i] 236 } 237 RightHandSide(yTemp, ak4) ; / 238 239 for(i=0; i<numberOfVariables; ++i) 240 { 241 yTemp[i] = yIn[i] + Step*(b51*DyDx[i] 242 b54*ak4[i]) 243 } 244 RightHandSide(yTemp, ak5) ; / 245 246 for(i=0; i<numberOfVariables; ++i) 247 { 248 yTemp[i] = yIn[i] + Step*(b61*DyDx[i] 249 b64*ak4[i] + 250 } 251 RightHandSide(yTemp, ak6) ; / 252 253 for(i=0; i<numberOfVariables; ++i) 254 { 255 yTemp[i] = yIn[i] + Step*(b71*DyDx[i] 256 b74*ak4[i] + 257 } 258 RightHandSide(yTemp, ak7); / 259 260 for(i=0; i<numberOfVariables; ++i) 261 { 262 yOut[i] = yIn[i] + Step*(b81*DyDx[i] + 263 b84*ak4[i] + 264 b87*ak7[i]); 265 } 266 RightHandSide(yOut, ak8); / 267 268 269 for(i=0; i<numberOfVariables; ++i) 270 { 271 272 yErr[i] = Step*(dc1*DyDx[i] + dc2*ak2[ 273 dc5*ak5[i] + dc6*ak6[i 274 275 276 // FSAL stepper : Must pass the last D 277 // 278 nextDydx[i] = ak8[i]; 279 280 // Store Input and Final values, for p 281 // 282 fLastInitialVector[i] = yIn[i] ; 283 fLastFinalVector[i] = yOut[i]; 284 fLastDyDx[i] = DyDx[i]; 285 } 286 fLastStepLength = Step; 287 288 return; 289 } 290 291 // DistChord 292 // 293 G4double G4FSALBogackiShampine45::DistChord() 294 { 295 G4double distLine, distChord; 296 G4ThreeVector initialPoint, finalPoint, mi 297 298 // Store last initial and final points 299 // (they will be overwritten in self-Stepp 300 // 301 initialPoint = G4ThreeVector( fLastInitial 302 fLastInitialV 303 finalPoint = G4ThreeVector( fLastFinalVe 304 fLastFinalVec 305 306 // Do half a step using StepNoErr 307 308 fAuxStepper->Stepper( fLastInitialVector, 309 fMidVector, fMidErro 310 311 midPoint = G4ThreeVector( fMidVector[0], f 312 313 // Use stored values of Initial and Endpoi 314 // distance of Chord 315 // 316 if (initialPoint != finalPoint) 317 { 318 distLine = G4LineSection::Distline(mid 319 distChord = distLine; 320 } 321 else 322 { 323 distChord = (midPoint-initialPoint).ma 324 } 325 return distChord; 326 } 327 328 // PrepareConstants 329 // 330 void G4FSALBogackiShampine45::PrepareConstants 331 { 332 // -------------------------------------- 333 // COEFFICIENTS FOR INTERPOLANT bi WITH 334 // -------------------------------------- 335 336 // Initialise all values of G4double bi[12 337 // 338 for(auto i=1; i<12; ++i) 339 { 340 for(auto j=1; j<7; ++j) 341 { 342 bi[i][j] = 0.0 ; 343 } 344 } 345 346 bi[1][6] = -12134338393.0/1050809760.0 , 347 bi[1][5] = -1620741229.0/50038560.0 , 348 bi[1][4] = -2048058893.0/59875200.0 , 349 bi[1][3] = -87098480009.0/5254048800.0 , 350 bi[1][2] = -11513270273.0/3502699200.0 , 351 // 352 bi[3][6] = -33197340367.0/1218433216.0 , 353 bi[3][5] = -539868024987.0/6092166080.0 , 354 bi[3][4] = -39991188681.0/374902528.0 , 355 bi[3][3] = -69509738227.0/1218433216.0 , 356 bi[3][2] = -29327744613.0/2436866432.0 , 357 // 358 bi[4][6] = -284800997201.0/19905339168.0 , 359 bi[4][5] = -7896875450471.0/165877826400.0 360 bi[4][4] = -333945812879.0/5671036800.0 , 361 bi[4][3] = -16209923456237.0/497633479200. 362 bi[4][2] = -2382590741699.0/331755652800.0 363 // 364 bi[5][6] = -540919.0/741312.0 , 365 bi[5][5] = -103626067.0/43243200.0 , 366 bi[5][4] = -633779.0/211200.0 , 367 bi[5][3] = -32406787.0/18532800.0 , 368 bi[5][2] = -36591193.0/86486400.0 , 369 // 370 bi[6][6] = 7157998304.0/374350977.0 , 371 bi[6][5] = 30405842464.0/623918295.0 , 372 bi[6][4] = 183022264.0/5332635.0 , 373 bi[6][3] = -3357024032.0/1871754885.0 , 374 bi[6][2] = -611586736.0/89131185.0 , 375 // 376 bi[7][6] = -138073.0/9408.0 , 377 bi[7][5] = -719433.0/15680.0 , 378 bi[7][4] = -1620541.0/31360.0 , 379 bi[7][3] = -385151.0/15680.0 , 380 bi[7][2] = -65403.0/15680.0 , 381 // 382 bi[8][6] = 1245.0/64.0 , 383 bi[8][5] = 3991.0/64.0 , 384 bi[8][4] = 4715.0/64.0 , 385 bi[8][3] = 2501.0/64.0 , 386 bi[8][2] = 149.0/16.0 , 387 bi[8][1] = 1.0 , 388 // 389 bi[9][6] = 55.0/3.0 , 390 bi[9][5] = 71.0 , 391 bi[9][4] = 103.0 , 392 bi[9][3] = 199.0/3.0 , 393 bi[9][2] = 16.0 , 394 // 395 bi[10][6] = -1774004627.0/75810735.0 , 396 bi[10][5] = -1774004627.0/25270245.0 , 397 bi[10][4] = -26477681.0/359975.0 , 398 bi[10][3] = -11411880511.0/379053675.0 , 399 bi[10][2] = -423642896.0/126351225.0 , 400 // 401 bi[11][6] = 35.0 , 402 bi[11][5] = 105.0 , 403 bi[11][4] = 117.0 , 404 bi[11][3] = 59.0 , 405 bi[11][2] = 12.0 ; 406 } 407 408 // ------------------------------------------- 409 410 void G4FSALBogackiShampine45::interpolate( con 411 con 412 413 414 415 { 416 const G4double a91 = 455.0/6144.0 , 417 a92 = 0.0 , 418 a93 = 10256301.0/35409920.0 419 a94 = 2307361.0/17971200.0 420 a95 = -387.0/102400.0 , 421 a96 = 73.0/5130.0 , 422 a97 = -7267.0/215040.0 , 423 a98 = 1.0/32.0 , 424 425 a101 = -837888343715.0/1317 426 a102 = 30409415.0/52955362. 427 a103 = -48321525963.0/75916 428 a104 = 8530738453321.0/1976 429 a105 = 1361640523001.0/1626 430 a106 = -13143060689.0/38604 431 a107 = 18700221969.0/379584 432 a108 = -5831595.0/847285792 433 a109 = -5183640.0/26477681. 434 435 a111 = 98719073263.0/155196 436 a112 = 1307.0/123552.0 , 437 a113 = 4632066559387.0/7018 438 a114 = 7828594302389.0/3821 439 a115 = 40763687.0/110702592 440 a116 = 34872732407.0/224610 441 a117 = -2561897.0/30105600. 442 a118 = 1.0/10.0 , 443 a119 = -1.0/10.0 , 444 a1110 = -1403317093.0/11371 445 446 const G4int numberOfVariables = GetNumberO 447 448 // Saving yInput because yInput and yOut c 449 // 450 for(G4int i=0; i<numberOfVariables; ++i) 451 { 452 yIn[i]=yInput[i]; 453 } 454 455 // The number of variables to be integrate 456 // 457 yOut[7] = yTemp[7] = yIn[7]; 458 459 // Calculating extra stages 460 // 461 for(G4int i=0; i<numberOfVariables; ++i) 462 { 463 yTemp[i] = yIn[i] + Step*(a91*dydx[i] 464 a94*ak4[i] + 465 a97*ak7[i] + 466 } 467 468 RightHandSide(yTemp, ak9); 469 470 for(G4int i=0; i<numberOfVariables; ++i) 471 { 472 yTemp[i] = yIn[i] + Step*(a101*dydx[i] 473 a104*ak4[i] 474 a107*ak7[i] 475 } 476 477 RightHandSide(yTemp, ak10); 478 479 for(G4int i=0; i<numberOfVariables; ++i) 480 { 481 yTemp[i] = yIn[i] + Step*(a111*dydx[i] 482 a114*ak4[i] 483 a117*ak7[i] 484 a1110*ak10[i 485 } 486 487 RightHandSide(yTemp, ak11); 488 489 G4double tau0 = tau; 490 491 // Calculating the polynomials 492 // 493 for(auto i=1; i<=11; ++i) // i is NOT the 494 { 495 b[i] = 0.0; 496 tau = tau0; 497 for(auto j=1; j<=6; ++j) 498 { 499 b[i] += bi[i][j]*tau; 500 tau*=tau0; 501 } 502 } 503 504 for(G4int i=0; i<numberOfVariables; ++i) 505 { 506 yOut[i] = yIn[i] + Step*(b[1]*dydx[i] 507 b[4]*ak4[i] + 508 b[7]*ak7[i] + 509 b[10]*ak10[i] 510 } 511 } 512