8 #ifndef HIPACE_UPDATE_MOMENTUM_PEREZ_ELASTIC_H_
9 #define HIPACE_UPDATE_MOMENTUM_PEREZ_ELASTIC_H_
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,
39 using namespace amrex::literals;
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);
48 if ( diffm < std::numeric_limits<T_R>::min() || diffm/summm < 1.0e-10 ) {
53 if (normalized_units) {
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;
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 );
80 T_R
const vcDv1 = (vcx*u1x + vcy*u1y + vcz*u1z) / g1;
81 T_R
const vcDv2 = (vcx*u2x + vcy*u2y + vcz*u2z) / g2;
87 if ( vcms > std::numeric_limits<T_R>::min() )
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;
101 const double p1sm = std::sqrt(
static_cast<double>(p1sx)*p1sx + p1sy*p1sy + p1sz*p1sz );
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;
109 if ( L > T_R(0.0) ) { lnLmd = L; }
115 const double b0 = amrex::Math::abs(q1*q2) * inv_c2_SI /
117 ( m1*g1s*m2*g2s/(p1sm*p1sm*inv_c2_SI) + T_R(1.0) );
124 T_R(0.5)*std::log(T_R(1.0)+lmdD*lmdD/(bmin*bmin)) );
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;
132 T_R s =
static_cast<double>(n1)*n2/n12 * dt*lnLmd*q1*q1*q2*q2*charge_fac * inv_c2_SI*inv_c2_SI /
134 m1*g1*m2*g2 ) * gc*p1sm/mass_g * tts2;
136 const auto cbrt_n1 = std::cbrt(n1);
137 const auto cbrt_n2 = std::cbrt(n2);
138 const auto coeff =
static_cast<T_R
>(
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) /
158 cosXs = T_R(1.0) + s * std::log(
r);
160 if ( cosXs >= T_R(-1.0) ) {
break; }
164 else if ( s > T_R(0.1) && s <= T_R(3.0) )
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) );
172 else if ( s > T_R(3.0) && s <= T_R(6.0) )
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) );
180 cosXs = T_R(2.0) *
r - T_R(1.0);
182 sinXs = std::sqrt(T_R(1.0) - cosXs*cosXs);
186 T_R
const cosphis = std::cos(phis);
187 T_R
const sinphis = std::sin(phis);
194 double p1sp = std::sqrt(
static_cast<double>(p1sx)*p1sx + p1sy*p1sy );
196 if ( p1sp > std::numeric_limits<T_R>::min() )
198 p1fsx = ( p1sx*p1sz/p1sp ) * sinXs*cosphis +
199 ( p1sy*p1sm/p1sp ) * sinXs*sinphis +
201 p1fsy = ( p1sy*p1sz/p1sp ) * sinXs*cosphis +
202 (-p1sx*p1sm/p1sp ) * sinXs*sinphis +
204 p1fsz = (-p1sp ) * sinXs*cosphis +
205 ( T_R(0.0) ) * sinXs*sinphis +
216 p1sp = std::sqrt(
static_cast<double>(p1sy)*p1sy + p1sz*p1sz );
217 p1fsy = ( p1sy*p1sx/p1sp ) * sinXs*cosphis +
218 ( p1sz*p1sm/p1sp ) * sinXs*sinphis +
220 p1fsz = ( p1sz*p1sx/p1sp ) * sinXs*cosphis +
221 (-p1sy*p1sm/p1sp ) * sinXs*sinphis +
223 p1fsx = (-p1sp ) * sinXs*cosphis +
224 ( T_R(0.0) ) * sinXs*sinphis +
229 T_R
const p2fsx = -p1fsx;
230 T_R
const p2fsy = -p1fsy;
231 T_R
const p2fsz = -p1fsz;
237 if ( vcms > std::numeric_limits<T_R>::min() )
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;
269 "NaNs detected in the updated momentum of the collision module. "
270 "Most likely due to a negative root in the calculation of gc");
280 "NaNs detected in the updated momentum of the collision module. "
281 "Most likely due to a negative root in the calculation of gc");
284 if (normalized_units) {
#define AMREX_ASSERT_WITH_MESSAGE(EX, MSG)
#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
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