include/Geant4/G4OldMagIntDriver.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 // G4OldMagIntDriver
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 G4OLD_MAGINT_DRIVER_HH
0038 #define G4OLD_MAGINT_DRIVER_HH
0039
0040 #include "G4VIntegrationDriver.hh"
0041 #include "G4MagIntegratorStepper.hh"
0042 #include "G4ChordFinderDelegate.hh"
0043
0044 class G4OldMagIntDriver : public G4VIntegrationDriver,
0045                           public G4ChordFinderDelegate<G4OldMagIntDriver>
0046 {
0047   public:
0048
0049     G4OldMagIntDriver(G4double hminimum,
0050                       G4MagIntegratorStepper* pItsStepper,
0051                       G4int numberOfComponents = 6,
0052                       G4int statisticsVerbosity = 0);
0053    ~G4OldMagIntDriver() override;
0054       // Constructor, destructor.
0055
0056     G4OldMagIntDriver(const G4OldMagIntDriver&) = delete;
0057     G4OldMagIntDriver& operator=(const G4OldMagIntDriver&) = 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     inline 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     G4bool QuickAdvance(      G4FieldTrack& y_posvel,   // INOUT
0086                         const G4double dydx[],
0087                               G4double hstep,           // IN
0088                               G4double& dchord_step,
0089                               G4double& dyerr_pos_sq,
0090                               G4double& dyerr_mom_rel_sq);
0091       // QuickAdvance that also just tries one Step (so also does not
0092       // ensure accuracy), but does return the errors in  position and
0093       // momentum (normalised: Delta_Integration(p^2)/(p^2) ).
0094
0095     inline G4double GetHmin() const;
0096     inline G4double Hmin() const;     // Obsolete
0097     inline G4double GetSafety() const;
0098     inline G4double GetPshrnk() const;
0099     inline G4double GetPgrow() const;
0100     inline G4double GetErrcon() const;
0101     void GetDerivatives(const G4FieldTrack& y_curr,            // INput
0102                               G4double dydx[]) const override; // OUTput
0103
0104     void GetDerivatives(const G4FieldTrack& track,
0105                               G4double dydx[],
0106                               G4double field[]) const override;
0107     // Accessors
0108
0109     G4EquationOfMotion* GetEquationOfMotion() override;
0110     void SetEquationOfMotion(G4EquationOfMotion* equation) override;
0111
0112     void RenewStepperAndAdjust(G4MagIntegratorStepper* pItsStepper) override;
0113       // Sets a new stepper pItsStepper for this driver. Then it calls
0114       // ReSetParameters to reset its parameters accordingly.
0115
0116     inline void ReSetParameters(G4double new_safety = 0.9);
0117       //  i) sets the exponents (pgrow & pshrnk),
0118       //     using the current Stepper's order,
0119       // ii) sets the safety
0120       // ii) calculates "errcon" according to the above values.
0121
0122     inline void SetSafety(G4double valS);
0123     inline void SetPshrnk(G4double valPs);
0124     inline void SetPgrow (G4double valPg);
0125     inline void SetErrcon(G4double valEc);
0126       // When setting safety or pgrow, errcon will be set to a compatible value.
0127
0128     inline G4double ComputeAndSetErrcon();
0129
0130     const G4MagIntegratorStepper* GetStepper() const override;
0131           G4MagIntegratorStepper* GetStepper() override;
0132
0133     void OneGoodStep(      G4double  ystart[], // Like old RKF45step()
0134                      const G4double  dydx[],
0135                            G4double& x,
0136                            G4double htry,
0137                            G4double  eps,      //  memb variables ?
0138                            G4double& hdid,
0139                            G4double& hnext) ;
0140       // This takes one Step that is as large as possible while
0141       // satisfying the accuracy criterion of:
0142       // yerr < eps * |y_end-y_start|
0143
0144     G4double ComputeNewStepSize(G4double errMaxNorm, // normalised
0145                                 G4double hstepCurrent) override;
0146       // Taking the last step's normalised error, calculate
0147       // a step size for the next step.
0148       // Do not limit the next step's size within a factor of the
0149       // current one.
0150
0151     void StreamInfo( std::ostream& os ) const override;
0152
0153     G4double ComputeNewStepSize_WithinLimits(G4double errMaxNorm, // normalised
0154                                              G4double hstepCurrent);
0155       // Taking the last step's normalised error, calculate
0156       // a step size for the next step.
0157       // Limit the next step's size within a range around the current one.
0158
0159     inline G4int GetMaxNoSteps() const;
0160     inline void SetMaxNoSteps(G4int val);
0161       // Modify and Get the Maximum number of Steps that can be
0162       // taken for the integration of a single segment -
0163       // (i.e. a single call to AccurateAdvance).
0164
0165     inline void SetHmin(G4double newval);
0166     void SetVerboseLevel(G4int newLevel) override;
0167     G4int GetVerboseLevel() const override;
0168
0169     inline G4double GetSmallestFraction() const;
0170     void SetSmallestFraction( G4double val );
0171
0172   protected:
0173
0174     void WarnSmallStepSize(G4double hnext, G4double hstep,
0175                            G4double h, G4double xDone,
0176                            G4int noSteps);
0177
0178     void WarnTooManyStep(G4double x1start, G4double x2end, G4double xCurrent);
0179     void WarnEndPointTooFar(G4double endPointDist,
0180                             G4double hStepSize ,
0181                             G4double epsilonRelative,
0182                             G4int debugFlag);
0183       // Issue warnings for undesirable situations
0184
0185     void PrintStatus(const G4double* StartArr,
0186                            G4double xstart,
0187                      const G4double* CurrentArr,
0188                            G4double xcurrent,
0189                            G4double requestStep,
0190                            G4int subStepNo);
0191     void PrintStatus(const G4FieldTrack& StartFT,
0192                      const G4FieldTrack& CurrentFT,
0193                            G4double requestStep,
0194                            G4int subStepNo);
0195     void PrintStat_Aux(const G4FieldTrack& aFieldTrack,
0196                              G4double requestStep,
0197                              G4double actualStep,
0198                              G4int subStepNo,
0199                              G4double subStepSize,
0200                              G4double dotVelocities);
0201       // Verbose output for debugging
0202
0203     void PrintStatisticsReport();
0204       // Report on the number of steps, maximum errors etc.
0205
0206 #ifdef QUICK_ADV_TWO
0207      G4bool QuickAdvance(      G4double  yarrin[],     // In
0208                          const G4double  dydx[],
0209                                G4double  hstep,
0210                                G4double  yarrout[],    // Out
0211                                G4double& dchord_step,  // Out
0212                                G4double& dyerr );      // in length
0213 #endif
0214
0215   private:
0216
0217     // ---------------------------------------------------------------
0218     //  INVARIANTS
0219
0220     G4double fMinimumStep = 0.0;
0221       // Minimum Step allowed in a Step (in absolute units)
0222     G4double fSmallestFraction = 1.0e-12;    // Expected range 1e-12 to 5e-15
0223       // Smallest fraction of (existing) curve length - in relative units
0224       //  below this fraction the current step will be the last
0225
0226     const G4int fNoIntegrationVariables = 0; // Variables in integration
0227     const G4int fMinNoVars = 12;             // Minimum number for FieldTrack
0228     const G4int fNoVars = 0;                 // Full number of variable
0229
0230     G4int fMaxNoSteps;
0231     G4int fMaxStepBase = 250;  // was 5000
0232       // Default maximum number of steps is Base divided by the order of Stepper
0233
0234     G4double safety;
0235     G4double pshrnk;   //  exponent for shrinking
0236     G4double pgrow;    //  exponent for growth
0237     G4double errcon;
0238     // Parameters used to grow and shrink trial stepsize.
0239
0240     G4int fStatisticsVerboseLevel = 0;
0241
0242     // ---------------------------------------------------------------
0243     // DEPENDENT Objects
0244
0245     G4MagIntegratorStepper* pIntStepper = nullptr;
0246
0247     // ---------------------------------------------------------------
0248     //  STATE
0249
0250     unsigned long fNoTotalSteps=0, fNoBadSteps=0;
0251     unsigned long fNoSmallSteps=0, fNoInitialSmallSteps=0, fNoCalls=0;
0252     G4double fDyerr_max=0.0, fDyerr_mx2=0.0;
0253     G4double fDyerrPos_smTot=0.0, fDyerrPos_lgTot=0.0, fDyerrVel_lgTot=0.0;
0254     G4double fSumH_sm=0.0, fSumH_lg=0.0;
0255       // Step Statistics
0256
0257     G4int fVerboseLevel = 0;   // Verbosity level for printing (debug, ..)
0258       // Could be varied during tracking - to help identify issues
0259
0260     using ChordFinderDelegate = G4ChordFinderDelegate<G4OldMagIntDriver>;
0261 };
0262
0263 #include "G4OldMagIntDriver.icc"
0264
0265 #endif