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