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 // G4MagInt_Driver 27 // 28 // Class description: 29 // 30 // Provides a driver that talks to the Integrator Stepper, and insures that 31 // the error is within acceptable bounds. 32 33 // V.Grichine, 07.10.1996 - Created 34 // W.Wander, 28.01.1998 - Added ability for low order integrators 35 // J.Apostolakis, 08.11.2001 - Respect minimum step in AccurateAdvance 36 // -------------------------------------------------------------------- 37 #ifndef G4MAGINT_DRIVER_HH 38 #define G4MAGINT_DRIVER_HH 39 40 #include "G4VIntegrationDriver.hh" 41 #include "G4MagIntegratorStepper.hh" 42 #include "G4ChordFinderDelegate.hh" 43 44 class G4MagInt_Driver : public G4VIntegrationDriver, 45 public G4ChordFinderDelegate<G4MagInt_Driver> 46 { 47 public: 48 49 G4MagInt_Driver(G4double hminimum, 50 G4MagIntegratorStepper* pItsStepper, 51 G4int numberOfComponents = 6, 52 G4int statisticsVerbosity = 0); 53 ~G4MagInt_Driver() override; 54 // Constructor, destructor. 55 56 G4MagInt_Driver(const G4MagInt_Driver&) = delete; 57 G4MagInt_Driver& operator=(const G4MagInt_Driver&) = delete; 58 59 inline G4double AdvanceChordLimited(G4FieldTrack& track, 60 G4double stepMax, 61 G4double epsStep, 62 G4double chordDistance) override; 63 64 inline void OnStartTracking() override; 65 inline void OnComputeStep(const G4FieldTrack* = nullptr) override {} 66 G4bool DoesReIntegrate() const override { return true; } 67 68 G4bool AccurateAdvance(G4FieldTrack& y_current, 69 G4double hstep, 70 G4double eps, // Requested y_err/hstep 71 G4double hinitial = 0.0) override; 72 // Above drivers for integrator (Runge-Kutta) with stepsize control. 73 // Integrates ODE starting values y_current 74 // from current s (s=s0) to s=s0+h with accuracy eps. 75 // On output ystart is replaced by value at end of interval. 76 // The concept is similar to the odeint routine from NRC p.721-722. 77 78 G4bool QuickAdvance( G4FieldTrack& y_val, // INOUT 79 const G4double dydx[], 80 G4double hstep, 81 G4double& dchord_step, 82 G4double& dyerr) override; 83 // QuickAdvance just tries one Step - it does not ensure accuracy. 84 85 void StreamInfo( std::ostream& os ) const override; 86 // Write out the parameters / state of the driver 87 88 G4bool QuickAdvance( G4FieldTrack& y_posvel, // INOUT 89 const G4double dydx[], 90 G4double hstep, // IN 91 G4double& dchord_step, 92 G4double& dyerr_pos_sq, 93 G4double& dyerr_mom_rel_sq ); 94 // New QuickAdvance that also just tries one Step 95 // (so also does not ensure accuracy) 96 // but does return the errors in position and 97 // momentum (normalised: Delta_Integration(p^2)/(p^2) ) 98 99 inline G4double GetHmin() const; 100 inline G4double Hmin() const; // Obsolete 101 inline G4double GetSafety() const; 102 inline G4double GetPshrnk() const; 103 inline G4double GetPgrow() const; 104 inline G4double GetErrcon() const; 105 void GetDerivatives(const G4FieldTrack& y_curr, // INput 106 G4double dydx[]) const override; // OUTput 107 108 void GetDerivatives(const G4FieldTrack& track, 109 G4double dydx[], 110 G4double field[]) const override; 111 // Accessors 112 113 G4EquationOfMotion* GetEquationOfMotion() override; 114 void SetEquationOfMotion(G4EquationOfMotion* equation) override; 115 116 void RenewStepperAndAdjust(G4MagIntegratorStepper* pItsStepper) override; 117 // Sets a new stepper pItsStepper for this driver. Then it calls 118 // ReSetParameters to reset its parameters accordingly. 119 120 inline void ReSetParameters(G4double new_safety = 0.9); 121 // i) sets the exponents (pgrow & pshrnk), 122 // using the current Stepper's order, 123 // ii) sets the safety 124 // ii) calculates "errcon" according to the above values. 125 126 inline void SetSafety(G4double valS); 127 inline void SetPshrnk(G4double valPs); 128 inline void SetPgrow (G4double valPg); 129 inline void SetErrcon(G4double valEc); 130 // When setting safety or pgrow, errcon will be set to a compatible value. 131 132 inline G4double ComputeAndSetErrcon(); 133 134 const G4MagIntegratorStepper* GetStepper() const override; 135 G4MagIntegratorStepper* GetStepper() override; 136 137 void OneGoodStep( G4double ystart[], // Like old RKF45step() 138 const G4double dydx[], 139 G4double& x, 140 G4double htry, 141 G4double eps, // memb variables ? 142 G4double& hdid, 143 G4double& hnext ) ; 144 // This takes one Step that is as large as possible while 145 // satisfying the accuracy criterion of: 146 // yerr < eps * |y_end-y_start| 147 148 G4double ComputeNewStepSize(G4double errMaxNorm, // normalised 149 G4double hstepCurrent) override; 150 // Taking the last step's normalised error, calculate 151 // a step size for the next step. 152 // Does it limit the next step's size within a factor of the current? 153 // -- DOES NOT limit for very bad steps 154 // -- DOES limit for very good (x5) 155 156 G4double 157 ComputeNewStepSize_WithoutReductionLimit(G4double errMaxNorm, 158 G4double hstepCurrent); 159 // Taking the last step's normalised error, calculate 160 // a step size for the next step. 161 // Do not limit the next step's size within a factor of the 162 // current one when *reducing* the size, i.e. for badly failing steps. 163 164 G4double ComputeNewStepSize_WithinLimits(G4double errMaxNorm, // normalised 165 G4double hstepCurrent); 166 // Taking the last step's normalised error, calculate 167 // a step size for the next step. 168 // Limit the next step's size within a range around the current one. 169 170 inline G4int GetMaxNoSteps() const; 171 inline void SetMaxNoSteps(G4int val); 172 // Modify and Get the Maximum number of Steps that can be 173 // taken for the integration of a single segment - 174 // (i.e. a single call to AccurateAdvance). 175 176 inline void SetHmin(G4double newval); 177 void SetVerboseLevel(G4int newLevel) override; 178 G4int GetVerboseLevel() const override; 179 180 inline G4double GetSmallestFraction() const; 181 void SetSmallestFraction( G4double val ); 182 183 protected: 184 185 void WarnSmallStepSize(G4double hnext, G4double hstep, 186 G4double h, G4double xDone, 187 G4int noSteps); 188 189 void WarnTooManyStep(G4double x1start, G4double x2end, G4double xCurrent); 190 void WarnEndPointTooFar(G4double endPointDist, 191 G4double hStepSize , 192 G4double epsilonRelative, 193 G4int debugFlag); 194 // Issue warnings for undesirable situations 195 196 void PrintStatus(const G4double* StartArr, 197 G4double xstart, 198 const G4double* CurrentArr, 199 G4double xcurrent, 200 G4double requestStep, 201 G4int subStepNo); 202 void PrintStatus(const G4FieldTrack& StartFT, 203 const G4FieldTrack& CurrentFT, 204 G4double requestStep, 205 G4int subStepNo); 206 void PrintStat_Aux(const G4FieldTrack& aFieldTrack, 207 G4double requestStep, 208 G4double actualStep, 209 G4int subStepNo, 210 G4double subStepSize, 211 G4double dotVelocities); 212 // Verbose output for debugging 213 214 void PrintStatisticsReport(); 215 // Report on the number of steps, maximum errors etc. 216 217 #ifdef QUICK_ADV_TWO 218 G4bool QuickAdvance( G4double yarrin[], // In 219 const G4double dydx[], 220 G4double hstep, 221 G4double yarrout[], // Out 222 G4double& dchord_step, // Out 223 G4double& dyerr ); // in length 224 #endif 225 226 private: 227 228 // --------------------------------------------------------------- 229 // INVARIANTS 230 231 G4double fMinimumStep = 0.0; 232 // Minimum Step allowed in a Step (in absolute units) 233 G4double fSmallestFraction = 1.0e-12; // Expected range 1e-12 to 5e-15 234 // Smallest fraction of (existing) curve length - in relative units 235 // below this fraction the current step will be the last 236 237 const G4int fNoIntegrationVariables = 0; // Variables in integration 238 const G4int fMinNoVars = 12; // Minimum number for FieldTrack 239 const G4int fNoVars = 0; // Full number of variable 240 241 G4int fMaxNoSteps; 242 G4int fMaxStepBase = 250; // was 5000 243 // Default maximum number of steps is Base divided by the order of Stepper 244 245 G4double safety; 246 G4double pshrnk; // exponent for shrinking 247 G4double pgrow; // exponent for growth 248 G4double errcon; 249 // Parameters used to grow and shrink trial stepsize. 250 251 G4int fStatisticsVerboseLevel = 0; 252 253 // --------------------------------------------------------------- 254 // DEPENDENT Objects 255 256 G4MagIntegratorStepper* pIntStepper = nullptr; 257 258 // --------------------------------------------------------------- 259 // STATE 260 261 unsigned long fNoTotalSteps=0, fNoBadSteps=0; 262 unsigned long fNoSmallSteps=0, fNoInitialSmallSteps=0, fNoCalls=0; 263 G4double fDyerr_max=0.0, fDyerr_mx2=0.0; 264 G4double fDyerrPos_smTot=0.0, fDyerrPos_lgTot=0.0, fDyerrVel_lgTot=0.0; 265 G4double fSumH_sm=0.0, fSumH_lg=0.0; 266 // Step Statistics 267 268 G4int fVerboseLevel = 0; // Verbosity level for printing (debug, ..) 269 // Could be varied during tracking - to help identify issues 270 271 using ChordFinderDelegate = G4ChordFinderDelegate<G4MagInt_Driver>; 272 }; 273 274 #include "G4MagIntegratorDriver.icc" 275 276 #endif 277