Hipace
UpdateMomentumPerez.H
Go to the documentation of this file.
1 /* Copyright 2019 Yinjian Zhao
2  *
3  * This file is part of WarpX.
4  *
5  * License: BSD-3-Clause-LBNL
6  */
7 
8 #ifndef HIPACE_UPDATE_MOMENTUM_PEREZ_ELASTIC_H_
9 #define HIPACE_UPDATE_MOMENTUM_PEREZ_ELASTIC_H_
10 
11 #include "utils/Constants.H"
12 
13 #include <AMReX_Random.H>
14 
15 #include <cmath> // isnan() isinf()
16 #include <limits> // numeric_limits<float>::min()
17 
18 /* \brief Update particle velocities according to
19  * F. Perez et al., Phys.Plasmas.19.083104 (2012),
20  * which is based on Nanbu's method, PhysRevE.55.4642 (1997).
21  * @param[in] LmdD is max(Debye length, minimal interparticle distance).
22  * @param[in] L is the Coulomb log. A fixed L will be used if L > 0,
23  * otherwise L will be calculated based on the algorithm.
24  * To see if there are nan or inf updated velocities,
25  * compile with USE_ASSERTION=TRUE.
26 */
27 
28 template <typename T_R>
31  T_R& u1x, T_R& u1y, T_R& u1z, T_R const g1, T_R& u2x, T_R& u2y, T_R& u2z, T_R const g2,
32  T_R const n1, T_R const n2, T_R const n12,
33  T_R const q1, T_R m1, T_R const w1,
34  T_R const q2, T_R m2, T_R const w2,
35  T_R const dt, T_R const L, T_R const lmdD,
36  T_R const inv_c_SI, T_R const inv_c2_SI, const bool normalized_units,
37  amrex::RandomEngine const& engine)
38 {
39  using namespace amrex::literals;
40 
41  T_R const diffx = amrex::Math::abs(u1x-u2x);
42  T_R const diffy = amrex::Math::abs(u1y-u2y);
43  T_R const diffz = amrex::Math::abs(u1z-u2z);
44  T_R const diffm = std::sqrt(diffx*diffx+diffy*diffy+diffz*diffz);
45  T_R const summm = std::sqrt(u1x*u1x+u1y*u1y+u1z*u1z) + std::sqrt(u2x*u2x+u2y*u2y+u2z*u2z);
46  // If g = u1 - u2 = 0, do not collide.
47  // Or if the relative difference is less than 1.0e-10.
48  if ( diffm < std::numeric_limits<T_R>::min() || diffm/summm < 1.0e-10 ) {
49  return;
50  }
51 
52  // Transform to SI units (u is the proper velocity here)
53  if (normalized_units) {
54  m1 *= PhysConstSI::m_e;
55  m2 *= PhysConstSI::m_e;
56  u1x *= PhysConstSI::c;
57  u1y *= PhysConstSI::c;
58  u1z *= PhysConstSI::c;
59  u2x *= PhysConstSI::c;
60  u2y *= PhysConstSI::c;
61  u2z *= PhysConstSI::c;
62  }
63  // Compute momenta
64  T_R const p1x = u1x * m1;
65  T_R const p1y = u1y * m1;
66  T_R const p1z = u1z * m1;
67  T_R const p2x = u2x * m2;
68  T_R const p2y = u2y * m2;
69  T_R const p2z = u2z * m2;
70 
71  // Compute center-of-mass (COM) velocity and gamma
72  T_R const mass_g = m1 * g1 + m2 * g2;
73  T_R const vcx = (p1x+p2x) / mass_g;
74  T_R const vcy = (p1y+p2y) / mass_g;
75  T_R const vcz = (p1z+p2z) / mass_g;
76  T_R const vcms = vcx*vcx + vcy*vcy + vcz*vcz;
77  T_R const gc = T_R(1.0) / std::sqrt( T_R(1.0) - vcms*inv_c2_SI );
78 
79  // Compute vc dot v1 and v2
80  T_R const vcDv1 = (vcx*u1x + vcy*u1y + vcz*u1z) / g1;
81  T_R const vcDv2 = (vcx*u2x + vcy*u2y + vcz*u2z) / g2;
82 
83  // Compute p1 star
84  T_R p1sx;
85  T_R p1sy;
86  T_R p1sz;
87  if ( vcms > std::numeric_limits<T_R>::min() )
88  {
89  T_R const lorentz_transform_factor =
90  ( (gc-T_R(1.0))/vcms*vcDv1 - gc )*m1*g1;
91  p1sx = p1x + vcx*lorentz_transform_factor;
92  p1sy = p1y + vcy*lorentz_transform_factor;
93  p1sz = p1z + vcz*lorentz_transform_factor;
94  }
95  else // If vcms = 0, don't do Lorentz-transform.
96  {
97  p1sx = p1x;
98  p1sy = p1y;
99  p1sz = p1z;
100  }
101  const double p1sm = std::sqrt( static_cast<double>(p1sx)*p1sx + p1sy*p1sy + p1sz*p1sz );
102 
103  // Compute gamma star
104  T_R const g1s = ( T_R(1.0) - vcDv1*inv_c2_SI )*gc*g1;
105  T_R const g2s = ( T_R(1.0) - vcDv2*inv_c2_SI )*gc*g2;
106 
107  // Compute the Coulomb log lnLmd
108  double lnLmd;
109  if ( L > T_R(0.0) ) { lnLmd = L; }
110  else
111  {
112  // Compute b0 according to eq (22) from Perez et al., Phys.Plasmas.19.083104 (2012)
113  // Note: there is a typo in the equation, the last square is incorrect!
114  // See the SMILEI documentation: https://smileipic.github.io/Smilei/Understand/collisions.html
115  const double b0 = amrex::Math::abs(q1*q2) * inv_c2_SI /
116  (T_R(4.0)*MathConst::pi*PhysConstSI::ep0) * gc/mass_g *
117  ( m1*g1s*m2*g2s/(p1sm*p1sm*inv_c2_SI) + T_R(1.0) );
118 
119  // Compute the minimal impact parameter
120  double bmin = amrex::max(PhysConstSI::hbar*MathConst::pi/p1sm,b0);
121 
122  // Compute the Coulomb log lnLmd according to eq (23) from Perez et al., Phys.Plasmas.19.083104 (2012)
123  lnLmd = amrex::max( static_cast<double>(2.0),
124  T_R(0.5)*std::log(T_R(1.0)+lmdD*lmdD/(bmin*bmin)) );
125  }
126 
127  // Compute s according to eq (17) from Perez et al., Phys.Plasmas.19.083104 (2012)
128  const double tts = static_cast<double>(m1)*g1s*m2*g2s/(inv_c2_SI*p1sm*p1sm) + T_R(1.0);
129  const double tts2 = tts*tts;
130  const double charge_fac = normalized_units ? static_cast<double>(PhysConstSI::q_e)*PhysConstSI::q_e*
132  T_R s = static_cast<double>(n1)*n2/n12 * dt*lnLmd*q1*q1*q2*q2*charge_fac * inv_c2_SI*inv_c2_SI /
133  ( static_cast<double>(4.0) * MathConst::pi * PhysConstSI::ep0 * PhysConstSI::ep0 *
134  m1*g1*m2*g2 ) * gc*p1sm/mass_g * tts2;
135 
136  const auto cbrt_n1 = std::cbrt(n1);
137  const auto cbrt_n2 = std::cbrt(n2);
138  const auto coeff = static_cast<T_R>(
139  std::pow(4.0*MathConst::pi/3.0,1.0/3.0));
140  T_R const vrel = static_cast<double>(mass_g)*p1sm/(m1*g1s*m2*g2s*gc);
141  T_R const sp = static_cast<double>(coeff) * n1*n2/n12 * dt * vrel * (m1+m2) /
142  amrex::max( m1*cbrt_n1*cbrt_n1,
143  m2*cbrt_n2*cbrt_n2);
144 
145  // Determine s
146  s = amrex::min(s,sp);
147 
148  // Get random numbers
149  T_R r = amrex::Random(engine);
150 
151  // Compute scattering angle according to eq (10) and following from Perez et al., Phys.Plasmas.19.083104 (2012)
152  T_R cosXs;
153  T_R sinXs;
154  if ( s <= T_R(0.1) )
155  {
156  while ( true )
157  {
158  cosXs = T_R(1.0) + s * std::log(r);
159  // Avoid the bug when r is too small such that cosXs < -1
160  if ( cosXs >= T_R(-1.0) ) { break; }
161  r = amrex::Random(engine);
162  }
163  }
164  else if ( s > T_R(0.1) && s <= T_R(3.0) )
165  {
166  T_R const Ainv = static_cast<T_R>(
167  0.0056958 + 0.9560202*s - 0.508139*s*s +
168  0.47913906*s*s*s - 0.12788975*s*s*s*s + 0.02389567*s*s*s*s*s);
169  cosXs = Ainv * std::log( std::exp(T_R(-1.0)/Ainv) +
170  T_R(2.0) * r * std::sinh(T_R(1.0)/Ainv) );
171  }
172  else if ( s > T_R(3.0) && s <= T_R(6.0) )
173  {
174  T_R const A = T_R(3.0) * std::exp(-s);
175  cosXs = T_R(1.0)/A * std::log( std::exp(-A) +
176  T_R(2.0) * r * std::sinh(A) );
177  }
178  else
179  {
180  cosXs = T_R(2.0) * r - T_R(1.0);
181  }
182  sinXs = std::sqrt(T_R(1.0) - cosXs*cosXs);
183 
184  // Get random azimuthal angle
185  T_R const phis = amrex::Random(engine) * T_R(2.0) * MathConst::pi;
186  T_R const cosphis = std::cos(phis);
187  T_R const sinphis = std::sin(phis);
188 
189  // Compute post-collision momenta pfs in COM
190  T_R p1fsx;
191  T_R p1fsy;
192  T_R p1fsz;
193  // p1sp is the p1s perpendicular
194  double p1sp = std::sqrt( static_cast<double>(p1sx)*p1sx + p1sy*p1sy );
195  // Make sure p1sp is not almost zero
196  if ( p1sp > std::numeric_limits<T_R>::min() )
197  {
198  p1fsx = ( p1sx*p1sz/p1sp ) * sinXs*cosphis +
199  ( p1sy*p1sm/p1sp ) * sinXs*sinphis +
200  ( p1sx ) * cosXs;
201  p1fsy = ( p1sy*p1sz/p1sp ) * sinXs*cosphis +
202  (-p1sx*p1sm/p1sp ) * sinXs*sinphis +
203  ( p1sy ) * cosXs;
204  p1fsz = (-p1sp ) * sinXs*cosphis +
205  ( T_R(0.0) ) * sinXs*sinphis +
206  ( p1sz ) * cosXs;
207  // Note a negative sign is different from
208  // Eq. (12) in Perez's paper,
209  // but they are the same due to the random nature of phis.
210  }
211  else
212  {
213  // If the previous p1sp is almost zero
214  // x->y y->z z->x
215  // This set is equivalent to the one in Nanbu's paper
216  p1sp = std::sqrt( static_cast<double>(p1sy)*p1sy + p1sz*p1sz );
217  p1fsy = ( p1sy*p1sx/p1sp ) * sinXs*cosphis +
218  ( p1sz*p1sm/p1sp ) * sinXs*sinphis +
219  ( p1sy ) * cosXs;
220  p1fsz = ( p1sz*p1sx/p1sp ) * sinXs*cosphis +
221  (-p1sy*p1sm/p1sp ) * sinXs*sinphis +
222  ( p1sz ) * cosXs;
223  p1fsx = (-p1sp ) * sinXs*cosphis +
224  ( T_R(0.0) ) * sinXs*sinphis +
225  ( p1sx ) * cosXs;
226 
227  }
228 
229  T_R const p2fsx = -p1fsx;
230  T_R const p2fsy = -p1fsy;
231  T_R const p2fsz = -p1fsz;
232 
233  // Transform from COM to lab frame
234  T_R p1fx; T_R p2fx;
235  T_R p1fy; T_R p2fy;
236  T_R p1fz; T_R p2fz;
237  if ( vcms > std::numeric_limits<T_R>::min() )
238  {
239  T_R const vcDp1fs = vcx*p1fsx + vcy*p1fsy + vcz*p1fsz;
240  T_R const vcDp2fs = vcx*p2fsx + vcy*p2fsy + vcz*p2fsz;
241  T_R const factor = (gc-T_R(1.0))/vcms;
242  T_R const factor1 = factor*vcDp1fs + m1*g1s*gc;
243  T_R const factor2 = factor*vcDp2fs + m2*g2s*gc;
244  p1fx = p1fsx + vcx * factor1;
245  p1fy = p1fsy + vcy * factor1;
246  p1fz = p1fsz + vcz * factor1;
247  p2fx = p2fsx + vcx * factor2;
248  p2fy = p2fsy + vcy * factor2;
249  p2fz = p2fsz + vcz * factor2;
250  }
251  else // If vcms = 0, don't do Lorentz-transform.
252  {
253  p1fx = p1fsx;
254  p1fy = p1fsy;
255  p1fz = p1fsz;
256  p2fx = p2fsx;
257  p2fy = p2fsy;
258  p2fz = p2fsz;
259  }
260 
261  // Rejection method according to eq (14) and following from Perez et al., Phys.Plasmas.19.083104 (2012)
262  r = amrex::Random(engine);
263  if ( w2 > r*amrex::max(w1, w2) )
264  {
265  u1x = p1fx / m1;
266  u1y = p1fy / m1;
267  u1z = p1fz / m1;
268  AMREX_ASSERT_WITH_MESSAGE(!std::isnan(u1x+u1y+u1z+u2x+u2y+u2z),
269  "NaNs detected in the updated momentum of the collision module. "
270  "Most likely due to a negative root in the calculation of gc");
271  AMREX_ASSERT(!std::isinf(u1x+u1y+u1z+u2x+u2y+u2z));
272  }
273  r = amrex::Random(engine);
274  if ( w1 > r*amrex::max(w1, w2) )
275  {
276  u2x = p2fx / m2;
277  u2y = p2fy / m2;
278  u2z = p2fz / m2;
279  AMREX_ASSERT_WITH_MESSAGE(!std::isnan(u1x+u1y+u1z+u2x+u2y+u2z),
280  "NaNs detected in the updated momentum of the collision module. "
281  "Most likely due to a negative root in the calculation of gc");
282  AMREX_ASSERT(!std::isinf(u1x+u1y+u1z+u2x+u2y+u2z));
283  }
284  if (normalized_units) { // backtransform to normalized units
285  u1x *= inv_c_SI;
286  u1y *= inv_c_SI;
287  u1z *= inv_c_SI;
288  u2x *= inv_c_SI;
289  u2y *= inv_c_SI;
290  u2z *= inv_c_SI;
291  }
292 }
293 
294 #endif // HIPACE_UPDATE_MOMENTUM_PEREZ_ELASTIC_H_
#define AMREX_ASSERT_WITH_MESSAGE(EX, MSG)
#define AMREX_ASSERT(EX)
#define AMREX_INLINE
#define AMREX_GPU_HOST_DEVICE
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real pow(amrex::Real base)
calculate low integer powers base^exp
Definition: OpenBoundary.H:18
AMREX_GPU_HOST_DEVICE AMREX_INLINE void UpdateMomentumPerezElastic(T_R &u1x, T_R &u1y, T_R &u1z, T_R const g1, T_R &u2x, T_R &u2y, T_R &u2z, T_R const g2, T_R const n1, T_R const n2, T_R const n12, T_R const q1, T_R m1, T_R const w1, T_R const q2, T_R m2, T_R const w2, T_R const dt, T_R const L, T_R const lmdD, T_R const inv_c_SI, T_R const inv_c2_SI, const bool normalized_units, amrex::RandomEngine const &engine)
Definition: UpdateMomentumPerez.H:30
static constexpr amrex::Real pi
Definition: Constants.H:30
static constexpr auto q_e
Definition: Constants.H:20
static constexpr auto hbar
Definition: Constants.H:23
static constexpr auto m_e
Definition: Constants.H:21
static constexpr auto c
Definition: Constants.H:17
static constexpr auto ep0
Definition: Constants.H:18
Real Random()
AMREX_GPU_HOST_DEVICE constexpr AMREX_FORCE_INLINE const T & max(const T &a, const T &b) noexcept
AMREX_GPU_HOST_DEVICE constexpr AMREX_FORCE_INLINE const T & min(const T &a, const T &b) noexcept