OPAL (Object Oriented Parallel Accelerator Library)  2024.1
OPAL
RK4.hpp
Go to the documentation of this file.
1 //
2 // Class RK4
3 // Fourth order Runge-Kutta time integrator
4 //
5 // Copyright (c) 2008 - 2020, Paul Scherrer Institut, Villigen PSI, Switzerland
6 // All rights reserved
7 //
8 // This file is part of OPAL.
9 //
10 // OPAL is free software: you can redistribute it and/or modify
11 // it under the terms of the GNU General Public License as published by
12 // the Free Software Foundation, either version 3 of the License, or
13 // (at your option) any later version.
14 //
15 // You should have received a copy of the GNU General Public License
16 // along with OPAL. If not, see <https://www.gnu.org/licenses/>.
17 //
18 
19 template <typename FieldFunction, typename ... Arguments>
21  const size_t& i,
22  const double& t,
23  const double dt,
24  Arguments& ... args) const
25 {
26  // Fourth order Runge-Kutta integrator
27  // arguments:
28  // x Current value of dependent variable
29  // t Independent variable (usually time)
30  // dt Step size (usually time step)
31  // i index of particle
32 
33  double x[6];
34 
35  this->copyTo(bunch->R[i], bunch->P[i], &x[0]);
36 
37  double deriv1[6];
38  double deriv2[6];
39  double deriv3[6];
40  double deriv4[6];
41  double xtemp[6];
42 
43  // Evaluate f1 = f(x,t).
44  bool outOfBound = derivate_m(bunch, x, t, deriv1 , i, args ...);
45  if (outOfBound) {
46  return false;
47  }
48 
49  // Evaluate f2 = f( x+dt*f1/2, t+dt/2 ).
50  const double half_dt = 0.5 * dt;
51  const double t_half = t + half_dt;
52 
53  for(int j = 0; j < 6; ++j)
54  xtemp[j] = x[j] + half_dt * deriv1[j];
55 
56  outOfBound = derivate_m(bunch, xtemp, t_half, deriv2 , i, args ...);
57  if (outOfBound) {
58  return false;
59  }
60 
61  // Evaluate f3 = f( x+dt*f2/2, t+dt/2 ).
62  for (int j = 0; j < 6; ++j) {
63  xtemp[j] = x[j] + half_dt * deriv2[j];
64  }
65 
66  outOfBound = derivate_m(bunch, xtemp, t_half, deriv3 , i, args ...);
67  if (outOfBound) {
68  return false;
69  }
70 
71  // Evaluate f4 = f( x+dt*f3, t+dt ).
72  double t_full = t + dt;
73  for (int j = 0; j < 6; ++j) {
74  xtemp[j] = x[j] + dt * deriv3[j];
75  }
76 
77  outOfBound = derivate_m(bunch, xtemp, t_full, deriv4 , i, args ...);
78  if (outOfBound) {
79  return false;
80  }
81 
82  // Return x(t+dt) computed from fourth-order R-K.
83  for (int j = 0; j < 6; ++j) {
84  x[j] += dt / 6.*(deriv1[j] + deriv4[j] + 2.*(deriv2[j] + deriv3[j]));
85  }
86  this->copyFrom(bunch->R[i], bunch->P[i], &x[0]);
87 
88  return true;
89 }
90 
91 
92 template <typename FieldFunction, typename ... Arguments>
94  double* y,
95  const double& t,
96  double* yp,
97  const size_t& i,
98  Arguments& ... args) const
99 {
100  // New for OPAL 2.0: Changing variables to m, T, s
101  // Currently: m, ns, kG
102 
103  Vector_t externalE, externalB, tempR;
104 
105  externalB = Vector_t(0.0, 0.0, 0.0);
106  externalE = Vector_t(0.0, 0.0, 0.0);
107 
108  for (int j = 0; j < 3; ++j) {
109  tempR(j) = y[j];
110  }
111  bunch->R[i] = tempR;
112 
113  bool outOfBound = this->fieldfunc_m(t, i, externalE, externalB, args ...);
114 
115  double qtom = bunch->Q[i] / (bunch->M[i] * mass_coeff); // m^2/s^2/GV
116 
117  double tempgamma = std::sqrt(1 + (y[3] * y[3] + y[4] * y[4] + y[5] * y[5]));
118 
119  yp[0] = Physics::c / tempgamma * y[3];
120  yp[1] = Physics::c / tempgamma * y[4];
121  yp[2] = Physics::c / tempgamma * y[5];
122 
123  yp[3] = (externalE(0) / Physics::c + (externalB(2) * y[4] - externalB(1) * y[5]) / tempgamma) * qtom; // [1/ns]
124  yp[4] = (externalE(1) / Physics::c - (externalB(2) * y[3] - externalB(0) * y[5]) / tempgamma) * qtom; // [1/ns];
125  yp[5] = (externalE(2) / Physics::c + (externalB(1) * y[3] - externalB(0) * y[4]) / tempgamma) * qtom; // [1/ns];
126 
127  yp[3] /= Units::ns2s;
128  yp[4] /= Units::ns2s;
129  yp[5] /= Units::ns2s;
130 
131  return outOfBound;
132 }
133 
134 
135 template <typename FieldFunction, typename ... Arguments>
137  const Vector_t& P,
138  double* x) const
139 {
140  for (int j = 0; j < 3; j++) {
141  x[j] = R(j); // [x,y,z] (mm)
142  x[j+3] = P(j); // [px,py,pz] (beta*gamma)
143  }
144 }
145 
146 
147 template <typename FieldFunction, typename ... Arguments>
149  Vector_t& P,
150  double* x) const
151 {
152  for (int j = 0; j < 3; j++) {
153  R(j) = x[j]; // [x,y,z] (mm)
154  P(j) = x[j+3]; // [px,py,pz] (beta*gamma)
155  }
156 }
Tps< T > sqrt(const Tps< T > &x)
Square root.
Definition: TpsMath.h:91
constexpr double c
The velocity of light in m/s.
Definition: Physics.h:45
ParticleAttrib< Vector_t > P
Vektor< double, 3 > Vector_t
Definition: Vektor.h:6
void copyFrom(Vector_t &R, Vector_t &P, double *x) const
Definition: RK4.hpp:148
ParticleAttrib< double > M
void copyTo(const Vector_t &R, const Vector_t &P, double *x) const
Definition: RK4.hpp:136
constexpr double ns2s
Definition: Units.h:47
ParticleAttrib< double > Q
ParticlePos_t & R
bool derivate_m(PartBunchBase< double, 3 > *bunch, double *y, const double &t, double *yp, const size_t &i, Arguments &...args) const
Definition: RK4.hpp:93
bool doAdvance_m(PartBunchBase< double, 3 > *bunch, const size_t &i, const double &t, const double dt, Arguments &...args) const
Definition: RK4.hpp:20