8 #ifndef WARPX_PROTON_BORON_FUSION_INITIALIZE_MOMENTUM_H
9 #define WARPX_PROTON_BORON_FUSION_INITIALIZE_MOMENTUM_H
60 void ProtonBoronFusionInitializeMomentum (
61 const SoaData_type& soa_1,
const SoaData_type& soa_2,
62 SoaData_type& soa_alpha,
65 const amrex::ParticleReal& m1,
const amrex::ParticleReal& m2,
72 using namespace amrex::literals;
74 constexpr amrex::ParticleReal mev_to_joule =
PhysConst::q_e*1.e6_prt;
77 constexpr amrex::ParticleReal E_fusion = 8.59009_prt*mev_to_joule;
80 constexpr amrex::ParticleReal E_decay = 0.0918984_prt*mev_to_joule;
88 constexpr
double mBe_sq = m_beryllium*m_beryllium;
93 amrex::ParticleReal ux_alpha1 = 0.0, uy_alpha1 = 0.0, uz_alpha1 = 0.0;
94 amrex::ParticleReal ux_Be = 0.0, uy_Be = 0.0, uz_Be = 0.0;
95 TwoProductFusionComputeProductMomenta (
102 ux_alpha1, uy_alpha1, uz_alpha1,
static_cast<amrex::ParticleReal
>(m_alpha),
103 ux_Be, uy_Be, uz_Be,
static_cast<amrex::ParticleReal
>(m_beryllium),
107 const amrex::ParticleReal px_Be =
static_cast<amrex::ParticleReal
>(m_beryllium) * ux_Be;
108 const amrex::ParticleReal py_Be =
static_cast<amrex::ParticleReal
>(m_beryllium) * uy_Be;
109 const amrex::ParticleReal pz_Be =
static_cast<amrex::ParticleReal
>(m_beryllium) * uz_Be;
113 constexpr amrex::ParticleReal gamma_Bestar = (1._prt + 0.5_prt*E_decay/(m_alpha*c_sq));
114 constexpr amrex::ParticleReal gamma_Bestar_sq_minus_one = gamma_Bestar*gamma_Bestar - 1._prt;
115 const amrex::ParticleReal p_Bestar =
static_cast<amrex::ParticleReal
>(m_alpha)*
PhysConst::c*std::sqrt(gamma_Bestar_sq_minus_one);
118 amrex::ParticleReal px_Bestar, py_Bestar, pz_Bestar;
122 amrex::ParticleReal px_alpha2, py_alpha2, pz_alpha2;
124 const amrex::ParticleReal p_Be_sq = px_Be*px_Be + py_Be*py_Be + pz_Be*pz_Be;
125 const amrex::ParticleReal g_Be = std::sqrt(1._prt + p_Be_sq / (
static_cast<amrex::ParticleReal
>(mBe_sq)*c_sq));
126 const amrex::ParticleReal mg_Be =
static_cast<amrex::ParticleReal
>(m_beryllium)*g_Be;
127 const amrex::ParticleReal v_Bex = px_Be / mg_Be;
128 const amrex::ParticleReal v_Bey = py_Be / mg_Be;
129 const amrex::ParticleReal v_Bez = pz_Be / mg_Be;
130 const amrex::ParticleReal v_Be_sq = v_Bex*v_Bex + v_Bey*v_Bey + v_Bez*v_Bez;
134 if ( v_Be_sq > std::numeric_limits<amrex::ParticleReal>::min() )
136 const amrex::ParticleReal vcDps = v_Bex*px_Bestar + v_Bey*py_Bestar + v_Bez*pz_Bestar;
137 const amrex::ParticleReal factor0 = (g_Be-1._prt)/v_Be_sq;
138 const amrex::ParticleReal factor = factor0*vcDps +
static_cast<amrex::ParticleReal
>(m_alpha)*gamma_Bestar*g_Be;
139 px_alpha2 = px_Bestar + v_Bex * factor;
140 py_alpha2 = py_Bestar + v_Bey * factor;
141 pz_alpha2 = pz_Bestar + v_Bez * factor;
145 px_alpha2 = px_Bestar;
146 py_alpha2 = py_Bestar;
147 pz_alpha2 = pz_Bestar;
151 const amrex::ParticleReal px_alpha3 = px_Be - px_alpha2;
152 const amrex::ParticleReal py_alpha3 = py_Be - py_alpha2;
153 const amrex::ParticleReal pz_alpha3 = pz_Be - pz_alpha2;
158 soa_alpha.m_rdata[
PIdx::ux][idx_alpha_start] = ux_alpha1;
159 soa_alpha.m_rdata[
PIdx::uy][idx_alpha_start] = uy_alpha1;
160 soa_alpha.m_rdata[
PIdx::uz][idx_alpha_start] = uz_alpha1;
161 soa_alpha.m_rdata[
PIdx::ux][idx_alpha_start + 1] = ux_alpha1;
162 soa_alpha.m_rdata[
PIdx::uy][idx_alpha_start + 1] = uy_alpha1;
163 soa_alpha.m_rdata[
PIdx::uz][idx_alpha_start + 1] = uz_alpha1;
164 soa_alpha.m_rdata[
PIdx::ux][idx_alpha_start + 2] = px_alpha2/
static_cast<amrex::ParticleReal
>(m_alpha);
165 soa_alpha.m_rdata[
PIdx::uy][idx_alpha_start + 2] = py_alpha2/
static_cast<amrex::ParticleReal
>(m_alpha);
166 soa_alpha.m_rdata[
PIdx::uz][idx_alpha_start + 2] = pz_alpha2/
static_cast<amrex::ParticleReal
>(m_alpha);
167 soa_alpha.m_rdata[
PIdx::ux][idx_alpha_start + 3] = px_alpha2/
static_cast<amrex::ParticleReal
>(m_alpha);
168 soa_alpha.m_rdata[
PIdx::uy][idx_alpha_start + 3] = py_alpha2/
static_cast<amrex::ParticleReal
>(m_alpha);
169 soa_alpha.m_rdata[
PIdx::uz][idx_alpha_start + 3] = pz_alpha2/
static_cast<amrex::ParticleReal
>(m_alpha);
170 soa_alpha.m_rdata[
PIdx::ux][idx_alpha_start + 4] = px_alpha3/
static_cast<amrex::ParticleReal
>(m_alpha);
171 soa_alpha.m_rdata[
PIdx::uy][idx_alpha_start + 4] = py_alpha3/
static_cast<amrex::ParticleReal
>(m_alpha);
172 soa_alpha.m_rdata[
PIdx::uz][idx_alpha_start + 4] = pz_alpha3/
static_cast<amrex::ParticleReal
>(m_alpha);
173 soa_alpha.m_rdata[
PIdx::ux][idx_alpha_start + 5] = px_alpha3/
static_cast<amrex::ParticleReal
>(m_alpha);
174 soa_alpha.m_rdata[
PIdx::uy][idx_alpha_start + 5] = py_alpha3/
static_cast<amrex::ParticleReal
>(m_alpha);
175 soa_alpha.m_rdata[
PIdx::uz][idx_alpha_start + 5] = pz_alpha3/
static_cast<amrex::ParticleReal
>(m_alpha);
#define AMREX_GPU_HOST_DEVICE
ParticleTile< ParticleType, NArrayReal, NArrayInt, Allocator > ParticleTileType
T_ParticleType ParticleType
AMREX_GPU_HOST_DEVICE AMREX_INLINE void RandomizeVelocity(amrex::ParticleReal &ux, amrex::ParticleReal &uy, amrex::ParticleReal &uz, const amrex::ParticleReal vp, amrex::RandomEngine const &engine)
Function to perform scattering of a particle that results in a random velocity vector with given magn...
Definition: ParticleUtils.H:162
typename WarpXParticleContainer::ParticleTileType ParticleTileType
Definition: ParticleUtils.cpp:32
typename ParticleTileType::ParticleTileDataType ParticleTileDataType
Definition: ParticleUtils.cpp:33
typename ParticleBins::index_type index_type
Definition: ParticleUtils.cpp:35
DenseBins< ParticleTileDataType > ParticleBins
Definition: ParticleUtils.cpp:34
typename WarpXParticleContainer::ParticleType ParticleType
Definition: ParticleUtils.cpp:31
static constexpr auto c
vacuum speed of light [m/s]
Definition: constant.H:44
static constexpr auto m_p
proton mass [kg]
Definition: constant.H:54
static constexpr auto m_u
dalton: unified atomic mass unit [kg]
Definition: constant.H:56
static constexpr auto q_e
elementary charge [C]
Definition: constant.H:50
@ uz
Definition: NamedComponentParticleContainer.H:34
@ uy
Definition: NamedComponentParticleContainer.H:34
@ ux
Definition: NamedComponentParticleContainer.H:34
ParticleTileData< StorageParticleType, NArrayReal, NArrayInt > ParticleTileDataType