Hipace
ElasticCollisionPerez.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 HIPACE_ELASTIC_COLLISION_PEREZ_H_
8 #define HIPACE_ELASTIC_COLLISION_PEREZ_H_
9 
10 #include "utils/Constants.H"
11 #include "UpdateMomentumPerez.H"
12 #include "ComputeTemperature.H"
13 
14 #include <AMReX_Random.H>
15 
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,
67  T_R q1, T_R q2,
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,
74  amrex::RandomEngine const& engine)
75 {
76  using namespace amrex::literals;
77 
78  const int NI1 = I1e - I1s;
79  const int NI2 = I2e - I2s;
80 
81  // get local T1t and T2t
82  T_R T1t; T_R T2t;
83  if ( T1 <= T_R(0.0) && L <= T_R(0.0) )
84  {
85  T1t = ComputeTemperature(I1s,I1e,I1,u1x,u1y,psi1,m1,clight,inv_c2,is_beam_coll);
86  }
87  else { T1t = T1; }
88  if ( T2 <= T_R(0.0) && L <= T_R(0.0) )
89  {
90  // second species is always plasma thus is_beam_coll is always false
91  T2t = ComputeTemperature(I2s,I2e,I2,u2x,u2y,psi2,m2,clight,inv_c2,false);
92  }
93  else { T2t = T2; }
94 
95  // local density
96  T_R n1 = T_R(0.0);
97  T_R n2 = T_R(0.0);
98  T_R n12 = 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) {
102  n1 = n1 + n2;
103  n2 = n1;
104  }
105  if (n1 == 0 || n2 == 0) return;
106  // compute n12 according to eq. 16 in Perez et al., Phys. Plasmas 19 (8) (2012) 083104
107  {
108  int i1 = I1s; int i2 = I2s;
109  for (int k = 0; k < amrex::max(NI1,NI2); ++k)
110  {
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; }
114  }
115  if (is_same_species) n12*= T_R(2.0);
116  }
117 
118  // in normalized units, the weights already represent a density, so it must only be rescaled
119  // to SI units. In SI units, the weights represent the amount of particles, so it must be
120  // divided by the volume
121  if (normalized_units) {
122  n1 *= background_density_SI;
123  n2 *= background_density_SI;
124  n12 *= background_density_SI;
125  } else {
126  n1 *= inv_dV;
127  n2 *= inv_dV;
128  n12 *= inv_dV;
129  }
130 
131  // compute Debye length lmdD
132  T_R lmdD;
133  if ( T1t <= T_R(0.0) || T2t <= T_R(0.0) ) {
134  lmdD = T_R(0.0);
135  }
136  else {
137  lmdD = T_R(1.0)/std::sqrt( n1*q1*q1/(T1t*PhysConstSI::ep0) +
138  n2*q2*q2/(T2t*PhysConstSI::ep0) );
139  }
140  // minimum mean interatomic distance rmin (see Perez et al., Phys. Plasmas 19 (8) (2012) 083104)
141  T_R rmin = std::pow( T_R(4.0) * MathConst::pi / T_R(3.0) *
142  amrex::max(n1,n2), T_R(-1.0/3.0) );
143  lmdD = amrex::max(lmdD, rmin);
144 
145  // call UpdateMomentumPerezElastic()
146  {
147  int i1 = I1s; int i2 = I2s;
148  for (int k = 0; k < amrex::max(NI1,NI2); ++k)
149  {
150  // adjust the charge to the ionization level of the ion. This assumes that the impact
151  // parameter is much larger than the atomic radius, as the bound electrons are not
152  // treated separately
153  if (can_ionize1) q1 *= ion_lev1[I1[i1]];
154  if (can_ionize2) q2 *= ion_lev2[I2[i2]];
155  // particle's Lorentz factor
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]] );
160  // particle's Lorentz factor
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]] );
163 
164  // Convert from pseudo-potential to momentum
165  amrex::Real u1z = is_beam_coll ? psi1[I1[i1]] : clight * (g1 - psi1[I1[i1]]);
166  amrex::Real u2z = clight * (g2 - psi2[I2[i2]]);
167 
168  // In the longitudinal push of plasma particles, the dt is different for each particle.
169  // The dt applied for collision probability is the average (in the lab frame) of these
170  // dts. This is NOT clean. TODO FIXME.
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);
177 
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;
184 
185  ++i1; if ( i1 == static_cast<int>(I1e) ) { i1 = I1s; }
186  ++i2; if ( i2 == static_cast<int>(I2e) ) { i2 = I2s; }
187  }
188  }
189 }
190 
191 #endif // HIPACE_ELASTIC_COLLISION_PEREZ_H_
#define AMREX_INLINE
#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