include/Geant4/G4MagIntegratorDriver.hh

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 // G4MagInt_Driver
0027 //
0028 // Class description:
0029 //
0030 // Provides a driver that talks to the Integrator Stepper, and insures that
0031 // the error is within acceptable bounds.
0032
0033 // V.Grichine, 07.10.1996 - Created
0034 // W.Wander, 28.01.1998 - Added ability for low order integrators
0035 // J.Apostolakis, 08.11.2001 - Respect minimum step in AccurateAdvance
0036 // --------------------------------------------------------------------
0037 #ifndef G4MAGINT_DRIVER_HH
0038 #define G4MAGINT_DRIVER_HH
0039
0040 #include "G4VIntegrationDriver.hh"
0041 #include "G4MagIntegratorStepper.hh"
0042 #include "G4ChordFinderDelegate.hh"
0043
0044 class G4MagInt_Driver : public G4VIntegrationDriver,
0045                         public G4ChordFinderDelegate<G4MagInt_Driver>
0046 {
0047   public:
0048
0049     G4MagInt_Driver(G4double hminimum,
0050                     G4MagIntegratorStepper* pItsStepper,
0051                     G4int numberOfComponents = 6,
0052                     G4int statisticsVerbosity = 0);
0053    ~G4MagInt_Driver() override;
0054       // Constructor, destructor.
0055
0056     G4MagInt_Driver(const G4MagInt_Driver&) = delete;
0057     G4MagInt_Driver& operator=(const G4MagInt_Driver&) = delete;
0058
0059     inline G4double AdvanceChordLimited(G4FieldTrack& track,
0060                                         G4double stepMax,
0061                                         G4double epsStep,
0062                                         G4double chordDistance) override;
0063
0064     inline void OnStartTracking() override;
0065     inline void OnComputeStep(const G4FieldTrack* = nullptr) override {}
0066     G4bool DoesReIntegrate() const override { return true; }
0067
0068     G4bool AccurateAdvance(G4FieldTrack& y_current,
0069                            G4double hstep,
0070                            G4double eps,  // Requested y_err/hstep
0071                            G4double hinitial = 0.0) override;
0072       // Above drivers for integrator (Runge-Kutta) with stepsize control.
0073       // Integrates ODE starting values y_current
0074       // from current s (s=s0) to s=s0+h with accuracy eps.
0075       // On output ystart is replaced by value at end of interval.
0076       // The concept is similar to the odeint routine from NRC p.721-722.
0077
0078     G4bool QuickAdvance(      G4FieldTrack& y_val,      // INOUT
0079                         const G4double dydx[],
0080                               G4double hstep,
0081                               G4double& dchord_step,
0082                               G4double& dyerr) override;
0083       // QuickAdvance just tries one Step - it does not ensure accuracy.
0084
0085     void  StreamInfo( std::ostream& os ) const override;
0086      // Write out the parameters / state of the driver
0087
0088     G4bool QuickAdvance(      G4FieldTrack& y_posvel,   // INOUT
0089                         const G4double dydx[],
0090                               G4double hstep,           // IN
0091                               G4double& dchord_step,
0092                               G4double& dyerr_pos_sq,
0093                               G4double& dyerr_mom_rel_sq );
0094       // New QuickAdvance that also just tries one Step
0095       //    (so also does not ensure accuracy)
0096       //    but does return the errors in  position and
0097       //        momentum (normalised: Delta_Integration(p^2)/(p^2) )
0098
0099     inline G4double GetHmin() const;
0100     inline G4double Hmin() const;     // Obsolete
0101     inline G4double GetSafety() const;
0102     inline G4double GetPshrnk() const;
0103     inline G4double GetPgrow() const;
0104     inline G4double GetErrcon() const;
0105     void GetDerivatives(const G4FieldTrack& y_curr,            // INput
0106                               G4double dydx[]) const override; // OUTput
0107
0108     void GetDerivatives(const G4FieldTrack& track,
0109                               G4double dydx[],
0110                               G4double field[]) const override;
0111     // Accessors
0112
0113     G4EquationOfMotion* GetEquationOfMotion() override;
0114     void SetEquationOfMotion(G4EquationOfMotion* equation) override;
0115
0116     void RenewStepperAndAdjust(G4MagIntegratorStepper* pItsStepper) override;
0117       // Sets a new stepper pItsStepper for this driver. Then it calls
0118       // ReSetParameters to reset its parameters accordingly.
0119
0120     inline void ReSetParameters(G4double new_safety = 0.9);
0121       //  i) sets the exponents (pgrow & pshrnk),
0122       //     using the current Stepper's order,
0123       // ii) sets the safety
0124       // ii) calculates "errcon" according to the above values.
0125
0126     inline void SetSafety(G4double valS);
0127     inline void SetPshrnk(G4double valPs);
0128     inline void SetPgrow (G4double valPg);
0129     inline void SetErrcon(G4double valEc);
0130       // When setting safety or pgrow, errcon will be set to a compatible value.
0131
0132     inline G4double ComputeAndSetErrcon();
0133
0134     const G4MagIntegratorStepper* GetStepper() const override;
0135           G4MagIntegratorStepper* GetStepper() override;
0136
0137     void OneGoodStep(      G4double  ystart[], // Like old RKF45step()
0138                      const G4double  dydx[],
0139                            G4double& x,
0140                            G4double htry,
0141                            G4double  eps,      //  memb variables ?
0142                            G4double& hdid,
0143                            G4double& hnext ) ;
0144       // This takes one Step that is as large as possible while
0145       // satisfying the accuracy criterion of:
0146       // yerr < eps * |y_end-y_start|
0147
0148     G4double ComputeNewStepSize(G4double errMaxNorm, // normalised
0149                                 G4double hstepCurrent) override;
0150       // Taking the last step's normalised error, calculate
0151       // a step size for the next step.
0152       // Does it limit the next step's size within a factor of the current?
0153       // --  DOES NOT limit for very bad steps
0154       // --  DOES     limit for very good (x5)
0155
0156     G4double
0157     ComputeNewStepSize_WithoutReductionLimit(G4double  errMaxNorm,
0158                                              G4double hstepCurrent);
0159       // Taking the last step's normalised error, calculate
0160       // a step size for the next step.
0161       // Do not limit the next step's size within a factor of the
0162       //  current one when *reducing* the size, i.e. for badly failing steps.
0163
0164     G4double ComputeNewStepSize_WithinLimits(G4double errMaxNorm, // normalised
0165                                              G4double hstepCurrent);
0166       // Taking the last step's normalised error, calculate
0167       // a step size for the next step.
0168       // Limit the next step's size within a range around the current one.
0169
0170     inline G4int GetMaxNoSteps() const;
0171     inline void SetMaxNoSteps(G4int val);
0172       // Modify and Get the Maximum number of Steps that can be
0173       // taken for the integration of a single segment -
0174       // (i.e. a single call to AccurateAdvance).
0175
0176     inline void SetHmin(G4double newval);
0177     void SetVerboseLevel(G4int newLevel) override;
0178     G4int GetVerboseLevel() const override;
0179
0180     inline G4double GetSmallestFraction() const;
0181     void SetSmallestFraction( G4double val );
0182
0183   protected:
0184
0185     void WarnSmallStepSize(G4double hnext, G4double hstep,
0186                            G4double h, G4double xDone,
0187                            G4int noSteps);
0188
0189     void WarnTooManyStep(G4double x1start, G4double x2end, G4double xCurrent);
0190     void WarnEndPointTooFar(G4double endPointDist,
0191                             G4double hStepSize ,
0192                             G4double epsilonRelative,
0193                             G4int debugFlag);
0194       // Issue warnings for undesirable situations
0195
0196     void PrintStatus(const G4double* StartArr,
0197                            G4double xstart,
0198                      const G4double* CurrentArr,
0199                            G4double xcurrent,
0200                            G4double requestStep,
0201                            G4int subStepNo);
0202     void PrintStatus(const G4FieldTrack& StartFT,
0203                      const G4FieldTrack& CurrentFT,
0204                            G4double requestStep,
0205                            G4int subStepNo);
0206     void PrintStat_Aux(const G4FieldTrack& aFieldTrack,
0207                              G4double requestStep,
0208                              G4double actualStep,
0209                              G4int subStepNo,
0210                              G4double subStepSize,
0211                              G4double dotVelocities);
0212       // Verbose output for debugging
0213
0214     void PrintStatisticsReport();
0215       // Report on the number of steps, maximum errors etc.
0216
0217 #ifdef QUICK_ADV_TWO
0218     G4bool QuickAdvance(      G4double  yarrin[],     // In
0219                         const G4double  dydx[],
0220                               G4double  hstep,
0221                               G4double  yarrout[],    // Out
0222                               G4double& dchord_step,  // Out
0223                               G4double& dyerr );      // in length
0224 #endif
0225
0226   private:
0227
0228     // ---------------------------------------------------------------
0229     //  INVARIANTS
0230
0231     G4double fMinimumStep = 0.0;
0232       // Minimum Step allowed in a Step (in absolute units)
0233     G4double fSmallestFraction = 1.0e-12;    // Expected range 1e-12 to 5e-15
0234       // Smallest fraction of (existing) curve length - in relative units
0235       //  below this fraction the current step will be the last
0236
0237     const G4int fNoIntegrationVariables = 0; // Variables in integration
0238     const G4int fMinNoVars = 12;             // Minimum number for FieldTrack
0239     const G4int fNoVars = 0;                 // Full number of variable
0240
0241     G4int fMaxNoSteps;
0242     G4int fMaxStepBase = 250;  // was 5000
0243       // Default maximum number of steps is Base divided by the order of Stepper
0244
0245     G4double safety;
0246     G4double pshrnk;   //  exponent for shrinking
0247     G4double pgrow;    //  exponent for growth
0248     G4double errcon;
0249     // Parameters used to grow and shrink trial stepsize.
0250
0251     G4int fStatisticsVerboseLevel = 0;
0252
0253     // ---------------------------------------------------------------
0254     // DEPENDENT Objects
0255
0256     G4MagIntegratorStepper* pIntStepper = nullptr;
0257
0258     // ---------------------------------------------------------------
0259     //  STATE
0260
0261     unsigned long fNoTotalSteps=0, fNoBadSteps=0;
0262     unsigned long fNoSmallSteps=0, fNoInitialSmallSteps=0, fNoCalls=0;
0263     G4double fDyerr_max=0.0, fDyerr_mx2=0.0;
0264     G4double fDyerrPos_smTot=0.0, fDyerrPos_lgTot=0.0, fDyerrVel_lgTot=0.0;
0265     G4double fSumH_sm=0.0, fSumH_lg=0.0;
0266       // Step Statistics
0267
0268     G4int fVerboseLevel = 0;   // Verbosity level for printing (debug, ..)
0269       // Could be varied during tracking - to help identify issues
0270
0271     using ChordFinderDelegate = G4ChordFinderDelegate<G4MagInt_Driver>;
0272 };
0273
0274 #include "G4MagIntegratorDriver.icc"
0275
0276 #endif