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