include/Geant4/G4BorisDriver.hh

0001 // ********************************************************************
0002 // * License and Disclaimer                                           *
0003 // *                                                                  *
0004 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0005 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0006 // * conditions of the Geant4 Software License,  included in the file *
0007 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0008 // * include a list of copyright holders.                             *
0009 // *                                                                  *
0010 // * Neither the authors of this software system, nor their employing *
0011 // * institutes,nor the agencies providing financial support for this *
0012 // * work  make  any representation or  warranty, express or implied, *
0013 // * regarding  this  software system or assume any liability for its *
0014 // * use.  Please see the license in the file  LICENSE  and URL above *
0015 // * for the full disclaimer and the limitation of liability.         *
0016 // *                                                                  *
0017 // * This  code  implementation is the result of  the  scientific and *
0018 // * technical work of the GEANT4 collaboration.                      *
0019 // * By using,  copying,  modifying or  distributing the software (or *
0020 // * any work based  on the software)  you  agree  to acknowledge its *
0021 // * use  in  resulting  scientific  publications,  and indicate your *
0022 // * acceptance of all terms of the Geant4 Software license.          *
0023 // ********************************************************************
0024 //
0025 // G4BorisDriver
0026 //
0027 // Class description:
0028 //
0029 //   G4BorisDriver is a driver class using the second order Boris
0030 // method to integrate the equation of motion.
0031 //
0032 //
0033 // Author: Divyansh Tiwari, Google Summer of Code 2022
0034 // Supervision: John Apostolakis,Renee Fatemi, Soon Yung Jun
0035 // --------------------------------------------------------------------
0036 #ifndef G4BORIS_DRIVER_HH
0037 #define G4BORIS_DRIVER_HH
0038
0039 #include "G4VIntegrationDriver.hh"
0040 #include "G4BorisScheme.hh"
0041 #include "G4ChordFinderDelegate.hh"
0042
0043
0044 class G4BorisDriver : public G4VIntegrationDriver,
0045                       public G4ChordFinderDelegate<G4BorisDriver>
0046 {
0047   public:
0048
0049     G4BorisDriver( G4double hminimum,
0050                    G4BorisScheme* Boris,
0051                    G4int numberOfComponents = 6,
0052                    G4bool verbosity = false);
0053
0054     inline ~G4BorisDriver() override = default;
0055
0056     inline G4BorisDriver(const G4BorisDriver&) = delete;
0057     inline G4BorisDriver& operator=(const G4BorisDriver&) = delete;
0058
0059     // 1. Core methods that advance the integration
0060
0061     G4bool AccurateAdvance( G4FieldTrack& track,
0062                             G4double stepLen,
0063                             G4double epsilon,
0064                             G4double beginStep = 0) override;
0065       // Advance integration accurately
0066       // - by relative accuracy better than 'epsilon'
0067
0068     G4bool QuickAdvance(       G4FieldTrack& y_val,   // In/Out
0069                          const G4double dydx[],
0070                                G4double       hstep,
0071                                G4double&   missDist,  // Out: estimated sagitta
0072                                G4double&   dyerr   ) override;
0073       // Attempt one integration step, and return estimated error 'dyerr'
0074
0075     void OneGoodStep(G4double  yCurrentState[],  // In/Out: state ('y')
0076                      G4double& curveLength,      // In/Out: 'x'
0077                      G4double  htry,             // step to attempt
0078                      G4double  epsilon_rel,      // relative accuracy
0079                      G4double  restMass,
0080                      G4double  charge,
0081                      G4double& hdid,             // Out: step achieved
0082                      G4double& hnext);           // Out: proposed next step
0083       // Method to implement Accurate Advance
0084
0085     // 2. Methods needed to co-work with G4ChordFinder
0086
0087     G4double AdvanceChordLimited(G4FieldTrack& track,
0088                                  G4double hstep,
0089                                  G4double eps,
0090                                  G4double chordDistance) override
0091     {
0092       return ChordFinderDelegate::
0093              AdvanceChordLimitedImpl(track, hstep, eps, chordDistance);
0094     }
0095
0096     void OnStartTracking() override
0097     {
0098       ChordFinderDelegate::ResetStepEstimate();
0099     }
0100
0101     void OnComputeStep(const G4FieldTrack*) override {}
0102
0103     // 3. Does the method redo integrations when called to obtain values for
0104     //    internal, smaller intervals? (when needed to identify an intersection)
0105
0106     G4bool DoesReIntegrate() const override { return true; }
0107       // It would be no if it just used interpolation to provide a result.
0108
0109     // 4. Relevant for calculating a new step size to achieve required accuracy
0110
0111     inline G4double ComputeNewStepSize(G4double  errMaxNorm, // normalised error
0112                                        G4double  hstepCurrent) override; // current step size
0113
0114     G4double ShrinkStepSize2(G4double h, G4double error2) const;
0115     G4double GrowStepSize2(G4double h, G4double error2) const;
0116       // Calculate the next step size given the square of the relative error
0117
0118     // 5. Auxiliary Methods ...
0119
0120     void GetDerivatives( const G4FieldTrack& track,
0121                                G4double dydx[] ) const override;
0122
0123     void GetDerivatives( const G4FieldTrack& track,
0124                                G4double dydx[],
0125                                G4double field[] ) const override;
0126
0127     inline void SetVerboseLevel(G4int level) override;
0128     inline G4int GetVerboseLevel() const override;
0129
0130     inline G4EquationOfMotion* GetEquationOfMotion() override;
0131     inline const G4EquationOfMotion* GetEquationOfMotion() const;
0132     void SetEquationOfMotion(G4EquationOfMotion* equation) override;
0133
0134     void  StreamInfo( std::ostream& os ) const override;
0135      // Write out the parameters / state of the driver
0136
0137     // 6. Not relevant for Boris and other non-RK methods
0138
0139     inline const G4MagIntegratorStepper* GetStepper() const override;
0140     inline G4MagIntegratorStepper* GetStepper() override;
0141
0142   private:
0143
0144     inline G4int GetNumberOfVariables() const;
0145
0146     inline void CheckStep(const G4ThreeVector& posIn,
0147                           const G4ThreeVector& posOut,
0148                                 G4double hdid) const;
0149
0150   private:
0151
0152     // INVARIANTS -- remain unchanged during tracking / integration
0153     // Parameters
0154     G4double fMinimumStep;
0155     G4bool   fVerbosity;
0156
0157     // State -- The core stepping algorithm
0158     G4BorisScheme* boris;
0159
0160     // STATE -- intermediate state (to avoid creation / churn )
0161     G4double yIn[G4FieldTrack::ncompSVEC],
0162              yMid[G4FieldTrack::ncompSVEC],
0163              yOut[G4FieldTrack::ncompSVEC],
0164              yError[G4FieldTrack::ncompSVEC];
0165
0166     G4double yCurrent[G4FieldTrack::ncompSVEC];
0167
0168     // - Unused 2022.11.03:
0169     // G4double derivs[2][6][G4FieldTrack::ncompSVEC];
0170     // const G4int interval_sequence[2];
0171
0172     // INVARIANTS -- Parameters for ensuring that one call has finite number of integration steps
0173     static constexpr G4int    fMaxNoSteps = 300;
0174     static constexpr G4double fSmallestFraction= 1e-12; // To avoid FP underflow !  ( 1.e-6 for single prec)
0175
0176     static constexpr G4int    fIntegratorOrder= 2; //  2nd order method -- needed for error control
0177     static constexpr G4double fSafetyFactor = 0.9; //
0178
0179     static constexpr G4double fMaxSteppingIncrease= 10.0; //  Increase no more than 10x
0180     static constexpr G4double fMaxSteppingDecrease= 0.1;  //  Reduce   no more than 10x
0181     static constexpr G4double fPowerShrink = -1.0 / fIntegratorOrder;
0182     static constexpr G4double fPowerGrow   = -1.0 / (1.0 + fIntegratorOrder);
0183
0184     static const G4double fErrorConstraintShrink;
0185     static const G4double fErrorConstraintGrow;
0186
0187     using ChordFinderDelegate = G4ChordFinderDelegate<G4BorisDriver>;
0188 };
0189
0190 #include "G4BorisDriver.icc"
0191
0192 #endif