File indexing completed on 2024-04-06 12:30:26
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
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);
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
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
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;
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
0085
0086
0087
0088 G4double fLastStepLength = -1.0;
0089 T_Equation* fEquation_Rhs;
0090 };
0091
0092
0093
0094
0095
0096
0097
0098
0099
0100
0101
0102
0103
0104
0105
0106
0107
0108
0109
0110
0111
0112
0113
0114
0115
0116
0117
0118
0119 #include "G4LineSection.hh"
0120
0121 #include <cstring>
0122
0123
0124
0125
0126
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
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
0141
0142
0143
0144
0145
0146
0147
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
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
0176
0177
0178
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
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
0198
0199
0200
0201
0202
0203
0204
0205
0206
0207
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
0213
0214 field_utils::ShortState<N8> yTemp;
0215
0216 yOut[7] = yTemp[7] = fyIn[7] = yInput[7];
0217
0218
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);
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);
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);
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);
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);
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);
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
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
0280
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
0301
0302
0303
0304
0305 template <class T_Equation, unsigned int N>
0306 void CMSTDormandPrince45<T_Equation, N>::Interpolate4thOrder(G4double yOut[], G4double tau) const {
0307
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
0337
0338
0339
0340
0341
0342
0343 template <class T_Equation, unsigned int N>
0344 void CMSTDormandPrince45<T_Equation, N>::SetupInterpolation5thOrder() {
0345
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
0355 field_utils::ShortState<N> yTemp;
0356
0357
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);
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);
0372 }
0373
0374
0375
0376 template <class T_Equation, unsigned int N>
0377 void CMSTDormandPrince45<T_Equation, N>::Interpolate5thOrder(G4double yOut[], G4double tau) const {
0378
0379
0380 G4double bi[10][5];
0381
0382
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
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
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
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
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
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
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
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
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
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
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