WarpX
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 WARPX_PARTICLES_COLLISION_ELASTIC_COLLISION_PEREZ_H_
8 #define WARPX_PARTICLES_COLLISION_ELASTIC_COLLISION_PEREZ_H_
9 
10 #include "ComputeTemperature.H"
13 #include "Utils/WarpXConst.H"
14 
15 #include <AMReX_Random.H>
16 
17 
43 template <typename T_index, typename T_PR, typename T_R, typename SoaData_type>
46  T_index const I1s, T_index const I1e,
47  T_index const I2s, T_index const I2e,
48  T_index const* AMREX_RESTRICT I1,
49  T_index const* AMREX_RESTRICT I2,
50  SoaData_type soa_1, SoaData_type soa_2,
51  T_PR const q1, T_PR const q2,
52  T_PR const m1, T_PR const m2,
53  T_PR const T1, T_PR const T2,
54  T_R const dt, T_PR const L, T_R const dV,
55  amrex::RandomEngine const& engine,
56  bool const isSameSpecies, T_index coll_idx)
57 {
58  const T_index NI1 = I1e - I1s;
59  const T_index NI2 = I2e - I2s;
60  const T_index max_N = amrex::max(NI1,NI2);
61  const T_index min_N = amrex::max(NI1,NI2);
62 
63  T_PR * const AMREX_RESTRICT w1 = soa_1.m_rdata[PIdx::w];
64  T_PR * const AMREX_RESTRICT u1x = soa_1.m_rdata[PIdx::ux];
65  T_PR * const AMREX_RESTRICT u1y = soa_1.m_rdata[PIdx::uy];
66  T_PR * const AMREX_RESTRICT u1z = soa_1.m_rdata[PIdx::uz];
67 
68  T_PR * const AMREX_RESTRICT w2 = soa_2.m_rdata[PIdx::w];
69  T_PR * const AMREX_RESTRICT u2x = soa_2.m_rdata[PIdx::ux];
70  T_PR * const AMREX_RESTRICT u2y = soa_2.m_rdata[PIdx::uy];
71  T_PR * const AMREX_RESTRICT u2z = soa_2.m_rdata[PIdx::uz];
72 
73  // get local T1t and T2t
74  T_PR T1t; T_PR T2t;
75  if ( T1 <= T_PR(0.0) && L <= T_PR(0.0) )
76  {
77  T1t = ComputeTemperature(I1s,I1e,I1,u1x,u1y,u1z,m1);
78  }
79  else { T1t = T1; }
80  if ( T2 <= T_PR(0.0) && L <= T_PR(0.0) )
81  {
82  T2t = ComputeTemperature(I2s,I2e,I2,u2x,u2y,u2z,m2);
83  }
84  else { T2t = T2; }
85 
86  // local density
87  T_PR n1 = T_PR(0.0);
88  T_PR n2 = T_PR(0.0);
89  T_PR n12 = T_PR(0.0);
90  for (T_index i1=I1s; i1<I1e; ++i1) { n1 += w1[ I1[i1] ]; }
91  for (T_index i2=I2s; i2<I2e; ++i2) { n2 += w2[ I2[i2] ]; }
92  // Intra-species: the density is in fact the sum of the density of
93  // each sub-group of particles
94  if (isSameSpecies) {
95  n1 = n1 + n2;
96  n2 = n1;
97  }
98  n1 = n1 / dV; n2 = n2 / dV;
99  {
100  T_index i1 = I1s; T_index i2 = I2s;
101  for (T_index k = 0; k < max_N; ++k)
102  {
103  n12 += amrex::min( w1[ I1[i1] ], w2[ I2[i2] ] );
104  ++i1; if ( i1 == I1e ) { i1 = I1s; }
105  ++i2; if ( i2 == I2e ) { i2 = I2s; }
106  }
107  n12 = n12 / dV;
108  }
109  // Intra-species: Apply correction in collision rate
110  if (isSameSpecies) { n12 *= T_PR(2.0); }
111 
112  // compute Debye length lmdD
113  T_PR lmdD;
114  if ( T1t < T_PR(0.0) || T2t < T_PR(0.0) ) {
115  lmdD = T_PR(0.0);
116  }
117  else {
118  lmdD = T_PR(1.0)/std::sqrt( n1*q1*q1/(T1t*PhysConst::ep0) +
119  n2*q2*q2/(T2t*PhysConst::ep0) );
120  }
121  T_PR rmin = std::pow( T_PR(4.0) * MathConst::pi / T_PR(3.0) *
122  amrex::max(n1,n2), T_PR(-1.0/3.0) );
123  lmdD = amrex::max(lmdD, rmin);
124 
125 #if (defined WARPX_DIM_RZ)
126  T_PR * const AMREX_RESTRICT theta1 = soa_1.m_rdata[PIdx::theta];
127  T_PR * const AMREX_RESTRICT theta2 = soa_2.m_rdata[PIdx::theta];
128 #endif
129 
130  // call UpdateMomentumPerezElastic()
131  {
132  T_index i1 = I1s + coll_idx;
133  T_index i2 = I2s + coll_idx;
134 
135  // we will start from collision number = coll_idx and then add
136  // stride (smaller set size) until we do all collisions (larger set size)
137  for (T_index k = coll_idx; k < max_N; k += min_N)
138  {
139 #if (defined WARPX_DIM_RZ)
140  /* In RZ geometry, macroparticles can collide with other macroparticles
141  * in the same *cylindrical* cell. For this reason, collisions between macroparticles
142  * are actually not local in space. In this case, the underlying assumption is that
143  * particles within the same cylindrical cell represent a cylindrically-symmetry
144  * momentum distribution function. Therefore, here, we temporarily rotate the
145  * momentum of one of the macroparticles in agreement with this cylindrical symmetry.
146  * (This is technically only valid if we use only the m=0 azimuthal mode in the simulation;
147  * there is a corresponding assert statement at initialization.) */
148  T_PR const theta = theta2[I2[i2]]-theta1[I1[i1]];
149  T_PR const u1xbuf = u1x[I1[i1]];
150  u1x[I1[i1]] = u1xbuf*std::cos(theta) - u1y[I1[i1]]*std::sin(theta);
151  u1y[I1[i1]] = u1xbuf*std::sin(theta) + u1y[I1[i1]]*std::cos(theta);
152 #endif
153 
155  u1x[ I1[i1] ], u1y[ I1[i1] ], u1z[ I1[i1] ],
156  u2x[ I2[i2] ], u2y[ I2[i2] ], u2z[ I2[i2] ],
157  n1, n2, n12,
158  q1, m1, w1[ I1[i1] ], q2, m2, w2[ I2[i2] ],
159  dt, L, lmdD,
160  engine);
161 
162 #if (defined WARPX_DIM_RZ)
163  T_PR const u1xbuf_new = u1x[I1[i1]];
164  u1x[I1[i1]] = u1xbuf_new*std::cos(-theta) - u1y[I1[i1]]*std::sin(-theta);
165  u1y[I1[i1]] = u1xbuf_new*std::sin(-theta) + u1y[I1[i1]]*std::cos(-theta);
166 #endif
167 
168  if (max_N == NI1) {
169  i1 += min_N;
170  }
171  if (max_N == NI2) {
172  i2 += min_N;
173  }
174  }
175  }
176 
177 }
178 
179 #endif // WARPX_PARTICLES_COLLISION_ELASTIC_COLLISION_PEREZ_H_
#define AMREX_RESTRICT
#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 uz, T_R const m)
Definition: ComputeTemperature.H:15
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 const *AMREX_RESTRICT I1, T_index const *AMREX_RESTRICT I2, SoaData_type soa_1, SoaData_type soa_2, T_PR const q1, T_PR const q2, T_PR const m1, T_PR const m2, T_PR const T1, T_PR const T2, T_R const dt, T_PR const L, T_R const dV, amrex::RandomEngine const &engine, bool const isSameSpecies, T_index coll_idx)
Definition: ElasticCollisionPerez.H:45
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 ep0
vacuum permittivity: dielectric permittivity of vacuum [F/m]
Definition: constant.H:46
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
@ theta
RZ needs all three position components.
Definition: NamedComponentParticleContainer.H:36
@ uz
Definition: NamedComponentParticleContainer.H:34
@ w
weight
Definition: NamedComponentParticleContainer.H:33
@ uy
Definition: NamedComponentParticleContainer.H:34
@ ux
Definition: NamedComponentParticleContainer.H:34