include/Geant4/G4TClassicalRK4.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 // G4TClassicalRK4
0027 //
0028 // Class description:
0029 //
0030 // Templated version of G4ClassicalRK4
0031 //
0032 //
0033 // Created: Josh Xie  (supported by Google Summer of Code 2014 )
0034 // Supervisors:  Sandro Wenzel, John Apostolakis (CERN)
0035 // Adapted from G4G4TClassicalRK4 class
0036 // --------------------------------------------------------------------
0037 #include "G4ThreeVector.hh"
0038 #include "G4MagIntegratorStepper.hh"
0039 #include "G4TMagErrorStepper.hh"
0040
0041 template <class T_Equation, unsigned int N>
0042 class G4TClassicalRK4
0043   : public G4TMagErrorStepper<G4TClassicalRK4<T_Equation, N>, T_Equation, N>
0044 {
0045  public:  // with description
0046   static constexpr G4double IntegratorCorrection = 1. / ((1 << 4) - 1);
0047
0048   G4TClassicalRK4(T_Equation* EqRhs, G4int numberOfVariables = 8);
0049
0050   virtual ~G4TClassicalRK4() { ; }
0051
0052   void RightHandSideInl(G4double y[],
0053                      G4double dydx[])
0054   {
0055     fEquation_Rhs->T_Equation::RightHandSide(y, dydx);
0056   }
0057
0058   // A stepper that does not know about errors.
0059   // It is used by the MagErrorStepper stepper.
0060
0061   inline  // __attribute__((always_inline))
0062   void DumbStepper(const G4double yIn[],
0063                    const G4double dydx[],
0064                    G4double h,
0065                    G4double yOut[]);
0066
0067  public:  // without description
0068   G4int IntegratorOrder() const { return 4; }
0069
0070   G4TClassicalRK4(const G4TClassicalRK4&) = delete;
0071   G4TClassicalRK4& operator=(const G4TClassicalRK4&) = delete;
0072   // No copy constructor and assignment operator.
0073
0074  private:
0075   // G4int fNumberOfVariables ; // is set default to 6 in constructor
0076   G4double dydxm[N < 8 ? 8 : N];
0077   G4double dydxt[N < 8 ? 8 : N];
0078   G4double yt[N < 8 ? 8 : N];
0079   // scratch space - not state
0080
0081   T_Equation* fEquation_Rhs;
0082 };
0083
0084 template <class T_Equation, unsigned int N >
0085 G4TClassicalRK4<T_Equation,N>::
0086 G4TClassicalRK4(T_Equation* EqRhs, G4int numberOfVariables)
0087    : G4TMagErrorStepper<G4TClassicalRK4<T_Equation, N>, T_Equation, N>(
0088       EqRhs, numberOfVariables > 8 ? numberOfVariables : 8 )
0089    , fEquation_Rhs(EqRhs)
0090 {
0091   // unsigned int noVariables = std::max(numberOfVariables, 8);  // For Time .. 7+1
0092   if( dynamic_cast<G4EquationOfMotion*>(EqRhs) == nullptr )
0093   {
0094     G4Exception("G4TClassicalRK4: constructor", "GeomField0001",
0095                 FatalException, "Equation is not an G4EquationOfMotion.");
0096   }
0097 }
0098
0099 template <class T_Equation, unsigned int N >
0100 void
0101 G4TClassicalRK4<T_Equation,N>::DumbStepper(const G4double yIn[],
0102                                            const G4double dydx[],
0103                                            G4double h,
0104                                            G4double yOut[])
0105 // Given values for the variables y[0,..,n-1] and their derivatives
0106 // dydx[0,...,n-1] known at x, use the classical 4th Runge-Kutta
0107 // method to advance the solution over an interval h and return the
0108 // incremented variables as yout[0,...,n-1], which not be a distinct
0109 // array from y. The user supplies the routine RightHandSide(x,y,dydx),
0110 // which returns derivatives dydx at x. The source is routine rk4 from
0111 // NRC p. 712-713 .
0112 {
0113   G4double hh = h * 0.5, h6 = h / 6.0;
0114
0115   // Initialise time to t0, needed when it is not updated by the integration.
0116   //        [ Note: Only for time dependent fields (usually electric)
0117   //                  is it neccessary to integrate the time.]
0118   yt[7]   = yIn[7];
0119   yOut[7] = yIn[7];
0120
0121   for(unsigned int i = 0; i < N; ++i)
0122   {
0123      yt[i] = yIn[i] + hh * dydx[i];  // 1st Step K1=h*dydx
0124   }
0125   this->RightHandSideInl(yt, dydxt);  // 2nd Step K2=h*dydxt
0126
0127   for(unsigned int i = 0; i < N; ++i)
0128   {
0129      yt[i] = yIn[i] + hh * dydxt[i];
0130   }
0131   this->RightHandSideInl(yt, dydxm);  // 3rd Step K3=h*dydxm
0132
0133   for(unsigned int i = 0; i < N; ++i)
0134   {
0135      yt[i] = yIn[i] + h * dydxm[i];
0136      dydxm[i] += dydxt[i];  // now dydxm=(K2+K3)/h
0137   }
0138   this->RightHandSideInl(yt, dydxt);  // 4th Step K4=h*dydxt
0139
0140   for(unsigned int i = 0; i < N; ++i)  // Final RK4 output
0141   {
0142      yOut[i] = yIn[i] + h6 * (dydx[i] + dydxt[i] +
0143                               2.0 * dydxm[i]);  //+K1/6+K4/6+(K2+K3)/3
0144   }
0145   if(N == 12)
0146   {
0147      this->NormalisePolarizationVector(yOut);
0148   }
0149
0150 }  // end of DumbStepper ....................................................
0151
0152 // template <class T_Equation, unsigned int N >
0153 // G4TClassicalRK4<T_Equation,N>::
0154