WarpX
UpdateMomentumPerezElastic.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 #ifndef WARPX_PARTICLES_COLLISION_UPDATE_MOMENTUM_PEREZ_ELASTIC_H_
8 #define WARPX_PARTICLES_COLLISION_UPDATE_MOMENTUM_PEREZ_ELASTIC_H_
9 
10 #include "Utils/WarpXConst.H"
11 
12 #include <AMReX_Math.H>
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  * Updates and corrections to the original publication are documented in
28  * https://github.com/ECP-WarpX/WarpX/issues/429
29  * https://github.com/ECP-WarpX/WarpX/files/3799803/main.pdf
30  */
31 
32 template <typename T_PR, typename T_R>
35  T_PR& u1x, T_PR& u1y, T_PR& u1z, T_PR& u2x, T_PR& u2y, T_PR& u2z,
36  T_PR const n1, T_PR const n2, T_PR const n12,
37  T_PR const q1, T_PR const m1, T_PR const w1,
38  T_PR const q2, T_PR const m2, T_PR const w2,
39  T_R const dt, T_PR const L, T_PR const lmdD,
40  amrex::RandomEngine const& engine)
41 {
42 
43  T_PR const diffx = amrex::Math::abs(u1x-u2x);
44  T_PR const diffy = amrex::Math::abs(u1y-u2y);
45  T_PR const diffz = amrex::Math::abs(u1z-u2z);
46  T_PR const diffm = std::sqrt(diffx*diffx+diffy*diffy+diffz*diffz);
47  T_PR const summm = std::sqrt(u1x*u1x+u1y*u1y+u1z*u1z) + std::sqrt(u2x*u2x+u2y*u2y+u2z*u2z);
48  // If g = u1 - u2 = 0, do not collide.
49  // Or if the relative difference is less than 1.0e-10.
50  if ( diffm < std::numeric_limits<T_PR>::min() || diffm/summm < 1.0e-10 ) { return; }
51 
52  T_PR constexpr inv_c2 = T_PR(1.0) / ( PhysConst::c * PhysConst::c );
53 
54  // Compute Lorentz factor gamma
55  T_PR const g1 = std::sqrt( T_PR(1.0) + (u1x*u1x+u1y*u1y+u1z*u1z)*inv_c2 );
56  T_PR const g2 = std::sqrt( T_PR(1.0) + (u2x*u2x+u2y*u2y+u2z*u2z)*inv_c2 );
57 
58  // Compute momenta
59  T_PR const p1x = u1x * m1;
60  T_PR const p1y = u1y * m1;
61  T_PR const p1z = u1z * m1;
62  T_PR const p2x = u2x * m2;
63  T_PR const p2y = u2y * m2;
64  T_PR const p2z = u2z * m2;
65 
66  // Compute center-of-mass (COM) velocity and gamma
67  T_PR const mass_g = m1 * g1 + m2 * g2;
68  T_PR const vcx = (p1x+p2x) / mass_g;
69  T_PR const vcy = (p1y+p2y) / mass_g;
70  T_PR const vcz = (p1z+p2z) / mass_g;
71  T_PR const vcms = vcx*vcx + vcy*vcy + vcz*vcz;
72  T_PR const gc = T_PR(1.0) / std::sqrt( T_PR(1.0) - vcms*inv_c2 );
73 
74  // Compute vc dot v1 and v2
75  T_PR const vcDv1 = (vcx*u1x + vcy*u1y + vcz*u1z) / g1;
76  T_PR const vcDv2 = (vcx*u2x + vcy*u2y + vcz*u2z) / g2;
77 
78  // Compute p1 star
79  T_PR p1sx;
80  T_PR p1sy;
81  T_PR p1sz;
82  if ( vcms > std::numeric_limits<T_PR>::min() )
83  {
84  /* lorentz_transform_factor = ( (gc-1.0)/vcms*vcDv1 - gc )*m1*g1;
85  * Rewrite to avoid loss of precision from subtracting similar
86  * numbers when gc is close to 1 */
87  T_PR const lorentz_transform_factor =
88  ( (gc*gc*vcms*inv_c2/(T_PR(1.0) + gc))/vcms*vcDv1 - gc )*m1*g1;
89  p1sx = p1x + vcx*lorentz_transform_factor;
90  p1sy = p1y + vcy*lorentz_transform_factor;
91  p1sz = p1z + vcz*lorentz_transform_factor;
92  }
93  else // If vcms = 0, don't do Lorentz-transform.
94  {
95  p1sx = p1x;
96  p1sy = p1y;
97  p1sz = p1z;
98  }
99  T_PR const p1sm = std::sqrt( p1sx*p1sx + p1sy*p1sy + p1sz*p1sz );
100 
101  // Compute gamma star
102  T_PR const g1s = ( T_PR(1.0) - vcDv1*inv_c2 )*gc*g1;
103  T_PR const g2s = ( T_PR(1.0) - vcDv2*inv_c2 )*gc*g2;
104 
105  // Compute s
106  T_PR s = 0;
107  if (p1sm > std::numeric_limits<T_PR>::min()) {
108 
109  // s is non-zero (i.e. particles scatter) only if the relative
110  // motion between particles is not negligible (p1sm non-zero)
111 
112  // Compute the Coulomb log lnLmd first
113  T_PR lnLmd;
114  if ( L > T_PR(0.0) ) { lnLmd = L; }
115  else
116  {
117  // Compute b0 according to eq (22) from Perez et al., Phys.Plasmas.19.083104 (2012)
118  // Note: there is a typo in the equation, the last square is incorrect!
119  // See the SMILEI documentation: https://smileipic.github.io/Smilei/Understand/collisions.html
120  // and https://github.com/ECP-WarpX/WarpX/files/3799803/main.pdf from GitHub #429
121  T_PR const b0 = amrex::Math::abs(q1*q2) * inv_c2 /
122  (T_PR(4.0)*MathConst::pi*PhysConst::ep0) * gc/mass_g *
123  ( m1*g1s*m2*g2s/(p1sm*p1sm*inv_c2) + T_PR(1.0) );
124 
125  // Compute the minimal impact parameter
126  constexpr T_PR hbar_pi = static_cast<T_PR>(PhysConst::hbar*MathConst::pi);
127  const T_PR bmin = amrex::max(hbar_pi/p1sm, b0);
128 
129  // Compute the Coulomb log lnLmd
130  lnLmd = amrex::max( T_PR(2.0),
131  T_PR(0.5)*std::log(T_PR(1.0)+lmdD*lmdD/(bmin*bmin)) );
132  }
133 
134  // Compute s
135  const auto tts = m1*g1s*m2*g2s/(inv_c2*p1sm*p1sm) + T_PR(1.0);
136  const auto tts2 = tts*tts;
137  s = n1*n2/n12 * dt*lnLmd*q1*q1*q2*q2 /
138  ( T_PR(4.0) * MathConst::pi * PhysConst::ep0 * PhysConst::ep0 *
139  m1*g1*m2*g2/(inv_c2*inv_c2) ) * gc*p1sm/mass_g * tts2;
140 
141  // Compute s'
142  const auto cbrt_n1 = std::cbrt(n1);
143  const auto cbrt_n2 = std::cbrt(n2);
144  const auto coeff = static_cast<T_PR>(
145  std::pow(4.0*MathConst::pi/3.0,1.0/3.0));
146  T_PR const vrel = mass_g*p1sm/(m1*g1s*m2*g2s*gc);
147  T_PR const sp = coeff * n1*n2/n12 * dt * vrel * (m1+m2) /
148  amrex::max( m1*cbrt_n1*cbrt_n1,
149  m2*cbrt_n2*cbrt_n2);
150 
151  // Determine s
152  s = amrex::min(s,sp);
153  }
154 
155  // Only modify momenta if is s is non-zero
156  if (s > std::numeric_limits<T_PR>::min()) {
157 
158  // Get random numbers
159  T_PR r = amrex::Random(engine);
160 
161  // Compute scattering angle
162  T_PR cosXs;
163  T_PR sinXs;
164  if ( s <= T_PR(0.1) )
165  {
166  while ( true )
167  {
168  cosXs = T_PR(1.0) + s * std::log(r);
169  // Avoid the bug when r is too small such that cosXs < -1
170  if ( cosXs >= T_PR(-1.0) ) { break; }
171  r = amrex::Random(engine);
172  }
173  }
174  else if ( s > T_PR(0.1) && s <= T_PR(3.0) )
175  {
176  T_PR const Ainv = static_cast<T_PR>(
177  0.0056958 + 0.9560202*s - 0.508139*s*s +
178  0.47913906*s*s*s - 0.12788975*s*s*s*s + 0.02389567*s*s*s*s*s);
179  cosXs = Ainv * std::log( std::exp(T_PR(-1.0)/Ainv) +
180  T_PR(2.0) * r * std::sinh(T_PR(1.0)/Ainv) );
181  }
182  else if ( s > T_PR(3.0) && s <= T_PR(6.0) )
183  {
184  T_PR const A = T_PR(3.0) * std::exp(-s);
185  cosXs = T_PR(1.0)/A * std::log( std::exp(-A) +
186  T_PR(2.0) * r * std::sinh(A) );
187  }
188  else
189  {
190  cosXs = T_PR(2.0) * r - T_PR(1.0);
191  }
192  sinXs = std::sqrt(T_PR(1.0) - cosXs*cosXs);
193 
194  // Get random azimuthal angle
195  T_PR const phis = amrex::Random(engine) * T_PR(2.0) * MathConst::pi;
196  T_PR const cosphis = std::cos(phis);
197  T_PR const sinphis = std::sin(phis);
198 
199  // Compute post-collision momenta pfs in COM
200  T_PR p1fsx;
201  T_PR p1fsy;
202  T_PR p1fsz;
203  // p1sp is the p1s perpendicular
204  T_PR p1sp = std::sqrt( p1sx*p1sx + p1sy*p1sy );
205  // Make sure p1sp is not almost zero
206  if ( p1sp > std::numeric_limits<T_PR>::min() )
207  {
208  p1fsx = ( p1sx*p1sz/p1sp ) * sinXs*cosphis +
209  ( p1sy*p1sm/p1sp ) * sinXs*sinphis +
210  ( p1sx ) * cosXs;
211  p1fsy = ( p1sy*p1sz/p1sp ) * sinXs*cosphis +
212  (-p1sx*p1sm/p1sp ) * sinXs*sinphis +
213  ( p1sy ) * cosXs;
214  p1fsz = (-p1sp ) * sinXs*cosphis +
215  ( T_PR(0.0) ) * sinXs*sinphis +
216  ( p1sz ) * cosXs;
217  // Note a negative sign is different from
218  // Eq. (12) in Perez's paper,
219  // but they are the same due to the random nature of phis.
220  }
221  else
222  {
223  // If the previous p1sp is almost zero
224  // x->y y->z z->x
225  // This set is equivalent to the one in Nanbu's paper
226  p1sp = std::sqrt( p1sy*p1sy + p1sz*p1sz );
227  p1fsy = ( p1sy*p1sx/p1sp ) * sinXs*cosphis +
228  ( p1sz*p1sm/p1sp ) * sinXs*sinphis +
229  ( p1sy ) * cosXs;
230  p1fsz = ( p1sz*p1sx/p1sp ) * sinXs*cosphis +
231  (-p1sy*p1sm/p1sp ) * sinXs*sinphis +
232  ( p1sz ) * cosXs;
233  p1fsx = (-p1sp ) * sinXs*cosphis +
234  ( T_PR(0.0) ) * sinXs*sinphis +
235  ( p1sx ) * cosXs;
236  }
237 
238  T_PR const p2fsx = -p1fsx;
239  T_PR const p2fsy = -p1fsy;
240  T_PR const p2fsz = -p1fsz;
241 
242  // Transform from COM to lab frame
243  T_PR p1fx; T_PR p2fx;
244  T_PR p1fy; T_PR p2fy;
245  T_PR p1fz; T_PR p2fz;
246  if ( vcms > std::numeric_limits<T_PR>::min() )
247  {
248  T_PR const vcDp1fs = vcx*p1fsx + vcy*p1fsy + vcz*p1fsz;
249  T_PR const vcDp2fs = vcx*p2fsx + vcy*p2fsy + vcz*p2fsz;
250  /* factor = (gc-1.0)/vcms; Rewrite to avoid subtraction losing precision when gc is close to 1 */
251  T_PR const factor = gc*gc*inv_c2/(gc+T_PR(1.0));
252  T_PR const factor1 = factor*vcDp1fs + m1*g1s*gc;
253  T_PR const factor2 = factor*vcDp2fs + m2*g2s*gc;
254  p1fx = p1fsx + vcx * factor1;
255  p1fy = p1fsy + vcy * factor1;
256  p1fz = p1fsz + vcz * factor1;
257  p2fx = p2fsx + vcx * factor2;
258  p2fy = p2fsy + vcy * factor2;
259  p2fz = p2fsz + vcz * factor2;
260  }
261  else // If vcms = 0, don't do Lorentz-transform.
262  {
263  p1fx = p1fsx;
264  p1fy = p1fsy;
265  p1fz = p1fsz;
266  p2fx = p2fsx;
267  p2fy = p2fsy;
268  p2fz = p2fsz;
269  }
270 
271  // Rejection method
272  r = amrex::Random(engine);
273  if ( w2 > r*amrex::max(w1, w2) )
274  {
275  u1x = p1fx / m1;
276  u1y = p1fy / m1;
277  u1z = p1fz / m1;
278  }
279  r = amrex::Random(engine);
280  if ( w1 > r*amrex::max(w1, w2) )
281  {
282  u2x = p2fx / m2;
283  u2y = p2fy / m2;
284  u2z = p2fz / m2;
285  }
286 #ifndef AMREX_USE_DPCPP
287  AMREX_ASSERT(!std::isnan(u1x+u1y+u1z+u2x+u2y+u2z));
288  AMREX_ASSERT(!std::isinf(u1x+u1y+u1z+u2x+u2y+u2z));
289 #endif
290 
291  } // if s > std::numeric_limits<T_PR>::min()
292 
293 }
294 
295 #endif // WARPX_PARTICLES_COLLISION_UPDATE_MOMENTUM_PEREZ_ELASTIC_H_
#define AMREX_ASSERT(EX)
#define AMREX_INLINE
#define AMREX_GPU_HOST_DEVICE
AMREX_GPU_HOST_DEVICE AMREX_INLINE void UpdateMomentumPerezElastic(T_PR &u1x, T_PR &u1y, T_PR &u1z, T_PR &u2x, T_PR &u2y, T_PR &u2z, T_PR const n1, T_PR const n2, T_PR const n12, T_PR const q1, T_PR const m1, T_PR const w1, T_PR const q2, T_PR const m2, T_PR const w2, T_R const dt, T_PR const L, T_PR const lmdD, amrex::RandomEngine const &engine)
Definition: UpdateMomentumPerezElastic.H:34
static constexpr auto c
vacuum speed of light [m/s]
Definition: constant.H:44
static constexpr auto ep0
vacuum permittivity: dielectric permittivity of vacuum [F/m]
Definition: constant.H:46
static constexpr auto hbar
reduced Planck Constant = h / tau [J*s]
Definition: constant.H:59
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
float dt
Definition: stencil.py:442
sp
Definition: stencil.py:27