WarpX
SingleNuclearFusionEvent.H
Go to the documentation of this file.
1 /* Copyright 2021 Neil Zaim
2  *
3  * This file is part of WarpX.
4  *
5  * License: BSD-3-Clause-LBNL
6  */
7 
8 #ifndef WARPX_SINGLE_NUCLEAR_FUSION_EVENT_H_
9 #define WARPX_SINGLE_NUCLEAR_FUSION_EVENT_H_
10 
13 
15 #include "Utils/WarpXConst.H"
16 
17 #include <AMReX_Algorithm.H>
18 #include <AMReX_Random.H>
19 #include <AMReX_REAL.H>
20 
21 #include <cmath>
22 
23 
52 template <typename index_type>
54 void SingleNuclearFusionEvent (const amrex::ParticleReal& u1x, const amrex::ParticleReal& u1y,
55  const amrex::ParticleReal& u1z, const amrex::ParticleReal& u2x,
56  const amrex::ParticleReal& u2y, const amrex::ParticleReal& u2z,
57  const amrex::ParticleReal& m1, const amrex::ParticleReal& m2,
58  amrex::ParticleReal w1, amrex::ParticleReal w2,
59  const amrex::Real& dt, const amrex::ParticleReal& dV, const int& pair_index,
60  index_type* AMREX_RESTRICT p_mask,
61  amrex::ParticleReal* AMREX_RESTRICT p_pair_reaction_weight,
62  const amrex::ParticleReal& fusion_multiplier,
63  const int& multiplier_ratio,
64  const amrex::ParticleReal& probability_threshold,
65  const amrex::ParticleReal& probability_target_value,
66  const NuclearFusionType& fusion_type,
67  const amrex::RandomEngine& engine)
68 {
69  amrex::ParticleReal E_coll, v_coll, lab_to_COM_factor;
70 
72  u1x, u1y, u1z, u2x, u2y, u2z, m1, m2,
73  E_coll, v_coll, lab_to_COM_factor);
74 
75  using namespace amrex::literals;
76 
77  const amrex::ParticleReal w_min = amrex::min(w1, w2);
78  const amrex::ParticleReal w_max = amrex::max(w1, w2);
79 
80  // Compute fusion cross section as a function of kinetic energy in the center of mass frame
81  auto fusion_cross_section = amrex::ParticleReal(0.);
82  if (fusion_type == NuclearFusionType::ProtonBoronToAlphas)
83  {
84  fusion_cross_section = ProtonBoronFusionCrossSection(E_coll);
85  }
89  {
90  fusion_cross_section = BoschHaleFusionCrossSection(E_coll, fusion_type, m1, m2);
91  }
92 
93  // First estimate of probability to have fusion reaction
94  amrex::ParticleReal probability_estimate = multiplier_ratio * fusion_multiplier *
95  lab_to_COM_factor * w_max * fusion_cross_section * v_coll * dt / dV;
96 
97  // Effective fusion multiplier
98  amrex::ParticleReal fusion_multiplier_eff = fusion_multiplier;
99 
100  // If the fusion probability is too high and the fusion multiplier greater than one, we risk to
101  // systematically underestimate the fusion yield. In this case, we reduce the fusion multiplier
102  // to reduce the fusion probability
103  if (probability_estimate > probability_threshold)
104  {
105  // We aim for a fusion probability of probability_target_value but take into account
106  // the constraint that the fusion_multiplier cannot be smaller than one
107  fusion_multiplier_eff = amrex::max(fusion_multiplier *
108  probability_target_value / probability_estimate , 1._prt);
109  probability_estimate *= fusion_multiplier_eff/fusion_multiplier;
110  }
111 
112  // Compute actual fusion probability that is always between zero and one
113  // In principle this is obtained by computing 1 - exp(-probability_estimate)
114  // However, the computation of this quantity can fail numerically when probability_estimate is
115  // too small (e.g. exp(-probability_estimate) returns 1 and the computation returns 0).
116  // std::expm1 is used since it maintains correctness for small exponent.
117  const amrex::ParticleReal probability = -std::expm1(-probability_estimate);
118 
119  // Get a random number
120  const amrex::ParticleReal random_number = amrex::Random(engine);
121 
122  // If we have a fusion event, set the mask the true and fill the product weight array
123  if (random_number < probability)
124  {
125  p_mask[pair_index] = true;
126  p_pair_reaction_weight[pair_index] = w_min/fusion_multiplier_eff;
127  }
128  else
129  {
130  p_mask[pair_index] = false;
131  }
132 
133 }
134 #endif // WARPX_SINGLE_NUCLEAR_FUSION_EVENT_H_
#define AMREX_INLINE
#define AMREX_GPU_HOST_DEVICE
NuclearFusionType
Definition: BinaryCollisionUtils.H:26
@ DeuteriumDeuteriumToProtonTritium
@ DeuteriumDeuteriumToNeutronHelium
AMREX_GPU_HOST_DEVICE AMREX_INLINE amrex::ParticleReal BoschHaleFusionCrossSection(const amrex::ParticleReal &E_kin_star, const NuclearFusionType &fusion_type, const amrex::ParticleReal &m1, const amrex::ParticleReal &m2)
Computes the fusion cross section, using the analytical fits given in H.-S. Bosch and G....
Definition: BoschHaleFusionCrossSection.H:29
AMREX_GPU_HOST_DEVICE AMREX_INLINE amrex::ParticleReal ProtonBoronFusionCrossSection(const amrex::ParticleReal &E_kin_star)
Computes the total proton-boron fusion cross section. When E_kin_star < 3.5 MeV, we use the analytica...
Definition: ProtonBoronFusionCrossSection.H:139
AMREX_GPU_HOST_DEVICE AMREX_INLINE void SingleNuclearFusionEvent(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 w1, 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 amrex::ParticleReal &fusion_multiplier, const int &multiplier_ratio, const amrex::ParticleReal &probability_threshold, const amrex::ParticleReal &probability_target_value, const NuclearFusionType &fusion_type, const amrex::RandomEngine &engine)
This function computes whether the collision between two particles result in a nuclear fusion event,...
Definition: SingleNuclearFusionEvent.H:54
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
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
float dt
Definition: stencil.py:442