Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/geometry/magneticfield/include/G4MagIntegratorDriver.hh

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  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