# Project CMSSW displayed by LXR

File indexing completed on 2024-04-06 12:30:26

```0001 //
0002 // copied from G4TDormandPrince45
0003 //
0004 // Class desription:
0005 //
0006 //  An implementation of the 5th order embedded RK method from the paper:
0007 //  J. R. Dormand and P. J. Prince, "A family of embedded Runge-Kutta formulae"
0008 //  Journal of computational and applied Math., vol.6, no.1, pp.19-26, 1980.
0009 //
0010 //  DormandPrince7 - 5(4) embedded RK method
0011 //
0012 //
0013 // Created: Somnath Banerjee, Google Summer of Code 2015, 25 May 2015
0014 // Supervision: John Apostolakis, CERN
0015 //
0016 //  Modified by Vincenzo Innocente on 7/2022 to account for constant momentum magnitude
0017 //
0018 // --------------------------------------------------------------------
0019 #ifndef CMSTDORMAND_PRINCE_45_HH
0020 #define CMSTDORMAND_PRINCE_45_HH
0021
0022 #include <cassert>
0023 #include "G4MagIntegratorStepper.hh"
0024 #include "G4FieldUtils.hh"
0025
0026 template <class T_Equation, unsigned int N = 6>
0027 class CMSTDormandPrince45 : public G4MagIntegratorStepper {
0028 public:
0029   CMSTDormandPrince45(T_Equation* equation);
0030   CMSTDormandPrince45(T_Equation* equation, G4int numVar);  // must have numVar == N
0031
0032   inline void StepWithError(
0033       const G4double yInput[], const G4double dydx[], G4double hstep, G4double yOutput[], G4double yError[]);
0034
0035   void Stepper(
0036       const G4double yInput[], const G4double dydx[], G4double hstep, G4double yOutput[], G4double yError[]) final;
0037
0038   inline void StepWithFinalDerivate(const G4double yInput[],
0039                                     const G4double dydx[],
0040                                     G4double hstep,
0041                                     G4double yOutput[],
0042                                     G4double yError[],
0043                                     G4double dydxOutput[]);
0044
0045   inline void SetupInterpolation() {}
0046
0047   void Interpolate(G4double tau, G4double yOut[]) const { Interpolate4thOrder(yOut, tau); }
0048   // For calculating the output at the tau fraction of Step
0049
0050   G4double DistChord() const final;
0051
0052   G4int IntegratorOrder() const override { return 4; }
0053
0054   const field_utils::ShortState<N>& GetYOut() const { return fyOut; }
0055
0056   void Interpolate4thOrder(G4double yOut[], G4double tau) const;
0057
0058   void SetupInterpolation5thOrder();
0059   void Interpolate5thOrder(G4double yOut[], G4double tau) const;
0060
0061   // __attribute__((always_inline))
0062   void RightHandSideInl(const G4double y[], G4double inv_momentum_magnitude, G4double dydx[]) {
0063     fEquation_Rhs->T_Equation::TRightHandSide(y, inv_momentum_magnitude, dydx);
0064   }
0065
0066   inline void Stepper(const G4double yInput[],
0067                       const G4double dydx[],
0068                       G4double hstep,
0069                       G4double yOutput[],
0070                       G4double yError[],
0071                       G4double dydxOutput[]) {
0072     StepWithFinalDerivate(yInput, dydx, hstep, yOutput, yError, dydxOutput);
0073   }
0074
0075   T_Equation* GetSpecificEquation() { return fEquation_Rhs; }
0076
0077   static constexpr int N8 = N > 8 ? N : 8;  //  y[
0078
0079 private:
0080   field_utils::ShortState<N> ak2, ak3, ak4, ak5, ak6, ak7, ak8, ak9;
0081   field_utils::ShortState<N8> fyIn;
0082   field_utils::ShortState<N> fyOut, fdydxIn;
0083
0084   // - Simpler :
0085   // field_utils::State ak2, ak3, ak4, ak5, ak6, ak7, ak8, ak9;
0086   // field_utils::State fyIn, fyOut, fdydxIn;
0087
0088   G4double fLastStepLength = -1.0;
0089   T_Equation* fEquation_Rhs;
0090 };
0091
0092 // G4TDormandPrince745 implementation -- borrowed from G4DormandPrince745
0093 //
0094 // DormandPrince7 - 5(4) non-FSAL
0095 // definition of the stepper() method that evaluates one step in
0096 // field propagation.
0097 // The coefficients and the algorithm have been adapted from
0098 //
0099 // J. R. Dormand and P. J. Prince, "A family of embedded Runge-Kutta formulae"
0100 // Journal of computational and applied Math., vol.6, no.1, pp.19-26, 1980.
0101 //
0102 // The Butcher table of the Dormand-Prince-7-4-5 method is as follows :
0103 //
0104 //    0   |
0105 //    1/5 | 1/5
0106 //    3/10| 3/40       9/40
0107 //    4/5 | 44/45      56/15      32/9
0108 //    8/9 | 19372/6561 25360/2187 64448/6561  212/729
0109 //    1   | 9017/3168  355/33     46732/5247  49/176  5103/18656
0110 //    1   | 35/384     0          500/1113    125/192 2187/6784    11/84
0111 //    ------------------------------------------------------------------------
0112 //          35/384     0          500/1113    125/192 2187/6784    11/84    0
0113 //          5179/57600 0          7571/16695  393/640 92097/339200 187/2100 1/40
0114 //
0115 // Created: Somnath Banerjee, Google Summer of Code 2015, 25 May 2015
0116 // Supervision: John Apostolakis, CERN
0117 // --------------------------------------------------------------------
0118
0119 #include "G4LineSection.hh"
0120
0121 #include <cstring>
0122
0123 // using namespace field_utils;
0124
0125 /////////////////////////////////////////////////////////////////////
0126 // Constructor
0127 //
0128 template <class T_Equation, unsigned int N>
0129 CMSTDormandPrince45<T_Equation, N>::CMSTDormandPrince45(T_Equation* equation)
0130     : G4MagIntegratorStepper(dynamic_cast<G4EquationOfMotion*>(equation), N), fEquation_Rhs(equation) {
0131   // assert( dynamic_cast<G4EquationOfMotion*>(equation) != nullptr );
0132   if (dynamic_cast<G4EquationOfMotion*>(equation) == nullptr) {
0133     G4Exception("G4TDormandPrince745CMS: constructor",
0134                 "GeomField0001",
0135                 FatalException,
0136                 "T_Equation is not an G4EquationOfMotion.");
0137   }
0138
0139   /***
0140   assert( equation->GetNumberOfVariables == N );
0141   if( equation->GetNumberOfVariables != N ){
0142     G4ExceptionDescription msg;
0143     msg << "Equation has an incompatible number of variables." ;
0144     msg << "   template N = " << N << " equation-Nvar= "
0145         << equation->GetNumberOfVariables;
0146     G4Exception("G4TCashKarpRKF45: constructor", "GeomField0001",
0147                 FatalException, msg );
0148    } ****/
0149 }
0150
0151 template <class T_Equation, unsigned int N>
0152 CMSTDormandPrince45<T_Equation, N>::CMSTDormandPrince45(T_Equation* equation, G4int numVar)
0153     : CMSTDormandPrince45<T_Equation, N>(equation) {
0154   if (numVar != G4int(N)) {
0155     G4ExceptionDescription msg;
0156     msg << "Equation has an incompatible number of variables.";
0157     msg << "   template N = " << N << "   argument numVar = " << numVar;
0158     //    << " equation-Nvar= " << equation->GetNumberOfVariables(); // --> Expected later
0159     G4Exception("G4TCashKarpRKF45CMS: constructor", "GeomField0001", FatalErrorInArgument, msg);
0160   }
0161   assert(numVar == N);
0162 }
0163
0164 template <class T_Equation, unsigned int N>
0165 inline void CMSTDormandPrince45<T_Equation, N>::StepWithFinalDerivate(const G4double yInput[],
0166                                                                       const G4double dydx[],
0167                                                                       G4double hstep,
0168                                                                       G4double yOutput[],
0169                                                                       G4double yError[],
0170                                                                       G4double dydxOutput[]) {
0171   StepWithError(yInput, dydx, hstep, yOutput, yError);
0172   field_utils::copy(dydxOutput, ak7, N);
0173 }
0174
0175 // Stepper
0176 //
0177 // Passing in the value of yInput[],the first time dydx[] and Step length
0178 // Giving back yOut and yErr arrays for output and error respectively
0179 //
0180
0181 template <class T_Equation, unsigned int N>
0182 inline void CMSTDormandPrince45<T_Equation, N>::StepWithError(
0183     const G4double yInput[], const G4double dydx[], G4double hstep, G4double yOut[], G4double yErr[]) {
0184   // The parameters of the Butcher tableu
0185   //
0186   constexpr G4double b21 = 0.2, b31 = 3.0 / 40.0, b32 = 9.0 / 40.0, b41 = 44.0 / 45.0, b42 = -56.0 / 15.0,
0187                      b43 = 32.0 / 9.0,
0188
0189                      b51 = 19372.0 / 6561.0, b52 = -25360.0 / 2187.0, b53 = 64448.0 / 6561.0, b54 = -212.0 / 729.0,
0190
0191                      b61 = 9017.0 / 3168.0, b62 = -355.0 / 33.0, b63 = 46732.0 / 5247.0, b64 = 49.0 / 176.0,
0192                      b65 = -5103.0 / 18656.0,
0193
0194                      b71 = 35.0 / 384.0, b72 = 0., b73 = 500.0 / 1113.0, b74 = 125.0 / 192.0, b75 = -2187.0 / 6784.0,
0195                      b76 = 11.0 / 84.0,
0196
0197                      //Sum of columns, sum(bij) = ei
0198       //    e1 = 0. ,
0199       //    e2 = 1.0/5.0 ,
0200       //    e3 = 3.0/10.0 ,
0201       //    e4 = 4.0/5.0 ,
0202       //    e5 = 8.0/9.0 ,
0203       //    e6 = 1.0 ,
0204       //    e7 = 1.0 ,
0205
0206       // Difference between the higher and the lower order method coeff. :
0207       // b7j are the coefficients of higher order
0208
0209       dc1 = -(b71 - 5179.0 / 57600.0), dc2 = -(b72 - .0), dc3 = -(b73 - 7571.0 / 16695.0), dc4 = -(b74 - 393.0 / 640.0),
0210                      dc5 = -(b75 + 92097.0 / 339200.0), dc6 = -(b76 - 187.0 / 2100.0), dc7 = -(-1.0 / 40.0);
0211
0212   // const G4int numberOfVariables = GetNumberOfVariables();
0213   //   The number of variables to be integrated over
0214   field_utils::ShortState<N8> yTemp;
0215
0216   yOut[7] = yTemp[7] = fyIn[7] = yInput[7];  // Pass along the time - used in RightHandSide
0217
0218   //  Saving yInput because yInput and yOut can be aliases for same array
0219   //
0220   for (unsigned int i = 0; i < N; ++i) {
0221     fyIn[i] = yInput[i];
0222     yTemp[i] = yInput[i] + b21 * hstep * dydx[i];
0223   }
0224   G4double momentum_mag_square = yTemp[3] * yTemp[3] + yTemp[4] * yTemp[4] + yTemp[5] * yTemp[5];
0225   G4double inv_momentum_magnitude = 1.0 / std::sqrt(momentum_mag_square);
0226   RightHandSideInl(yTemp, inv_momentum_magnitude, ak2);  // 2nd stage
0227
0228   for (unsigned int i = 0; i < N; ++i) {
0229     yTemp[i] = fyIn[i] + hstep * (b31 * dydx[i] + b32 * ak2[i]);
0230   }
0231   RightHandSideInl(yTemp, inv_momentum_magnitude, ak3);  // 3rd stage
0232
0233   for (unsigned int i = 0; i < N; ++i) {
0234     yTemp[i] = fyIn[i] + hstep * (b41 * dydx[i] + b42 * ak2[i] + b43 * ak3[i]);
0235   }
0236   RightHandSideInl(yTemp, inv_momentum_magnitude, ak4);  // 4th stage
0237
0238   for (unsigned int i = 0; i < N; ++i) {
0239     yTemp[i] = fyIn[i] + hstep * (b51 * dydx[i] + b52 * ak2[i] + b53 * ak3[i] + b54 * ak4[i]);
0240   }
0241   RightHandSideInl(yTemp, inv_momentum_magnitude, ak5);  // 5th stage
0242
0243   for (unsigned int i = 0; i < N; ++i) {
0244     yTemp[i] = fyIn[i] + hstep * (b61 * dydx[i] + b62 * ak2[i] + b63 * ak3[i] + b64 * ak4[i] + b65 * ak5[i]);
0245   }
0246   RightHandSideInl(yTemp, inv_momentum_magnitude, ak6);  // 6th stage
0247
0248   for (unsigned int i = 0; i < N; ++i) {
0249     yOut[i] =
0250         fyIn[i] + hstep * (b71 * dydx[i] + b72 * ak2[i] + b73 * ak3[i] + b74 * ak4[i] + b75 * ak5[i] + b76 * ak6[i]);
0251   }
0252   RightHandSideInl(yOut, inv_momentum_magnitude, ak7);  // 7th and Final stage
0253
0254   for (unsigned int i = 0; i < N; ++i) {
0255     yErr[i] = hstep * (dc1 * dydx[i] + dc2 * ak2[i] + dc3 * ak3[i] + dc4 * ak4[i] + dc5 * ak5[i] + dc6 * ak6[i] +
0256                        dc7 * ak7[i]) +
0257               1.5e-18;
0258
0259     // Store Input and Final values, for possible use in calculating chord
0260     //
0261     fyOut[i] = yOut[i];
0262     fdydxIn[i] = dydx[i];
0263   }
0264
0265   fLastStepLength = hstep;
0266 }
0267
0268 template <class T_Equation, unsigned int N>
0269 inline void CMSTDormandPrince45<T_Equation, N>::Stepper(
0270     const G4double yInput[], const G4double dydx[], G4double Step, G4double yOutput[], G4double yError[]) {
0271   assert(yOutput != yInput);
0272   assert(yError != yInput);
0273
0274   StepWithError(yInput, dydx, Step, yOutput, yError);
0275 }
0276
0277 template <class T_Equation, unsigned int N>
0278 G4double CMSTDormandPrince45<T_Equation, N>::DistChord() const {
0279   // Coefficients were taken from Some Practical Runge-Kutta Formulas
0280   // by Lawrence F. Shampine, page 149, c*
0281   //
0282   const G4double hf1 = 6025192743.0 / 30085553152.0, hf3 = 51252292925.0 / 65400821598.0,
0283                  hf4 = -2691868925.0 / 45128329728.0, hf5 = 187940372067.0 / 1594534317056.0,
0284                  hf6 = -1776094331.0 / 19743644256.0, hf7 = 11237099.0 / 235043384.0;
0285
0286   G4ThreeVector mid;
0287
0288   for (unsigned int i = 0; i < 3; ++i) {
0289     mid[i] =
0290         fyIn[i] + 0.5 * fLastStepLength *
0291                       (hf1 * fdydxIn[i] + hf3 * ak3[i] + hf4 * ak4[i] + hf5 * ak5[i] + hf6 * ak6[i] + hf7 * ak7[i]);
0292   }
0293
0294   const G4ThreeVector begin = makeVector(fyIn, field_utils::Value3D::Position);
0295   const G4ThreeVector end = makeVector(fyOut, field_utils::Value3D::Position);
0296
0297   return G4LineSection::Distline(mid, begin, end);
0298 }
0299
0300 // The lower (4th) order interpolant given by Dormand and Prince:
0301 //        J. R. Dormand and P. J. Prince, "Runge-Kutta triples"
0302 //        Computers & Mathematics with Applications, vol. 12, no. 9,
0303 //        pp. 1007-1017, 1986.
0304 //
0305 template <class T_Equation, unsigned int N>
0306 void CMSTDormandPrince45<T_Equation, N>::Interpolate4thOrder(G4double yOut[], G4double tau) const {
0307   // const G4int numberOfVariables = GetNumberOfVariables();
0308
0309   const G4double tau2 = tau * tau, tau3 = tau * tau2, tau4 = tau2 * tau2;
0310
0311   const G4double bf1 =
0312       1.0 / 11282082432.0 *
0313       (157015080.0 * tau4 - 13107642775.0 * tau3 + 34969693132.0 * tau2 - 32272833064.0 * tau + 11282082432.0);
0314
0315   const G4double bf3 =
0316       -100.0 / 32700410799.0 * tau * (15701508.0 * tau3 - 914128567.0 * tau2 + 2074956840.0 * tau - 1323431896.0);
0317
0318   const G4double bf4 =
0319       25.0 / 5641041216.0 * tau * (94209048.0 * tau3 - 1518414297.0 * tau2 + 2460397220.0 * tau - 889289856.0);
0320
0321   const G4double bf5 =
0322       -2187.0 / 199316789632.0 * tau * (52338360.0 * tau3 - 451824525.0 * tau2 + 687873124.0 * tau - 259006536.0);
0323
0324   const G4double bf6 =
0325       11.0 / 2467955532.0 * tau * (106151040.0 * tau3 - 661884105.0 * tau2 + 946554244.0 * tau - 361440756.0);
0326
0327   const G4double bf7 = 1.0 / 29380423.0 * tau * (1.0 - tau) * (8293050.0 * tau2 - 82437520.0 * tau + 44764047.0);
0328
0329   for (unsigned int i = 0; i < N; ++i) {
0330     yOut[i] =
0331         fyIn[i] + fLastStepLength * tau *
0332                       (bf1 * fdydxIn[i] + bf3 * ak3[i] + bf4 * ak4[i] + bf5 * ak5[i] + bf6 * ak6[i] + bf7 * ak7[i]);
0333   }
0334 }
0335
0336 // Following interpolant of order 5 was given by Baker,Dormand,Gilmore, Prince :
0337 //        T. S. Baker, J. R. Dormand, J. P. Gilmore, and P. J. Prince,
0338 //        "Continuous approximation with embedded Runge-Kutta methods"
0339 //        Applied Numerical Mathematics, vol. 22, no. 1, pp. 51-62, 1996.
0340 //
0341 // Calculating the extra stages for the interpolant
0342 //
0343 template <class T_Equation, unsigned int N>
0344 void CMSTDormandPrince45<T_Equation, N>::SetupInterpolation5thOrder() {
0345   // Coefficients for the additional stages
0346   //
0347   const G4double b81 = 6245.0 / 62208.0, b82 = 0.0, b83 = 8875.0 / 103032.0, b84 = -125.0 / 1728.0,
0348                  b85 = 801.0 / 13568.0, b86 = -13519.0 / 368064.0, b87 = 11105.0 / 368064.0,
0349
0350                  b91 = 632855.0 / 4478976.0, b92 = 0.0, b93 = 4146875.0 / 6491016.0, b94 = 5490625.0 / 14183424.0,
0351                  b95 = -15975.0 / 108544.0, b96 = 8295925.0 / 220286304.0, b97 = -1779595.0 / 62938944.0,
0352                  b98 = -805.0 / 4104.0;
0353
0354   // const G4int numberOfVariables = GetNumberOfVariables();
0355   field_utils::ShortState<N> yTemp;
0356
0357   // Evaluate the extra stages
0358   //
0359   for (unsigned int i = 0; i < N; ++i) {
0360     yTemp[i] = fyIn[i] + fLastStepLength * (b81 * fdydxIn[i] + b82 * ak2[i] + b83 * ak3[i] + b84 * ak4[i] +
0361                                             b85 * ak5[i] + b86 * ak6[i] + b87 * ak7[i]);
0362   }
0363   G4double momentum_mag_square = yTemp[3] * yTemp[3] + yTemp[4] * yTemp[4] + yTemp[5] * yTemp[5];
0364   G4double inv_momentum_magnitude = 1.0 / std::sqrt(momentum_mag_square);
0365   RightHandSideInl(yTemp, inv_momentum_magnitude, ak8);  // 8th Stage
0366
0367   for (unsigned int i = 0; i < N; ++i) {
0368     yTemp[i] = fyIn[i] + fLastStepLength * (b91 * fdydxIn[i] + b92 * ak2[i] + b93 * ak3[i] + b94 * ak4[i] +
0369                                             b95 * ak5[i] + b96 * ak6[i] + b97 * ak7[i] + b98 * ak8[i]);
0370   }
0371   RightHandSideInl(yTemp, inv_momentum_magnitude, ak9);  // 9th Stage
0372 }
0373
0374 // Calculating the interpolated result yOut with the coefficients
0375 //
0376 template <class T_Equation, unsigned int N>
0377 void CMSTDormandPrince45<T_Equation, N>::Interpolate5thOrder(G4double yOut[], G4double tau) const {
0378   // Define the coefficients for the polynomials
0379   //
0380   G4double bi[10][5];
0381
0382   //  COEFFICIENTS OF   bi[1]
0383   bi[1][0] = 1.0, bi[1][1] = -38039.0 / 7040.0, bi[1][2] = 125923.0 / 10560.0, bi[1][3] = -19683.0 / 1760.0,
0384   bi[1][4] = 3303.0 / 880.0,
0385   //  --------------------------------------------------------
0386       //
0387       //  COEFFICIENTS OF  bi[2]
0388       bi[2][0] = 0.0, bi[2][1] = 0.0, bi[2][2] = 0.0, bi[2][3] = 0.0, bi[2][4] = 0.0,
0389   //  --------------------------------------------------------
0390       //
0391       //  COEFFICIENTS OF  bi[3]
0392       bi[3][0] = 0.0, bi[3][1] = -12500.0 / 4081.0, bi[3][2] = 205000.0 / 12243.0, bi[3][3] = -90000.0 / 4081.0,
0393   bi[3][4] = 36000.0 / 4081.0,
0394   //  --------------------------------------------------------
0395       //
0396       //  COEFFICIENTS OF  bi[4]
0397       bi[4][0] = 0.0, bi[4][1] = -3125.0 / 704.0, bi[4][2] = 25625.0 / 1056.0, bi[4][3] = -5625.0 / 176.0,
0398   bi[4][4] = 1125.0 / 88.0,
0399   //  --------------------------------------------------------
0400       //
0401       //  COEFFICIENTS OF  bi[5]
0402       bi[5][0] = 0.0, bi[5][1] = 164025.0 / 74624.0, bi[5][2] = -448335.0 / 37312.0, bi[5][3] = 295245.0 / 18656.0,
0403   bi[5][4] = -59049.0 / 9328.0,
0404   //  --------------------------------------------------------
0405       //
0406       //  COEFFICIENTS OF  bi[6]
0407       bi[6][0] = 0.0, bi[6][1] = -25.0 / 28.0, bi[6][2] = 205.0 / 42.0, bi[6][3] = -45.0 / 7.0, bi[6][4] = 18.0 / 7.0,
0408   //  --------------------------------------------------------
0409       //
0410       //  COEFFICIENTS OF  bi[7]
0411       bi[7][0] = 0.0, bi[7][1] = -2.0 / 11.0, bi[7][2] = 73.0 / 55.0, bi[7][3] = -171.0 / 55.0, bi[7][4] = 108.0 / 55.0,
0412   //  --------------------------------------------------------
0413       //
0414       //  COEFFICIENTS OF  bi[8]
0415       bi[8][0] = 0.0, bi[8][1] = 189.0 / 22.0, bi[8][2] = -1593.0 / 55.0, bi[8][3] = 3537.0 / 110.0,
0416   bi[8][4] = -648.0 / 55.0,
0417   //  --------------------------------------------------------
0418       //
0419       //  COEFFICIENTS OF  bi[9]
0420       bi[9][0] = 0.0, bi[9][1] = 351.0 / 110.0, bi[9][2] = -999.0 / 55.0, bi[9][3] = 2943.0 / 110.0,
0421   bi[9][4] = -648.0 / 55.0;
0422   //  --------------------------------------------------------
0423
0424   // Calculating the polynomials
0425
0426   G4double b[10];
0427   std::memset(b, 0.0, sizeof(b));
0428
0429   G4double tauPower = 1.0;
0430   for (G4int j = 0; j <= 4; ++j) {
0431     for (G4int iStage = 1; iStage <= 9; ++iStage) {
0432       b[iStage] += bi[iStage][j] * tauPower;
0433     }
0434     tauPower *= tau;
0435   }
0436
0437   // const G4int numberOfVariables = GetNumberOfVariables();
0438   const G4double stepLen = fLastStepLength * tau;
0439   for (G4int i = 0; i < N; ++i) {
0440     yOut[i] = fyIn[i] + stepLen * (b[1] * fdydxIn[i] + b[2] * ak2[i] + b[3] * ak3[i] + b[4] * ak4[i] + b[5] * ak5[i] +
0441                                    b[6] * ak6[i] + b[7] * ak7[i] + b[8] * ak8[i] + b[9] * ak9[i]);
0442   }
0443 }
0444
0445 #endif```