7 #ifndef HIPACE_ELASTIC_COLLISION_PEREZ_H_
8 #define HIPACE_ELASTIC_COLLISION_PEREZ_H_
57 template <
typename T_index,
typename T_R>
60 T_index
const I1s, T_index
const I1e,
61 T_index
const I2s, T_index
const I2e,
62 T_index *I1, T_index *I2,
63 T_R *u1x, T_R *u1y, T_R *psi1,
64 T_R *u2x, T_R *u2y, T_R *psi2,
65 T_R
const *w1, T_R
const *w2,
66 int const *ion_lev1,
int const *ion_lev2,
68 T_R
const m1, T_R
const m2,
69 T_R
const T1, T_R
const T2,
70 const bool can_ionize1,
const bool can_ionize2,
71 T_R
const dt, T_R
const L, T_R
const inv_dV, T_R
const clight, T_R
const inv_c,
72 T_R
const inv_c_SI, T_R
const inv_c2, T_R
const inv_c2_SI,
const bool normalized_units,
73 const amrex::Real background_density_SI,
const bool is_same_species,
bool is_beam_coll,
76 using namespace amrex::literals;
78 const int NI1 = I1e - I1s;
79 const int NI2 = I2e - I2s;
83 if ( T1 <= T_R(0.0) && L <= T_R(0.0) )
88 if ( T2 <= T_R(0.0) && L <= T_R(0.0) )
99 for (
int i1=I1s; i1<static_cast<int>(I1e); ++i1) { n1 += w1[ I1[i1] ]; }
100 for (
int i2=I2s; i2<static_cast<int>(I2e); ++i2) { n2 += w2[ I2[i2] ]; }
101 if (is_same_species) {
105 if (n1 == 0 || n2 == 0)
return;
108 int i1 = I1s;
int i2 = I2s;
111 n12 +=
amrex::min( w1[ I1[i1] ], w2[ I2[i2] ] );
112 ++i1;
if ( i1 ==
static_cast<int>(I1e) ) { i1 = I1s; }
113 ++i2;
if ( i2 ==
static_cast<int>(I2e) ) { i2 = I2s; }
115 if (is_same_species) n12*= T_R(2.0);
121 if (normalized_units) {
122 n1 *= background_density_SI;
123 n2 *= background_density_SI;
124 n12 *= background_density_SI;
133 if ( T1t <= T_R(0.0) || T2t <= T_R(0.0) ) {
147 int i1 = I1s;
int i2 = I2s;
153 if (can_ionize1) q1 *= ion_lev1[I1[i1]];
154 if (can_ionize2) q2 *= ion_lev2[I2[i2]];
156 amrex::Real g1 = is_beam_coll ? std::sqrt( 1._rt
157 + (u1x[I1[i1]]*u1x[I1[i1]] + u1y[I1[i1]]*u1y[I1[i1]] + psi1[I1[i1]]*psi1[I1[i1]])*inv_c2 )
158 : (1.0_rt + u1x[I1[i1]]*u1x[I1[i1]]*inv_c2 + u1y[I1[i1]]*u1y[I1[i1]]*inv_c2 +
159 psi1[I1[i1]]*psi1[I1[i1]]) / (2.0_rt * psi1[I1[i1]] );
161 amrex::Real g2 = (1.0_rt + u2x[I2[i2]]*u2x[I2[i2]]*inv_c2 + u2y[I2[i2]]*u2y[I2[i2]]*inv_c2 +
162 psi2[I2[i2]]*psi2[I2[i2]]) / (2.0_rt * psi2[I2[i2]] );
165 amrex::Real u1z = is_beam_coll ? psi1[I1[i1]] : clight * (g1 - psi1[I1[i1]]);
166 amrex::Real u2z = clight * (g2 - psi2[I2[i2]]);
171 const amrex::Real dt_fac = is_beam_coll ? 1.0_rt : 0.5_rt * (g1/psi1[I1[i1]] + g2/psi2[I2[i2]]);
173 u1x[ I1[i1] ], u1y[ I1[i1] ], u1z, g1,
174 u2x[ I2[i2] ], u2y[ I2[i2] ], u2z, g2,
175 n1, n2, n12, q1, m1, w1[ I1[i1] ], q2, m2, w2[ I2[i2] ],
176 dt * dt_fac, L, lmdD, inv_c_SI, inv_c2_SI, normalized_units, engine);
178 g1 = std::sqrt( T_R(1.0) +
179 (u1x[I1[i1]]*u1x[I1[i1]]+u1y[I1[i1]]*u1y[I1[i1]]+u1z*u1z)*inv_c2 );
180 psi1[I1[i1]] = is_beam_coll ? u1z : g1 - u1z*inv_c;
181 g2 = std::sqrt( T_R(1.0) +
182 (u2x[I2[i2]]*u2x[I2[i2]]+u2y[I2[i2]]*u2y[I2[i2]]+u2z*u2z)*inv_c2 );
183 psi2[I2[i2]] = g2 - u2z*inv_c;
185 ++i1;
if ( i1 ==
static_cast<int>(I1e) ) { i1 = I1s; }
186 ++i2;
if ( i2 ==
static_cast<int>(I2e) ) { i2 = I2s; }
#define AMREX_GPU_HOST_DEVICE
AMREX_GPU_HOST_DEVICE T_R ComputeTemperature(T_index const Is, T_index const Ie, T_index const *AMREX_RESTRICT I, T_R const *AMREX_RESTRICT ux, T_R const *AMREX_RESTRICT uy, T_R const *AMREX_RESTRICT psi, T_R const m, T_R const clight, T_R const inv_c2, bool is_beam_coll)
Definition: ComputeTemperature.H:13
AMREX_GPU_HOST_DEVICE AMREX_INLINE void ElasticCollisionPerez(T_index const I1s, T_index const I1e, T_index const I2s, T_index const I2e, T_index *I1, T_index *I2, T_R *u1x, T_R *u1y, T_R *psi1, T_R *u2x, T_R *u2y, T_R *psi2, T_R const *w1, T_R const *w2, int const *ion_lev1, int const *ion_lev2, T_R q1, T_R q2, T_R const m1, T_R const m2, T_R const T1, T_R const T2, const bool can_ionize1, const bool can_ionize2, T_R const dt, T_R const L, T_R const inv_dV, T_R const clight, T_R const inv_c, T_R const inv_c_SI, T_R const inv_c2, T_R const inv_c2_SI, const bool normalized_units, const amrex::Real background_density_SI, const bool is_same_species, bool is_beam_coll, amrex::RandomEngine const &engine)
Prepare information for and call UpdateMomentumPerezElastic().
Definition: ElasticCollisionPerez.H:59
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 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