include/Geant4/G4TCashKarpRKF45.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 // G4TCashKarpRKF45
0027 //
0028 // Class description:
0029 //
0030 // Templated version of Cash-Karp 4th/5th order embedded stepper
0031 //
0032 // Knowing the type (class) of the equation of motion enables a non-
0033 // virtual call of its methods.
0034 // As an embedded 5th order method, it requires fewer field evaluations
0035 // (1 initial + 5 others per step = 6 per step) than ClassicalRK4 and
0036 // also non-embedded methods of the same order.
0037 //
0038 // Can be used to enable use of non-virtual calls for field, equation,
0039 // and stepper - potentially with inlined methods.
0040 //
0041 // Adapted from G4CashKarpRKF45 class
0042 // --------------------------------------------------------------------
0043 // Original description (G4CashKarpRKF45):
0044 // The Cash-Karp Runge-Kutta-Fehlberg 4/5 method is an embedded fourth
0045 // order method (giving fifth-order accuracy) for the solution of an ODE.
0046 // Two different fourth order estimates are calculated; their difference
0047 // gives an error estimate. [ref. Numerical Recipes in C, 2nd Edition]
0048 // Used to integrate the equations of motion of a particle in a field.
0049
0050 // Author: Josh Xie (CERN, Google Summer of Code 2014), June 2014
0051 // Supervisors:  Sandro Wenzel, John Apostolakis (CERN)
0052 // --------------------------------------------------------------------
0053 #ifndef G4T_CASH_KARP_RKF45_HH
0054 #define G4T_CASH_KARP_RKF45_HH
0055
0056 #include <cassert>
0057
0058 #include "G4LineSection.hh"
0059 #include "G4MagIntegratorStepper.hh"
0060
0061 /**
0062  * @brief G4TCashKarpRKF45 is a templated version of Cash-Karp
0063  * 4th/5th order embedded stepper.
0064  */
0065
0066 template <class T_Equation, unsigned int N = 6 >
0067 class G4TCashKarpRKF45 : public G4MagIntegratorStepper
0068 {
0069   public:
0070
0071     G4TCashKarpRKF45(T_Equation* EqRhs, // G4int noIntegrationVariables = 6,
0072                      G4bool primary = true);
0073
0074     virtual ~G4TCashKarpRKF45();
0075
0076     G4TCashKarpRKF45(const G4TCashKarpRKF45&) = delete;
0077     G4TCashKarpRKF45& operator=(const G4TCashKarpRKF45&) = delete;
0078
0079     inline void
0080     StepWithError(const G4double yInput[], // * __restrict__ yInput,
0081                   const G4double dydx[],   // * __restrict__ dydx,
0082                   G4double Step,
0083                   G4double yOut[],         // * __restrict__ yOut,
0084                   G4double yErr[] );       // * __restrict__ yErr);
0085
0086     virtual void Stepper(const G4double yInput[],
0087                          const G4double dydx[],
0088                          G4double hstep,
0089                          G4double yOutput[],
0090                          G4double yError[]) override final;
0091
0092     // __attribute__((always_inline))
0093     void RightHandSideInl( const G4double y[],  // * __restrict__  y,
0094                                  G4double dydx[] ) // * __restrict__  dydx )
0095     {
0096       fEquation_Rhs->T_Equation::RightHandSide(y, dydx);
0097     }
0098
0099     inline G4double DistChord() const override;
0100
0101     inline G4int IntegratorOrder() const override { return 4; }
0102
0103     G4StepperType StepperType() const override { return kTCashKarpRKF45; }
0104
0105   private:
0106
0107     G4double ak2[N], ak3[N], ak4[N], ak5[N], ak6[N], ak7[N], yTemp[N], yIn[N];
0108     // scratch space
0109
0110     G4double fLastStepLength= 0.0;
0111     G4double* fLastInitialVector;
0112     G4double* fLastFinalVector;
0113     G4double* fLastDyDx;
0114     G4double* fMidVector;
0115     G4double* fMidError;
0116     // for DistChord calculations
0117
0118     G4TCashKarpRKF45* fAuxStepper = nullptr;
0119     // ... or G4TCashKarpRKF45<T_Equation, N>* fAuxStepper;
0120     T_Equation* fEquation_Rhs;
0121 };
0122
0123 /////////////////////////////////////////////////////////////////////
0124 //
0125 // Constructor
0126 //
0127 template <class T_Equation, unsigned int N >
0128 G4TCashKarpRKF45<T_Equation,N>::G4TCashKarpRKF45(T_Equation* EqRhs,
0129                                                  G4bool primary)
0130     : G4MagIntegratorStepper(dynamic_cast<G4EquationOfMotion*>(EqRhs), N )
0131     , fEquation_Rhs(EqRhs)
0132 {
0133   if( dynamic_cast<G4EquationOfMotion*>(EqRhs) == nullptr )
0134   {
0135     G4Exception("G4TCashKarpRKF45: constructor", "GeomField0001",
0136                 FatalException, "Equation is not an G4EquationOfMotion.");
0137   }
0138
0139   fLastInitialVector = new G4double[N];
0140   fLastFinalVector   = new G4double[N];
0141   fLastDyDx          = new G4double[N];
0142
0143   fMidVector = new G4double[N];
0144   fMidError  = new G4double[N];
0145
0146   if(primary)
0147   {
0148      fAuxStepper = new G4TCashKarpRKF45<T_Equation, N> (EqRhs, !primary);
0149   }
0150 }
0151
0152 template <class T_Equation, unsigned int N >
0153 G4TCashKarpRKF45<T_Equation,N>::~G4TCashKarpRKF45()
0154 {
0155   delete[] fLastInitialVector;
0156   delete[] fLastFinalVector;
0157   delete[] fLastDyDx;
0158   delete[] fMidVector;
0159   delete[] fMidError;
0160
0161   delete fAuxStepper;
0162 }
0163
0164 //////////////////////////////////////////////////////////////////////
0165 //
0166 // Given values for n = 6 variables yIn[0,...,n-1]
0167 // known  at x, use the fifth-order Cash-Karp Runge-
0168 // Kutta-Fehlberg-4-5 method to advance the solution over an interval
0169 // Step and return the incremented variables as yOut[0,...,n-1]. Also
0170 // return an estimate of the local truncation error yErr[] using the
0171 // embedded 4th-order method. The equation's method is called (inline)
0172 // via RightHandSideInl(y,dydx), which returns derivatives dydx for y .
0173 //
0174 template <class T_Equation, unsigned int N >
0175 inline void
0176 G4TCashKarpRKF45<T_Equation,N>::StepWithError(const G4double* yInput,
0177                                               const G4double* dydx,
0178                                               G4double Step,
0179                                               G4double * yOut,
0180                                               G4double * yErr)
0181 {
0182   // const G4double a2 = 0.2 , a3 = 0.3 , a4 = 0.6 , a5 = 1.0 , a6 = 0.875;
0183
0184   const G4double b21 = 0.2, b31 = 3.0 / 40.0, b32 = 9.0 / 40.0, b41 = 0.3,
0185                  b42 = -0.9, b43 = 1.2,
0186
0187                  b51 = -11.0 / 54.0, b52 = 2.5, b53 = -70.0 / 27.0,
0188                  b54 = 35.0 / 27.0,
0189
0190                  b61 = 1631.0 / 55296.0, b62 = 175.0 / 512.0,
0191                  b63 = 575.0 / 13824.0, b64 = 44275.0 / 110592.0,
0192                  b65 = 253.0 / 4096.0,
0193
0194                  c1 = 37.0 / 378.0, c3 = 250.0 / 621.0, c4 = 125.0 / 594.0,
0195                  c6 = 512.0 / 1771.0, dc5 = -277.0 / 14336.0;
0196
0197   const G4double dc1 = c1 - 2825.0 / 27648.0, dc3 = c3 - 18575.0 / 48384.0,
0198                  dc4 = c4 - 13525.0 / 55296.0, dc6 = c6 - 0.25;
0199
0200   // Initialise time to t0, needed when it is not updated by the integration.
0201   //       [ Note: Only for time dependent fields (usually electric)
0202   //                 is it neccessary to integrate the time.]
0203   // yOut[7] = yTemp[7]   = yIn[7];
0204
0205   //  Saving yInput because yInput and yOut can be aliases for same array
0206   for(unsigned int i = 0; i < N; ++i)
0207   {
0208     yIn[i] = yInput[i];
0209   }
0210   // RightHandSideInl(yIn, dydx) ;              // 1st Step
0211
0212   for(unsigned int i = 0; i < N; ++i)
0213   {
0214     yTemp[i] = yIn[i] + b21 * Step * dydx[i];
0215   }
0216   this->RightHandSideInl(yTemp, ak2);  // 2nd Step
0217
0218   for(unsigned int i = 0; i < N; ++i)
0219   {
0220     yTemp[i] = yIn[i] + Step * (b31 * dydx[i] + b32 * ak2[i]);
0221   }
0222   this->RightHandSideInl(yTemp, ak3);  // 3rd Step
0223
0224   for(unsigned int i = 0; i < N; ++i)
0225   {
0226     yTemp[i] = yIn[i] + Step * (b41 * dydx[i] + b42 * ak2[i] + b43 * ak3[i]);
0227   }
0228   this->RightHandSideInl(yTemp, ak4);  // 4th Step
0229
0230   for(unsigned int i = 0; i < N; ++i)
0231   {
0232     yTemp[i] = yIn[i] + Step * (b51 * dydx[i] + b52 * ak2[i] + b53 * ak3[i] +
0233                                 b54 * ak4[i]);
0234   }
0235   this->RightHandSideInl(yTemp, ak5);  // 5th Step
0236
0237   for(unsigned int i = 0; i < N; ++i)
0238   {
0239     yTemp[i] = yIn[i] + Step * (b61 * dydx[i] + b62 * ak2[i] + b63 * ak3[i] +
0240                                 b64 * ak4[i] + b65 * ak5[i]);
0241   }
0242   this->RightHandSideInl(yTemp, ak6);  // 6th Step
0243
0244   for(unsigned int i = 0; i < N; ++i)
0245   {
0246     // Accumulate increments with proper weights
0247
0248     yOut[i] = yIn[i] +
0249        Step * (c1 * dydx[i] + c3 * ak3[i] + c4 * ak4[i] + c6 * ak6[i]);
0250   }
0251   for(unsigned int i = 0; i < N; ++i)
0252   {
0253     // Estimate error as difference between 4th and
0254     // 5th order methods
0255
0256     yErr[i] = Step * (dc1 * dydx[i] + dc3 * ak3[i] + dc4 * ak4[i] +
0257                       dc5 * ak5[i] + dc6 * ak6[i]);
0258   }
0259   for(unsigned int i = 0; i < N; ++i)
0260   {
0261     // Store Input and Final values, for possible use in calculating chord
0262     fLastInitialVector[i] = yIn[i];
0263     fLastFinalVector[i]   = yOut[i];
0264     fLastDyDx[i]          = dydx[i];
0265   }
0266   // NormaliseTangentVector( yOut ); // Not wanted
0267
0268   fLastStepLength = Step;
0269
0270   return;
0271 }
0272
0273 template <class T_Equation, unsigned int N >
0274 inline G4double
0275 G4TCashKarpRKF45<T_Equation,N>::DistChord() const
0276 {
0277   G4double distLine, distChord;
0278   G4ThreeVector initialPoint, finalPoint, midPoint;
0279
0280   // Store last initial and final points (they will be overwritten in
0281   // self-Stepper call!)
0282   initialPoint = G4ThreeVector(fLastInitialVector[0], fLastInitialVector[1],
0283                                fLastInitialVector[2]);
0284   finalPoint   = G4ThreeVector(fLastFinalVector[0], fLastFinalVector[1],
0285                                fLastFinalVector[2]);
0286
0287   // Do half a step using StepNoErr
0288
0289   fAuxStepper->G4TCashKarpRKF45::Stepper(fLastInitialVector, fLastDyDx,
0290                                          0.5 * fLastStepLength, fMidVector,
0291                                          fMidError);
0292
0293   midPoint = G4ThreeVector(fMidVector[0], fMidVector[1], fMidVector[2]);
0294
0295   // Use stored values of Initial and Endpoint + new Midpoint to evaluate
0296   //  distance of Chord
0297
0298   if(initialPoint != finalPoint)
0299   {
0300     distLine  = G4LineSection::Distline(midPoint, initialPoint, finalPoint);
0301     distChord = distLine;
0302   }
0303   else
0304   {
0305     distChord = (midPoint - initialPoint).mag();
0306   }
0307   return distChord;
0308 }
0309
0310 template <class T_Equation, unsigned int N >
0311 inline void
0312 G4TCashKarpRKF45<T_Equation,N>::Stepper(const G4double yInput[],
0313                                         const G4double dydx[],
0314                                         G4double Step,
0315                                         G4double yOutput[],
0316                                         G4double yError[])
0317 {
0318   assert( yOutput != yInput );
0319   assert( yError  != yInput );
0320
0321   StepWithError( yInput, dydx, Step, yOutput, yError);
0322 }
0323
0324 #endif