Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** 25 // 25 // 26 // G4DoLoMcPriRK34 implementation 26 // G4DoLoMcPriRK34 implementation 27 // 27 // 28 // Created: Somnath Banerjee, Google Summer of 28 // Created: Somnath Banerjee, Google Summer of Code 2015, 7 July 2015 29 // Supervision: John Apostolakis, CERN 29 // Supervision: John Apostolakis, CERN 30 // ------------------------------------------- 30 // -------------------------------------------------------------------- 31 31 32 #include "G4DoLoMcPriRK34.hh" 32 #include "G4DoLoMcPriRK34.hh" 33 #include "G4LineSection.hh" 33 #include "G4LineSection.hh" 34 34 35 // Constructor 35 // Constructor 36 // 36 // 37 G4DoLoMcPriRK34::G4DoLoMcPriRK34(G4EquationOfM 37 G4DoLoMcPriRK34::G4DoLoMcPriRK34(G4EquationOfMotion* EqRhs, 38 G4int noInte 38 G4int noIntegrationVariables, 39 G4bool prima 39 G4bool primary) 40 : G4MagIntegratorStepper(EqRhs, noIntegrati 40 : G4MagIntegratorStepper(EqRhs, noIntegrationVariables) 41 { 41 { 42 const G4int numberOfVariables = noIntegrati 42 const G4int numberOfVariables = noIntegrationVariables; 43 43 44 // New Chunk of memory being created for us 44 // New Chunk of memory being created for use by the Stepper 45 45 46 // aki - for storing intermediate RHS 46 // aki - for storing intermediate RHS 47 // 47 // 48 ak2 = new G4double[numberOfVariables]; 48 ak2 = new G4double[numberOfVariables]; 49 ak3 = new G4double[numberOfVariables]; 49 ak3 = new G4double[numberOfVariables]; 50 ak4 = new G4double[numberOfVariables]; 50 ak4 = new G4double[numberOfVariables]; 51 ak5 = new G4double[numberOfVariables]; 51 ak5 = new G4double[numberOfVariables]; 52 ak6 = new G4double[numberOfVariables]; 52 ak6 = new G4double[numberOfVariables]; 53 53 54 yTemp = new G4double[numberOfVariables] ; 54 yTemp = new G4double[numberOfVariables] ; 55 yIn = new G4double[numberOfVariables] ; 55 yIn = new G4double[numberOfVariables] ; 56 56 57 fLastInitialVector = new G4double[numberOfV 57 fLastInitialVector = new G4double[numberOfVariables] ; 58 fLastFinalVector = new G4double[numberOfVar 58 fLastFinalVector = new G4double[numberOfVariables] ; 59 fLastDyDx = new G4double[numberOfVariables] 59 fLastDyDx = new G4double[numberOfVariables]; 60 60 61 fMidVector = new G4double[numberOfVariables 61 fMidVector = new G4double[numberOfVariables]; 62 fMidError = new G4double[numberOfVariables] 62 fMidError = new G4double[numberOfVariables]; 63 if( primary ) 63 if( primary ) 64 { 64 { 65 fAuxStepper = new G4DoLoMcPriRK34(EqRhs, 65 fAuxStepper = new G4DoLoMcPriRK34(EqRhs, numberOfVariables, !primary); 66 } 66 } 67 } 67 } 68 68 69 // Destructor 69 // Destructor 70 // 70 // 71 G4DoLoMcPriRK34::~G4DoLoMcPriRK34() 71 G4DoLoMcPriRK34::~G4DoLoMcPriRK34() 72 { 72 { 73 // clear all previously allocated memory fo 73 // clear all previously allocated memory for Stepper and DistChord 74 74 75 delete [] ak2; 75 delete [] ak2; 76 delete [] ak3; 76 delete [] ak3; 77 delete [] ak4; 77 delete [] ak4; 78 delete [] ak5; 78 delete [] ak5; 79 delete [] ak6; 79 delete [] ak6; 80 80 81 delete [] yTemp; 81 delete [] yTemp; 82 delete [] yIn; 82 delete [] yIn; 83 83 84 delete [] fLastInitialVector; 84 delete [] fLastInitialVector; 85 delete [] fLastFinalVector; 85 delete [] fLastFinalVector; 86 delete [] fLastDyDx; 86 delete [] fLastDyDx; 87 delete [] fMidVector; 87 delete [] fMidVector; 88 delete [] fMidError; 88 delete [] fMidError; 89 89 90 delete fAuxStepper; 90 delete fAuxStepper; 91 } 91 } 92 92 93 // Stepper 93 // Stepper 94 // 94 // 95 // Passing in the value of yInput[],the first 95 // Passing in the value of yInput[],the first time dydx[] and Step length 96 // Giving back yOut and yErr arrays for output 96 // Giving back yOut and yErr arrays for output and error respectively 97 // 97 // 98 void G4DoLoMcPriRK34::Stepper(const G4double y 98 void G4DoLoMcPriRK34::Stepper(const G4double yInput[], 99 const G4double D 99 const G4double DyDx[], 100 G4double S 100 G4double Step, 101 G4double y 101 G4double yOut[], 102 G4double y 102 G4double yErr[] ) 103 { 103 { 104 G4int i; 104 G4int i; 105 105 106 // The various constants defined on the bas 106 // The various constants defined on the basis of butcher tableu 107 // 107 // 108 const G4double b21 = 7.0/27.0 , 108 const G4double b21 = 7.0/27.0 , 109 b31 = 7.0/72.0 , 109 b31 = 7.0/72.0 , 110 b32 = 7.0/24.0 , 110 b32 = 7.0/24.0 , 111 111 112 b41 = 3043.0/3528.0 , 112 b41 = 3043.0/3528.0 , 113 b42 = -3757.0/1176.0 , 113 b42 = -3757.0/1176.0 , 114 b43 = 1445.0/441.0, 114 b43 = 1445.0/441.0, 115 115 116 b51 = 17617.0/11662.0 , 116 b51 = 17617.0/11662.0 , 117 b52 = -4023.0/686.0 , 117 b52 = -4023.0/686.0 , 118 b53 = 9372.0/1715.0 , 118 b53 = 9372.0/1715.0 , 119 b54 = -66.0/595.0 , 119 b54 = -66.0/595.0 , 120 120 121 b61 = 29.0/238.0 , 121 b61 = 29.0/238.0 , 122 b62 = 0.0 , 122 b62 = 0.0 , 123 b63 = 216.0/385.0 , 123 b63 = 216.0/385.0 , 124 b64 = 54.0/85.0 , 124 b64 = 54.0/85.0 , 125 b65 = -7.0/22.0 , 125 b65 = -7.0/22.0 , 126 126 127 dc1 = 363.0/2975.0 - b61 , 127 dc1 = 363.0/2975.0 - b61 , 128 dc2 = 0.0 - b62 , 128 dc2 = 0.0 - b62 , 129 dc3 = 981.0/1750.0 - b63, 129 dc3 = 981.0/1750.0 - b63, 130 dc4 = 2709.0/4250.0 - b64 , 130 dc4 = 2709.0/4250.0 - b64 , 131 dc5 = -3.0/10.0 - b65 , 131 dc5 = -3.0/10.0 - b65 , 132 dc6 = -1.0/50.0 ; // 132 dc6 = -1.0/50.0 ; // end of declaration 133 133 134 const G4int numberOfVariables = GetNumberOf 134 const G4int numberOfVariables = GetNumberOfVariables(); 135 135 136 // The number of variables to be integrated 136 // The number of variables to be integrated over 137 // 137 // 138 yOut[7] = yTemp[7] = yIn[7]; 138 yOut[7] = yTemp[7] = yIn[7]; 139 139 140 // Saving yInput because yInput and yOut ca 140 // Saving yInput because yInput and yOut can be aliases for same array 141 // 141 // 142 for(i=0; i<numberOfVariables; ++i) 142 for(i=0; i<numberOfVariables; ++i) 143 { 143 { 144 yIn[i]=yInput[i]; 144 yIn[i]=yInput[i]; 145 } 145 } 146 146 147 // RightHandSide(yIn, DyDx) ; // 1st stag 147 // RightHandSide(yIn, DyDx) ; // 1st stage - Not doing, getting passed 148 148 149 for(i=0; i<numberOfVariables; ++i) 149 for(i=0; i<numberOfVariables; ++i) 150 { 150 { 151 yTemp[i] = yIn[i] + b21*Step*DyDx[i] ; 151 yTemp[i] = yIn[i] + b21*Step*DyDx[i] ; 152 } 152 } 153 RightHandSide(yTemp, ak2) ; // 153 RightHandSide(yTemp, ak2) ; // 2nd stage 154 154 155 for(i=0; i<numberOfVariables; ++i) 155 for(i=0; i<numberOfVariables; ++i) 156 { 156 { 157 yTemp[i] = yIn[i] + Step*(b31*DyDx[i] + 157 yTemp[i] = yIn[i] + Step*(b31*DyDx[i] + b32*ak2[i]) ; 158 } 158 } 159 RightHandSide(yTemp, ak3) ; // 159 RightHandSide(yTemp, ak3) ; // 3rd stage 160 160 161 for(i=0; i<numberOfVariables; ++i) 161 for(i=0; i<numberOfVariables; ++i) 162 { 162 { 163 yTemp[i] = yIn[i] + Step*(b41*DyDx[i] + 163 yTemp[i] = yIn[i] + Step*(b41*DyDx[i] + b42*ak2[i] + b43*ak3[i]) ; 164 } 164 } 165 RightHandSide(yTemp, ak4) ; // 165 RightHandSide(yTemp, ak4) ; // 4th stage 166 166 167 for(i=0; i<numberOfVariables; ++i) 167 for(i=0; i<numberOfVariables; ++i) 168 { 168 { 169 yTemp[i] = yIn[i] + Step*(b51*DyDx[i] + 169 yTemp[i] = yIn[i] + Step*(b51*DyDx[i] + b52*ak2[i] 170 + b53*ak3[i] + 170 + b53*ak3[i] + b54*ak4[i]) ; 171 } 171 } 172 RightHandSide(yTemp, ak5) ; // 172 RightHandSide(yTemp, ak5) ; // 5th stage 173 173 174 for(i=0; i<numberOfVariables; ++i) 174 for(i=0; i<numberOfVariables; ++i) 175 { 175 { 176 yOut[i] = yIn[i] + Step*(b61*DyDx[i] + 176 yOut[i] = yIn[i] + Step*(b61*DyDx[i] + b62*ak2[i] + b63*ak3[i] 177 + b64*ak4[i] + b 177 + b64*ak4[i] + b65*ak5[i]) ; 178 } 178 } 179 RightHandSide(yOut, ak6) ; // 179 RightHandSide(yOut, ak6) ; // 6th and Final stage 180 180 181 for(i=0; i<numberOfVariables; ++i) 181 for(i=0; i<numberOfVariables; ++i) 182 { 182 { 183 183 184 yErr[i] = Step*(dc1*DyDx[i] + dc2*ak2[i 184 yErr[i] = Step*(dc1*DyDx[i] + dc2*ak2[i] + dc3*ak3[i] + dc4*ak4[i] 185 + dc5*ak5[i] + dc6*ak6[i] 185 + dc5*ak5[i] + dc6*ak6[i] ) ; 186 186 187 // Store Input and Final values, for po 187 // Store Input and Final values, for possible use in calculating chord 188 // 188 // 189 fLastInitialVector[i] = yIn[i] ; 189 fLastInitialVector[i] = yIn[i] ; 190 fLastFinalVector[i] = yOut[i]; 190 fLastFinalVector[i] = yOut[i]; 191 fLastDyDx[i] = DyDx[i]; 191 fLastDyDx[i] = DyDx[i]; 192 } 192 } 193 193 194 fLastStepLength = Step; 194 fLastStepLength = Step; 195 195 196 return ; 196 return ; 197 } 197 } 198 198 199 // DistChord 199 // DistChord 200 // 200 // 201 G4double G4DoLoMcPriRK34::DistChord() const 201 G4double G4DoLoMcPriRK34::DistChord() const 202 { 202 { 203 G4double distLine, distChord; 203 G4double distLine, distChord; 204 G4ThreeVector initialPoint, finalPoint, mid 204 G4ThreeVector initialPoint, finalPoint, midPoint; 205 205 206 // Store last initial and final points 206 // Store last initial and final points 207 // (they will be overwritten in self-Steppe 207 // (they will be overwritten in self-Stepper call!) 208 // 208 // 209 initialPoint = G4ThreeVector( fLastInitialV 209 initialPoint = G4ThreeVector( fLastInitialVector[0], 210 fLastInitialV 210 fLastInitialVector[1], fLastInitialVector[2] ); 211 finalPoint = G4ThreeVector( fLastFinalVec 211 finalPoint = G4ThreeVector( fLastFinalVector[0], 212 fLastFinalVec 212 fLastFinalVector[1], fLastFinalVector[2] ); 213 213 214 // Do half a Step using StepNoErr 214 // Do half a Step using StepNoErr 215 215 216 fAuxStepper->Stepper( fLastInitialVector, f 216 fAuxStepper->Stepper( fLastInitialVector, fLastDyDx, 0.5 * fLastStepLength, 217 fMidVector, fMidError 217 fMidVector, fMidError ); 218 218 219 midPoint = G4ThreeVector( fMidVector[0], fM 219 midPoint = G4ThreeVector( fMidVector[0], fMidVector[1], fMidVector[2]); 220 220 221 // Use stored values of Initial and Endpoin 221 // Use stored values of Initial and Endpoint + new Midpoint to evaluate 222 // distance of Chord 222 // distance of Chord 223 // 223 // 224 if (initialPoint != finalPoint) 224 if (initialPoint != finalPoint) 225 { 225 { 226 distLine = G4LineSection::Distline( midP 226 distLine = G4LineSection::Distline( midPoint, initialPoint, finalPoint ); 227 distChord = distLine; 227 distChord = distLine; 228 } 228 } 229 else 229 else 230 { 230 { 231 distChord = (midPoint-initialPoint).mag() 231 distChord = (midPoint-initialPoint).mag(); 232 } 232 } 233 return distChord; 233 return distChord; 234 } 234 } 235 235 236 void G4DoLoMcPriRK34::SetupInterpolation() 236 void G4DoLoMcPriRK34::SetupInterpolation() 237 { 237 { 238 } 238 } 239 239 240 void G4DoLoMcPriRK34::SetupInterpolate( const 240 void G4DoLoMcPriRK34::SetupInterpolate( const G4double /* yInput */ [] , 241 const 241 const G4double /* dydx */ [] , 242 const 242 const G4double /* Step */ ) 243 { 243 { 244 // Do Nothing 244 // Do Nothing 245 } 245 } 246 246 247 void G4DoLoMcPriRK34::Interpolate( G4double ta 247 void G4DoLoMcPriRK34::Interpolate( G4double tau, 248 G4double yO 248 G4double yOut[] ) 249 { 249 { 250 Interpolate( fLastInitialVector, fLastDyDx, 250 Interpolate( fLastInitialVector, fLastDyDx, fLastStepLength, yOut, tau ); 251 } 251 } 252 252 253 // Function to evaluate the interpolation at t 253 // Function to evaluate the interpolation at tau fraction of the step 254 // 254 // 255 void G4DoLoMcPriRK34::Interpolate( const G4dou 255 void G4DoLoMcPriRK34::Interpolate( const G4double yInput[], 256 const G4dou 256 const G4double dydx[], 257 const G4dou 257 const G4double Step, 258 G4dou 258 G4double yOut[], 259 G4dou 259 G4double tau ) 260 { 260 { 261 G4double bf1, bf2, bf3, bf4, bf5, bf6; 261 G4double bf1, bf2, bf3, bf4, bf5, bf6; 262 262 263 const G4int numberOfVariables = GetNumberOf 263 const G4int numberOfVariables = GetNumberOfVariables(); 264 264 265 for(G4int i=0; i<numberOfVariables; ++i) 265 for(G4int i=0; i<numberOfVariables; ++i) 266 { 266 { 267 yIn[i]=yInput[i]; 267 yIn[i]=yInput[i]; 268 } 268 } 269 269 270 G4double tau_2 = tau*tau, tau_3 = tau*tau_2 270 G4double tau_2 = tau*tau, tau_3 = tau*tau_2; 271 271 272 // Calculating the polynomials (coefficient 272 // Calculating the polynomials (coefficients for the respective stages) 273 // 273 // 274 bf1 = -(162.0*tau_3 - 504.0*tau_2 + 551.0*t 274 bf1 = -(162.0*tau_3 - 504.0*tau_2 + 551.0*tau - 238.0)/238.0 , 275 bf2 = 0.0 , 275 bf2 = 0.0 , 276 bf3 = 27.0*tau*(27.0*tau_2 - 70.0*tau + 51 276 bf3 = 27.0*tau*(27.0*tau_2 - 70.0*tau + 51.0 )/385.0 , 277 bf4 = -27*tau*(27.0*tau_2 - 50.0*tau + 21.0 277 bf4 = -27*tau*(27.0*tau_2 - 50.0*tau + 21.0)/85.0 , 278 bf5 = 7.0*tau*(2232.0*tau_2 - 4166.0*tau + 278 bf5 = 7.0*tau*(2232.0*tau_2 - 4166.0*tau + 1785.0 )/3278.0 , 279 bf6 = tau*(tau - 1.0)*(387.0*tau - 238.0)/1 279 bf6 = tau*(tau - 1.0)*(387.0*tau - 238.0)/149.0 ; 280 280 281 for( G4int i=0; i<numberOfVariables; ++i) 281 for( G4int i=0; i<numberOfVariables; ++i) 282 { 282 { 283 yOut[i] = yIn[i] + Step*tau*(bf1*dydx[i] 283 yOut[i] = yIn[i] + Step*tau*(bf1*dydx[i] + bf2*ak2[i] + bf3*ak3[i] 284 + bf4*ak4[i] + bf5*ak5[i] + bf6*a 284 + bf4*ak4[i] + bf5*ak5[i] + bf6*ak6[i] ) ; 285 } 285 } 286 } 286 } 287 287