WarpX
CollisionFilterFunc.H
Go to the documentation of this file.
1 /* Copyright 2023 The WarpX Community
2  *
3  * This file is part of WarpX.
4  *
5  * Authors: Roelof Groenewald (TAE Technologies), Neil Zaim
6  *
7  * License: BSD-3-Clause-LBNL
8  */
9 #ifndef WARPX_COLLISION_FILTER_FUNC_H_
10 #define WARPX_COLLISION_FILTER_FUNC_H_
11 
14 
15 #include <AMReX_Random.H>
16 
37 template <typename index_type>
39 void CollisionPairFilter (const amrex::ParticleReal u1x, const amrex::ParticleReal u1y,
40  const amrex::ParticleReal u1z, const amrex::ParticleReal u2x,
41  const amrex::ParticleReal u2y, const amrex::ParticleReal u2z,
42  const amrex::ParticleReal m1, const amrex::ParticleReal m2,
43  const amrex::ParticleReal w1, const amrex::ParticleReal w2,
44  const amrex::Real dt, const amrex::ParticleReal dV, const int pair_index,
45  index_type* AMREX_RESTRICT p_mask,
46  amrex::ParticleReal* AMREX_RESTRICT p_pair_reaction_weight,
47  const int multiplier,
48  const int process_count,
49  const ScatteringProcess::Executor* scattering_processes,
50  const amrex::RandomEngine& engine)
51 {
52  amrex::ParticleReal E_coll, v_coll, lab_to_COM_factor;
53 
55  u1x, u1y, u1z, u2x, u2y, u2z, m1, m2,
56  E_coll, v_coll, lab_to_COM_factor);
57 
58  using namespace amrex::literals;
59 
60  const amrex::ParticleReal w_min = amrex::min(w1, w2);
61  const amrex::ParticleReal w_max = amrex::max(w1, w2);
62 
63  // convert E_coll from Joule to eV
64  E_coll /= PhysConst::q_e;
65 
66  // Evaluate the cross-section for each scattering process to determine
67  // the total collision probability.
69  (process_count < 4), "Too many scattering processes in DSMC routine."
70  );
71  int coll_type[4] = {0, 0, 0, 0};
72  amrex::ParticleReal sigma_sums[4] = {0._prt, 0._prt, 0._prt, 0._prt};
73  for (int ii = 0; ii < process_count; ii++) {
74  auto const& scattering_process = scattering_processes[ii];
75  coll_type[ii] = int(scattering_process.m_type);
76  const amrex::ParticleReal sigma = scattering_process.getCrossSection(E_coll);
77  sigma_sums[ii] = sigma + ((ii == 0) ? 0._prt : sigma_sums[ii-1]);
78  }
79  const auto sigma_tot = sigma_sums[process_count-1];
80 
81  // calculate total collision probability
82  const amrex::ParticleReal exponent = (
83  lab_to_COM_factor * multiplier * w_max * sigma_tot * v_coll * dt / dV
84  );
85 
86  // Compute actual collision probability that is always between zero and one
87  // In principle this is obtained by computing 1 - exp(-probability_estimate)
88  // However, the computation of this quantity can fail numerically when probability_estimate is
89  // too small (e.g. exp(-probability_estimate) returns 1 and the computation returns 0).
90  // std::expm1 is used since it maintains correctness for small exponent.
91  const amrex::ParticleReal probability = -std::expm1(-exponent);
92 
93  // Now we determine if a collision should occur
94  if (amrex::Random(engine) < probability)
95  {
96  const amrex::ParticleReal random_number = amrex::Random(engine);
97  for (int ii = 0; ii < process_count; ii++) {
98  if (random_number <= sigma_sums[ii] / sigma_tot)
99  {
100  p_mask[pair_index] = coll_type[ii];
101  p_pair_reaction_weight[pair_index] = w_min;
102  break;
103  }
104  }
105  }
106  else
107  {
108  p_mask[pair_index] = false;
109  }
110 }
111 
112 #endif // WARPX_COLLISION_FILTER_FUNC_H_
#define AMREX_ASSERT_WITH_MESSAGE(EX, MSG)
#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
AMREX_GPU_HOST_DEVICE AMREX_INLINE void get_collision_parameters(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, amrex::ParticleReal &E_kin_COM, amrex::ParticleReal &v_rel_COM, amrex::ParticleReal &lab_to_COM_lorentz_factor)
Return (relativistic) collision energy, collision speed and Lorentz factor for transforming between t...
Definition: BinaryCollisionUtils.H:50
typename ParticleBins::index_type index_type
Definition: ParticleUtils.cpp:35
static constexpr auto q_e
elementary charge [C]
Definition: constant.H:50
Real Random()
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
const int[]
ii
Definition: check_interp_points_and_weights.py:148
float dt
Definition: stencil.py:442
Definition: ScatteringProcess.H:71