Geant4 Cross Reference |
1 // ******************************************* 1 2 // * License and Disclaimer 3 // * 4 // * The Geant4 software is copyright of th 5 // * the Geant4 Collaboration. It is provided 6 // * conditions of the Geant4 Software License 7 // * LICENSE and available at http://cern.ch/ 8 // * include a list of copyright holders. 9 // * 10 // * Neither the authors of this software syst 11 // * institutes,nor the agencies providing fin 12 // * work make any representation or warran 13 // * regarding this software system or assum 14 // * use. Please see the license in the file 15 // * for the full disclaimer and the limitatio 16 // * 17 // * This code implementation is the result 18 // * technical work of the GEANT4 collaboratio 19 // * By using, copying, modifying or distri 20 // * any work based on the software) you ag 21 // * use in resulting scientific publicati 22 // * acceptance of all terms of the Geant4 Sof 23 // ******************************************* 24 // 25 // G4BorisDriver 26 // 27 // Class description: 28 // 29 // G4BorisDriver is a driver class using the 30 // method to integrate the equation of motion. 31 // 32 // 33 // Author: Divyansh Tiwari, Google Summer of C 34 // Supervision: John Apostolakis,Renee Fatemi, 35 // ------------------------------------------- 36 #include <cassert> 37 38 #include "G4BorisDriver.hh" 39 40 #include "G4SystemOfUnits.hh" 41 #include "G4LineSection.hh" 42 #include "G4FieldUtils.hh" 43 44 const G4double G4BorisDriver::fErrorConstraint 45 fMaxSteppingDecrease / fSafetyFact 46 47 const G4double G4BorisDriver::fErrorConstraint 48 fMaxSteppingIncrease / fSafetyFact 49 50 // ------------------------------------------- 51 52 G4BorisDriver:: 53 G4BorisDriver( G4double hminimum, G4BorisSchem 54 G4int numberOfComponents, G4boo 55 : fMinimumStep(hminimum), 56 fVerbosity(verbosity), 57 boris(Boris) 58 // , interval_sequence{2,4} 59 { 60 assert(boris->GetNumberOfVariables() == nu 61 62 if(boris->GetNumberOfVariables() != number 63 { 64 std::ostringstream msg; 65 msg << "Disagreement in number of variab 66 << boris->GetNumberOfVariables() 67 << " vs no of components = " << numb 68 G4Exception("G4BorisDriver Constructor:" 69 "GeomField1001", FatalExcept 70 } 71 72 } 73 74 // ------------------------------------------- 75 76 G4bool G4BorisDriver::AccurateAdvance( G4Field 77 G4doubl 78 G4doubl 79 G4doubl 80 { 81 // Specification: Driver with adaptive step 82 // Integrate starting values at y_current o 83 // On output 'track' is replaced by values 84 85 // Ensure that hstep > 0 86 if(hstep == 0) 87 { 88 std::ostringstream message; 89 message << "Proposed step is zero; hstep 90 G4Exception("G4BorisDriver::AccurateAdva 91 "GeomField1001", JustWarning 92 return true; 93 } 94 if(hstep < 0) 95 { 96 std::ostringstream message; 97 message << "Invalid run condition." << G 98 << "Proposed step is negative; h 99 << "Requested step cannot be neg 100 G4Exception("G4BorisDriver::AccurateAdva 101 "GeomField0003", EventMustBe 102 return false; 103 } 104 105 if( hinitial == 0.0 ) { hinitial = hstep; } 106 if( hinitial < 0.0 ) { hinitial = std::fabs 107 // G4double htrial = std::min( hstep, hinit 108 G4double htrial = hstep; 109 // Decide first step size 110 111 // G4int noOfSteps = h/hstep; 112 113 // integration variables 114 // 115 track.DumpToArray(yCurrent); 116 117 const G4double restMass = track.GetRestMass 118 const G4double charge = track.GetCharge()*e 119 const G4int nvar= GetNumberOfVariables() 120 121 // copy non-integration variables to out ar 122 // 123 std::memcpy(yOut + nvar, 124 yCurrent + nvar, 125 sizeof(G4double)*(G4FieldTrack: 126 127 G4double curveLength = track.GetCurveLength 128 const G4double endCurveLength = curveLength 129 130 // -- Initial version: Did it in one step - 131 // G4FieldTrack yFldTrk(track); 132 // yFldTrk.LoadFromArray(yCurrent, G4FieldT 133 // yFldTrk.SetCurveLength(curveLength); 134 // G4double dchord_step, dyerr_len; 135 // QuickAdvance(yFldTrk, dydxCurrent, htria 136 137 const G4double hThreshold = 138 std::max(epsilon * hstep, fSmallestFra 139 140 G4double htry= htrial; 141 142 for (G4int nstp = 0; nstp < fMaxNoSteps; ++ 143 { 144 G4double hdid= 0.0, hnext=0.0; 145 146 OneGoodStep(yCurrent, curveLength, htry, 147 //********* 148 149 // Simple check: move (distance of displ 150 const G4ThreeVector StartPos = field_uti 151 const G4ThreeVector EndPos = field_utils 152 CheckStep(EndPos, StartPos, hdid); 153 154 // Check 1. for finish and 2. *av 155 if (curveLength >= endCurveLength || htr 156 { 157 break; 158 } 159 160 htry = std::max(hnext, fMinimumStep); 161 if (curveLength + htry > endCurveLength) 162 { 163 htry = endCurveLength - curveLength; 164 } 165 } 166 167 // upload new state 168 track.LoadFromArray(yCurrent, G4FieldTrack: 169 track.SetCurveLength(curveLength); 170 171 return true; 172 } 173 174 // ------------------------------------------- 175 176 void G4BorisDriver::OneGoodStep(G4double y[], 177 G4double& curv 178 G4double htry 179 G4double epsi 180 G4double rest 181 G4double char 182 G4double& hdid 183 G4double& hnex 184 { 185 G4double error2 = DBL_MAX; 186 G4double yerr[G4FieldTrack::ncompSVEC], yt 187 188 G4double h = htry; 189 190 const G4int max_trials = 100; 191 192 for (G4int iter = 0; iter < max_trials; ++ 193 { 194 boris->StepWithErrorEstimate(y, restMa 195 196 error2 = field_utils::relativeError2(y 197 e 198 if (error2 <= 1.0) 199 { 200 break; 201 } 202 203 h = ShrinkStepSize2(h, error2); 204 205 G4double xnew = curveLength + h; 206 if(xnew == curveLength) 207 { 208 std::ostringstream message; 209 message << "Stepsize underflow in 210 << "- Step's start x=" << 211 << " and end x= " << xnew 212 << " are equal !! " << G4e 213 << " Due to step-size= " 214 << ". Note that input step 215 G4Exception("G4IntegrationDriver:: 216 "GeomField1001", JustW 217 break; 218 } 219 } 220 221 hnext = GrowStepSize2(h, error2); 222 curveLength += (hdid = h); 223 224 field_utils::copy(y, ytemp, GetNumberOfVar 225 } 226 227 // ===========-------------------------------- 228 229 G4bool G4BorisDriver:: 230 QuickAdvance( G4FieldTrack& track, const G4dou 231 G4double hstep, G4double& missDi 232 { 233 const auto nvar = boris->GetNumberOfVariab 234 235 track.DumpToArray(yIn); 236 const G4double curveLength = track.GetCurv 237 238 // call the boris method for step length h 239 G4double restMass = track.GetRestMass(); 240 G4double charge = track.GetCharge()*e_SI; 241 242 // boris->DoStep(restMass, charge, yIn, y 243 // boris->DoStep(restMass, charge, yMid, y 244 boris->StepWithMidAndErrorEstimate(yIn, re 245 yMid, y 246 // Same, and also return mid-point evalu 247 248 // How to calculate chord length?? 249 const auto mid = field_utils::makeVector(y 250 field_utils::Value3D::Pos 251 const auto in = field_utils::makeVector(y 252 field_utils::Value3D::Pos 253 const auto out = field_utils::makeVector(y 254 field_utils::Value3D::Pos 255 256 missDist = G4LineSection::Distline(mid, in 257 258 dyerr = field_utils::absoluteError(yOut, y 259 260 // copy non-integrated variables to output 261 // 262 std::memcpy(yOut + nvar, yIn + nvar, 263 sizeof(G4double) * (G4FieldTra 264 265 // set new state 266 // 267 track.LoadFromArray(yOut, G4FieldTrack::nc 268 track.SetCurveLength(curveLength + hstep) 269 270 return true; 271 } 272 273 // ------------------------------------------- 274 275 void G4BorisDriver:: 276 GetDerivatives( const G4FieldTrack& yTrack, G4 277 { 278 G4double EBfieldValue[6]; 279 GetDerivatives(yTrack, dydx, EBfieldValue); 280 } 281 282 // ------------------------------------------- 283 284 void G4BorisDriver:: 285 GetDerivatives( const G4FieldTrack& yTrack, G4 286 G4double EBfieldValue[]) const 287 { 288 // G4Exception("G4BorisDriver::GetDerivativ 289 // "GeomField0003", FatalExcept 290 G4double ytemp[G4FieldTrack::ncompSVEC]; 291 yTrack.DumpToArray(ytemp); 292 GetEquationOfMotion()->EvaluateRhsReturnB(y 293 } 294 295 // ------------------------------------------- 296 297 G4double G4BorisDriver::ShrinkStepSize2(G4doub 298 { 299 if (error2 > fErrorConstraintShrink * fErro 300 { 301 return fMaxSteppingDecrease * h; 302 } 303 return fSafetyFactor * h * std::pow(error2, 304 } 305 306 // ------------------------------------------- 307 308 G4double G4BorisDriver::GrowStepSize2(G4double 309 // Given the square of the relative error, 310 { 311 if (error2 < fErrorConstraintGrow * fError 312 { 313 return fMaxSteppingIncrease * h; 314 } 315 return fSafetyFactor * h * std::pow(error2 316 } 317 318 // ------------------------------------------- 319 320 void G4BorisDriver::SetEquationOfMotion(G4Equa 321 { 322 G4Exception("G4BorisDriver::SetEquationOfMo 323 "This method is not implemented 324 } 325 326 // ------------------------------------------- 327 328 void 329 G4BorisDriver::StreamInfo( std::ostream& os ) 330 { 331 os << "State of G4BorisDriver: " << std::end 332 os << " Method is implemented, but gives n 333 } 334