WarpX
DSMCFunc.H
Go to the documentation of this file.
1 /* Copyright 2024 The WarpX Community
2  *
3  * This file is part of WarpX.
4  *
5  * Authors: Roelof Groenewald (TAE Technologies)
6  *
7  * License: BSD-3-Clause-LBNL
8  */
9 
10 #ifndef WARPX_DSMC_FUNC_H_
11 #define WARPX_DSMC_FUNC_H_
12 
13 #include "CollisionFilterFunc.H"
14 
24 #include "Utils/ParticleUtils.H"
26 
27 #include <AMReX_DenseBins.H>
28 #include <AMReX_ParmParse.H>
29 #include <AMReX_Random.H>
30 
37 class DSMCFunc
38 {
39  // Define shortcuts for frequently-used type names
46 
47 public:
48 
52  DSMCFunc () = default;
53 
61  DSMCFunc ( const std::string& collision_name,
62  MultiParticleContainer const * mypc,
63  bool isSameSpecies );
64 
65  struct Executor {
93  index_type const I1s, index_type const I1e,
94  index_type const I2s, index_type const I2e,
95  index_type const* AMREX_RESTRICT I1,
96  index_type const* AMREX_RESTRICT I2,
97  const SoaData_type& soa_1, const SoaData_type& soa_2,
98  GetParticlePosition<PIdx> /*get_position_1*/, GetParticlePosition<PIdx> /*get_position_2*/,
99  amrex::ParticleReal const /*q1*/, amrex::ParticleReal const /*q2*/,
100  amrex::ParticleReal const m1, amrex::ParticleReal const m2,
101  amrex::Real const dt, amrex::Real const dV, index_type coll_idx,
102  index_type const cell_start_pair, index_type* AMREX_RESTRICT p_mask,
103  index_type* AMREX_RESTRICT p_pair_indices_1, index_type* AMREX_RESTRICT p_pair_indices_2,
104  amrex::ParticleReal* AMREX_RESTRICT p_pair_reaction_weight,
105  amrex::RandomEngine const& engine) const
106  {
107  amrex::ParticleReal * const AMREX_RESTRICT w1 = soa_1.m_rdata[PIdx::w];
108  amrex::ParticleReal * const AMREX_RESTRICT u1x = soa_1.m_rdata[PIdx::ux];
109  amrex::ParticleReal * const AMREX_RESTRICT u1y = soa_1.m_rdata[PIdx::uy];
110  amrex::ParticleReal * const AMREX_RESTRICT u1z = soa_1.m_rdata[PIdx::uz];
111 
112  amrex::ParticleReal * const AMREX_RESTRICT w2 = soa_2.m_rdata[PIdx::w];
113  amrex::ParticleReal * const AMREX_RESTRICT u2x = soa_2.m_rdata[PIdx::ux];
114  amrex::ParticleReal * const AMREX_RESTRICT u2y = soa_2.m_rdata[PIdx::uy];
115  amrex::ParticleReal * const AMREX_RESTRICT u2z = soa_2.m_rdata[PIdx::uz];
116 
117  // Number of macroparticles of each species
118  const index_type NI1 = I1e - I1s;
119  const index_type NI2 = I2e - I2s;
120  const index_type max_N = amrex::max(NI1,NI2);
121  const index_type min_N = amrex::min(NI1,NI2);
122 
123  index_type pair_index = cell_start_pair + coll_idx;
124 
125  // Because the number of particles of each species is not always equal (NI1 != NI2
126  // in general), some macroparticles will be paired with multiple macroparticles of the
127  // other species and we need to decrease their weight accordingly.
128  // c1 corresponds to the minimum number of times a particle of species 1 will be paired
129  // with a particle of species 2. Same for c2.
130  // index_type(1): https://github.com/AMReX-Codes/amrex/pull/3684
131  const index_type c1 = amrex::max(NI2/NI1, index_type(1));
132  const index_type c2 = amrex::max(NI1/NI2, index_type(1));
133 
134 #if (defined WARPX_DIM_RZ)
135  amrex::ParticleReal * const AMREX_RESTRICT theta1 = soa_1.m_rdata[PIdx::theta];
136  amrex::ParticleReal * const AMREX_RESTRICT theta2 = soa_2.m_rdata[PIdx::theta];
137 #endif
138  index_type i1 = I1s + coll_idx;
139  index_type i2 = I2s + coll_idx;
140  // we will start from collision number = coll_idx and then add
141  // stride (smaller set size) until we do all collisions (larger set size)
142  for (index_type k = coll_idx; k < max_N; k += min_N)
143  {
144  // c1k : how many times the current particle of species 1 is paired with a particle
145  // of species 2. Same for c2k.
146  const index_type c1k = (k%NI1 < max_N%NI1) ? c1 + 1: c1;
147  const index_type c2k = (k%NI2 < max_N%NI2) ? c2 + 1: c2;
148 
149 #if (defined WARPX_DIM_RZ)
150  /* In RZ geometry, macroparticles can collide with other macroparticles
151  * in the same *cylindrical* cell. For this reason, collisions between macroparticles
152  * are actually not local in space. In this case, the underlying assumption is that
153  * particles within the same cylindrical cell represent a cylindrically-symmetry
154  * momentum distribution function. Therefore, here, we temporarily rotate the
155  * momentum of one of the macroparticles in agreement with this cylindrical symmetry.
156  * (This is technically only valid if we use only the m=0 azimuthal mode in the simulation;
157  * there is a corresponding assert statement at initialization.)
158  */
159  amrex::ParticleReal const theta = theta2[I2[i2]]-theta1[I1[i1]];
160  amrex::ParticleReal const u1xbuf = u1x[I1[i1]];
161  u1x[I1[i1]] = u1xbuf*std::cos(theta) - u1y[I1[i1]]*std::sin(theta);
162  u1y[I1[i1]] = u1xbuf*std::sin(theta) + u1y[I1[i1]]*std::cos(theta);
163 #endif
164 
166  u1x[ I1[i1] ], u1y[ I1[i1] ], u1z[ I1[i1] ],
167  u2x[ I2[i2] ], u2y[ I2[i2] ], u2z[ I2[i2] ],
168  m1, m2, w1[ I1[i1] ]/c1k, w2[ I2[i2] ]/c2k,
169  dt, dV, static_cast<int>(pair_index), p_mask,
170  p_pair_reaction_weight, static_cast<int>(max_N),
172 
173 #if (defined WARPX_DIM_RZ)
174  amrex::ParticleReal const u1xbuf_new = u1x[I1[i1]];
175  u1x[I1[i1]] = u1xbuf_new*std::cos(-theta) - u1y[I1[i1]]*std::sin(-theta);
176  u1y[I1[i1]] = u1xbuf_new*std::sin(-theta) + u1y[I1[i1]]*std::cos(-theta);
177 #endif
178 
179  p_pair_indices_1[pair_index] = I1[i1];
180  p_pair_indices_2[pair_index] = I2[i2];
181  if (max_N == NI1) {
182  i1 += min_N;
183  }
184  if (max_N == NI2) {
185  i2 += min_N;
186  }
187  pair_index += min_N;
188  }
189  }
190 
193  };
194 
195  [[nodiscard]] Executor const& executor () const { return m_exe; }
196 
197 private:
200 
202 };
203 
204 #endif // WARPX_DSMC_FUNC_H_
#define AMREX_RESTRICT
#define AMREX_INLINE
#define AMREX_GPU_HOST_DEVICE
AMREX_GPU_HOST_DEVICE AMREX_INLINE void CollisionPairFilter(const amrex::ParticleReal u1x, const amrex::ParticleReal u1y, const amrex::ParticleReal u1z, const amrex::ParticleReal u2x, const amrex::ParticleReal u2y, const amrex::ParticleReal u2z, const amrex::ParticleReal m1, const amrex::ParticleReal m2, const amrex::ParticleReal w1, const amrex::ParticleReal w2, const amrex::Real dt, const amrex::ParticleReal dV, const int pair_index, index_type *AMREX_RESTRICT p_mask, amrex::ParticleReal *AMREX_RESTRICT p_pair_reaction_weight, const int multiplier, const int process_count, const ScatteringProcess::Executor *scattering_processes, const amrex::RandomEngine &engine)
This function determines whether a collision occurs for a given pair of particles.
Definition: CollisionFilterFunc.H:39
This class performs DSMC (direct simulation Monte Carlo) collisions within a cell....
Definition: DSMCFunc.H:38
amrex::Vector< ScatteringProcess > m_scattering_processes
Definition: DSMCFunc.H:198
ParticleBins::index_type index_type
Definition: DSMCFunc.H:44
Executor m_exe
Definition: DSMCFunc.H:201
DSMCFunc()=default
Constructor of the DSMCFunc class.
DSMCFunc(const std::string &collision_name, MultiParticleContainer const *mypc, bool isSameSpecies)
Constructor of the DSMCFunc class.
amrex::Gpu::DeviceVector< ScatteringProcess::Executor > m_scattering_processes_exe
Definition: DSMCFunc.H:199
WarpXParticleContainer::ParticleType ParticleType
Definition: DSMCFunc.H:40
Executor const & executor() const
Definition: DSMCFunc.H:195
Definition: MultiParticleContainer.H:66
unsigned int index_type
ParticleTile< ParticleType, NArrayReal, NArrayInt, Allocator > ParticleTileType
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
Definition: DSMCFunc.H:65
int m_process_count
Definition: DSMCFunc.H:191
ScatteringProcess::Executor * m_scattering_processes_data
Definition: DSMCFunc.H:192
AMREX_GPU_HOST_DEVICE AMREX_INLINE void operator()(index_type const I1s, index_type const I1e, index_type const I2s, index_type const I2e, index_type const *AMREX_RESTRICT I1, index_type const *AMREX_RESTRICT I2, const SoaData_type &soa_1, const SoaData_type &soa_2, GetParticlePosition< PIdx >, GetParticlePosition< PIdx >, amrex::ParticleReal const, amrex::ParticleReal const, amrex::ParticleReal const m1, amrex::ParticleReal const m2, amrex::Real const dt, amrex::Real const dV, index_type coll_idx, index_type const cell_start_pair, index_type *AMREX_RESTRICT p_mask, index_type *AMREX_RESTRICT p_pair_indices_1, index_type *AMREX_RESTRICT p_pair_indices_2, amrex::ParticleReal *AMREX_RESTRICT p_pair_reaction_weight, amrex::RandomEngine const &engine) const
Executor of the DSMCFunc class. Performs DSMC collisions at the cell level. Note that this function d...
Definition: DSMCFunc.H:92
@ 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
Definition: ScatteringProcess.H:71
GpuArray< ParticleReal *, NAR > m_rdata
ParticleTileData< StorageParticleType, NArrayReal, NArrayInt > ParticleTileDataType