diff options
Diffstat (limited to 'Source/Particles/Collision/BinaryCollision')
| -rw-r--r-- | Source/Particles/Collision/BinaryCollision/NuclearFusion/SingleNuclearFusionEvent.H | 17 | ||||
| -rw-r--r-- | Source/Particles/Collision/BinaryCollision/NuclearFusion/TwoProductFusionUtil.H | 9 |
2 files changed, 14 insertions, 12 deletions
diff --git a/Source/Particles/Collision/BinaryCollision/NuclearFusion/SingleNuclearFusionEvent.H b/Source/Particles/Collision/BinaryCollision/NuclearFusion/SingleNuclearFusionEvent.H index 30c054bc4..8b3bd625b 100644 --- a/Source/Particles/Collision/BinaryCollision/NuclearFusion/SingleNuclearFusionEvent.H +++ b/Source/Particles/Collision/BinaryCollision/NuclearFusion/SingleNuclearFusionEvent.H @@ -75,7 +75,6 @@ void SingleNuclearFusionEvent (const amrex::ParticleReal& u1x, const amrex::Part const amrex::ParticleReal w_max = amrex::max(w1, w2); constexpr auto one_pr = amrex::ParticleReal(1.); - constexpr auto inv_two_pr = amrex::ParticleReal(1./2.); constexpr auto inv_four_pr = amrex::ParticleReal(1./4.); constexpr amrex::ParticleReal c_sq = PhysConst::c * PhysConst::c; constexpr amrex::ParticleReal inv_csq = one_pr / ( c_sq ); @@ -107,7 +106,8 @@ void SingleNuclearFusionEvent (const amrex::ParticleReal& u1x, const amrex::Part const amrex::ParticleReal E_star_sq = E_lab*E_lab - c_sq*p_total_sq; // Kinetic energy in the center of mass frame - const amrex::ParticleReal E_kin_star = std::sqrt(E_star_sq) - (m1 + m2)*c_sq; + const amrex::ParticleReal E_star = std::sqrt(E_star_sq); + const amrex::ParticleReal E_kin_star = E_star - (m1 + m2)*c_sq; // Compute fusion cross section as a function of kinetic energy in the center of mass frame auto fusion_cross_section = amrex::ParticleReal(0.); @@ -118,14 +118,15 @@ void SingleNuclearFusionEvent (const amrex::ParticleReal& u1x, const amrex::Part // Square of the norm of the momentum of one of the particles in the center of mass frame // Formula obtained by inverting E^2 = p^2*c^2 + m^2*c^4 in the COM frame for each particle - auto constexpr pow3 = [](double const x) { return x*x*x; }; - const amrex::ParticleReal p_star_sq = - E_star_sq*inv_four_pr*inv_csq - (m1_sq + m2_sq)*c_sq*inv_two_pr + - pow3(c_sq) * inv_four_pr * pow2(m1_sq - m2_sq) / E_star_sq; + // The expression below is specifically written in a form that avoids returning + // small negative numbers due to machine precision errors, for low-energy particles + const amrex::ParticleReal E_ratio = E_star/((m1 + m2)*c_sq); + const amrex::ParticleReal p_star_sq = m1*m2*c_sq * ( pow2(E_ratio) - one_pr ) + + pow2(m1 - m2)*c_sq*inv_four_pr * pow2( E_ratio - 1._prt/E_ratio ); // Lorentz factors in the center of mass frame - const amrex::ParticleReal g1_star = std::sqrt(1._prt + p_star_sq / (m1_sq*c_sq)); - const amrex::ParticleReal g2_star = std::sqrt(1._prt + p_star_sq / (m2_sq*c_sq)); + const amrex::ParticleReal g1_star = std::sqrt(one_pr + p_star_sq / (m1_sq*c_sq)); + const amrex::ParticleReal g2_star = std::sqrt(one_pr + p_star_sq / (m2_sq*c_sq)); // relative velocity in the center of mass frame const amrex::ParticleReal v_rel = std::sqrt(p_star_sq) * (one_pr/(m1*g1_star) + diff --git a/Source/Particles/Collision/BinaryCollision/NuclearFusion/TwoProductFusionUtil.H b/Source/Particles/Collision/BinaryCollision/NuclearFusion/TwoProductFusionUtil.H index a440109d0..e66425e2d 100644 --- a/Source/Particles/Collision/BinaryCollision/NuclearFusion/TwoProductFusionUtil.H +++ b/Source/Particles/Collision/BinaryCollision/NuclearFusion/TwoProductFusionUtil.H @@ -108,10 +108,11 @@ namespace { // Square of the norm of the momentum of the products in the center of mass frame // Formula obtained by inverting E^2 = p^2*c^2 + m^2*c^4 in the COM frame for each particle - auto constexpr pow3 = [](double const x) { return x*x*x; }; - const amrex::ParticleReal p_star_f_sq = - E_star_f_sq*0.25_prt*inv_csq - (m1_out*m1_out + m2_out*m2_out)*c_sq*0.5_prt + - pow3(c_sq)*0.25_prt * pow2(m1_out*m1_out - m2_out*m2_out) / E_star_f_sq; + // The expression below is specifically written in a form that avoids returning + // small negative numbers due to machine precision errors, for low-energy particles + const amrex::ParticleReal E_ratio = std::sqrt(E_star_f_sq)/((m1_out + m2_out)*c_sq); + const amrex::ParticleReal p_star_f_sq = m1_out*m2_out*c_sq * ( pow2(E_ratio) - 1._prt ) + + pow2(m1_out - m2_out)*c_sq*0.25_prt * pow2( E_ratio - 1._prt/E_ratio ); // Compute momentum of first product in the center of mass frame, assuming isotropic // distribution |