7 #ifndef WARPX_PARTICLES_COLLISION_UPDATE_MOMENTUM_PEREZ_ELASTIC_H_
8 #define WARPX_PARTICLES_COLLISION_UPDATE_MOMENTUM_PEREZ_ELASTIC_H_
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,
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);
50 if ( diffm < std::numeric_limits<T_PR>::min() || diffm/summm < 1.0e-10 ) {
return; }
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 );
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;
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 );
75 T_PR
const vcDv1 = (vcx*u1x + vcy*u1y + vcz*u1z) / g1;
76 T_PR
const vcDv2 = (vcx*u2x + vcy*u2y + vcz*u2z) / g2;
82 if ( vcms > std::numeric_limits<T_PR>::min() )
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;
99 T_PR
const p1sm = std::sqrt( p1sx*p1sx + p1sy*p1sy + p1sz*p1sz );
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;
107 if (p1sm > std::numeric_limits<T_PR>::min()) {
114 if ( L > T_PR(0.0) ) { lnLmd = L; }
121 T_PR
const b0 = amrex::Math::abs(q1*q2) * inv_c2 /
123 ( m1*g1s*m2*g2s/(p1sm*p1sm*inv_c2) + T_PR(1.0) );
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);
131 T_PR(0.5)*std::log(T_PR(1.0)+lmdD*lmdD/(bmin*bmin)) );
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 /
139 m1*g1*m2*g2/(inv_c2*inv_c2) ) * gc*p1sm/mass_g * tts2;
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) /
156 if (s > std::numeric_limits<T_PR>::min()) {
164 if ( s <= T_PR(0.1) )
168 cosXs = T_PR(1.0) + s * std::log(
r);
170 if ( cosXs >= T_PR(-1.0) ) {
break; }
174 else if ( s > T_PR(0.1) && s <= T_PR(3.0) )
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) );
182 else if ( s > T_PR(3.0) && s <= T_PR(6.0) )
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) );
190 cosXs = T_PR(2.0) *
r - T_PR(1.0);
192 sinXs = std::sqrt(T_PR(1.0) - cosXs*cosXs);
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);
204 T_PR p1sp = std::sqrt( p1sx*p1sx + p1sy*p1sy );
206 if ( p1sp > std::numeric_limits<T_PR>::min() )
208 p1fsx = ( p1sx*p1sz/p1sp ) * sinXs*cosphis +
209 ( p1sy*p1sm/p1sp ) * sinXs*sinphis +
211 p1fsy = ( p1sy*p1sz/p1sp ) * sinXs*cosphis +
212 (-p1sx*p1sm/p1sp ) * sinXs*sinphis +
214 p1fsz = (-p1sp ) * sinXs*cosphis +
215 ( T_PR(0.0) ) * sinXs*sinphis +
226 p1sp = std::sqrt( p1sy*p1sy + p1sz*p1sz );
227 p1fsy = ( p1sy*p1sx/p1sp ) * sinXs*cosphis +
228 ( p1sz*p1sm/p1sp ) * sinXs*sinphis +
230 p1fsz = ( p1sz*p1sx/p1sp ) * sinXs*cosphis +
231 (-p1sy*p1sm/p1sp ) * sinXs*sinphis +
233 p1fsx = (-p1sp ) * sinXs*cosphis +
234 ( T_PR(0.0) ) * sinXs*sinphis +
238 T_PR
const p2fsx = -p1fsx;
239 T_PR
const p2fsy = -p1fsy;
240 T_PR
const p2fsz = -p1fsz;
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() )
248 T_PR
const vcDp1fs = vcx*p1fsx + vcy*p1fsy + vcz*p1fsz;
249 T_PR
const vcDp2fs = vcx*p2fsx + vcy*p2fsy + vcz*p2fsz;
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;
286 #ifndef AMREX_USE_DPCPP
#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
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