diff options
Diffstat (limited to 'Source/Particles/Collision/BinaryCollision')
3 files changed, 19 insertions, 15 deletions
diff --git a/Source/Particles/Collision/BinaryCollision/NuclearFusion/ProtonBoronFusionCrossSection.H b/Source/Particles/Collision/BinaryCollision/NuclearFusion/ProtonBoronFusionCrossSection.H index 6cd4385f3..4a5d518b7 100644 --- a/Source/Particles/Collision/BinaryCollision/NuclearFusion/ProtonBoronFusionCrossSection.H +++ b/Source/Particles/Collision/BinaryCollision/NuclearFusion/ProtonBoronFusionCrossSection.H @@ -10,6 +10,7 @@ #include "Utils/WarpXConst.H" +#include <AMReX_Math.H> #include <AMReX_REAL.H> #include <cmath> @@ -29,6 +30,7 @@ AMREX_GPU_HOST_DEVICE AMREX_INLINE amrex::ParticleReal ProtonBoronFusionCrossSectionNevins (const amrex::ParticleReal& E_keV) { using namespace amrex::literals; + using namespace amrex::Math; // If kinetic energy is 0, return a 0 cross section and avoid later division by 0. if (E_keV == 0._prt) {return 0._prt;} @@ -61,7 +63,7 @@ amrex::ParticleReal ProtonBoronFusionCrossSectionNevins (const amrex::ParticleRe constexpr auto AL = 1.82e4_prt; constexpr auto EL = 148._prt; constexpr auto dEL_sq = 2.35_prt*2.35_prt; - astrophysical_factor = C0 + C1*E_keV + C2*E_keV*E_keV + + astrophysical_factor = C0 + C1*E_keV + C2*powi<2>(E_keV) + AL/((E_keV - EL)*(E_keV - EL) + dEL_sq); } else if (E_keV < E_lim2) @@ -71,7 +73,7 @@ amrex::ParticleReal ProtonBoronFusionCrossSectionNevins (const amrex::ParticleRe constexpr auto D2 = -20.3_prt; constexpr auto D5 = -1.58_prt; const amrex::ParticleReal E_norm = (E_keV-400._prt) * 1.e-2_prt; - astrophysical_factor = D0 + D1*E_norm + D2*E_norm*E_norm + D5*std::pow(E_norm, 5_prt); + astrophysical_factor = D0 + D1*E_norm + D2*powi<2>(E_norm) + D5*powi<5>(E_norm); } else { diff --git a/Source/Particles/Collision/BinaryCollision/NuclearFusion/SingleNuclearFusionEvent.H b/Source/Particles/Collision/BinaryCollision/NuclearFusion/SingleNuclearFusionEvent.H index e4704f746..780aa1f76 100644 --- a/Source/Particles/Collision/BinaryCollision/NuclearFusion/SingleNuclearFusionEvent.H +++ b/Source/Particles/Collision/BinaryCollision/NuclearFusion/SingleNuclearFusionEvent.H @@ -15,6 +15,7 @@ #include "Utils/WarpXConst.H" #include <AMReX_Algorithm.H> +#include <AMReX_Math.H> #include <AMReX_Random.H> #include <AMReX_REAL.H> @@ -71,6 +72,7 @@ void SingleNuclearFusionEvent (const amrex::ParticleReal& u1x, const amrex::Part // x_star denotes the value of x in the center of mass frame using namespace amrex::literals; + using namespace amrex::Math; const amrex::ParticleReal w_min = amrex::min(w1, w2); const amrex::ParticleReal w_max = amrex::max(w1, w2); @@ -95,10 +97,9 @@ void SingleNuclearFusionEvent (const amrex::ParticleReal& u1x, const amrex::Part const amrex::ParticleReal p2y = u2y * m2; const amrex::ParticleReal p2z = u2z * m2; // Square norm of the total (sum between the two particles) momenta in the lab frame - auto constexpr pow2 = [](double const x) { return x*x; }; - const amrex::ParticleReal p_total_sq = pow2(p1x + p2x) + - pow2(p1y+p2y) + - pow2(p1z+p2z); + const amrex::ParticleReal p_total_sq = powi<2>(p1x + p2x) + + powi<2>(p1y+p2y) + + powi<2>(p1z+p2z); // Total energy in the lab frame const amrex::ParticleReal E_lab = (m1 * g1 + m2 * g2) * c_sq; @@ -128,8 +129,8 @@ void SingleNuclearFusionEvent (const amrex::ParticleReal& u1x, const amrex::Part // 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 ); + const amrex::ParticleReal p_star_sq = m1*m2*c_sq * ( powi<2>(E_ratio) - one_pr ) + + powi<2>(m1 - m2)*c_sq*inv_four_pr * powi<2>( E_ratio - 1._prt/E_ratio ); // Lorentz factors in the center of mass frame const amrex::ParticleReal g1_star = std::sqrt(one_pr + p_star_sq / (m1_sq*c_sq)); diff --git a/Source/Particles/Collision/BinaryCollision/NuclearFusion/TwoProductFusionUtil.H b/Source/Particles/Collision/BinaryCollision/NuclearFusion/TwoProductFusionUtil.H index e66425e2d..09afe1e13 100644 --- a/Source/Particles/Collision/BinaryCollision/NuclearFusion/TwoProductFusionUtil.H +++ b/Source/Particles/Collision/BinaryCollision/NuclearFusion/TwoProductFusionUtil.H @@ -11,6 +11,7 @@ #include "Utils/ParticleUtils.H" #include "Utils/WarpXConst.H" +#include <AMReX_Math.H> #include <AMReX_Random.H> #include <AMReX_REAL.H> @@ -65,6 +66,7 @@ namespace { const amrex::RandomEngine& engine ) { using namespace amrex::literals; + using namespace amrex::Math; constexpr amrex::ParticleReal c_sq = PhysConst::c * PhysConst::c; constexpr amrex::ParticleReal inv_csq = 1._prt / ( c_sq ); @@ -87,10 +89,9 @@ namespace { const amrex::ParticleReal p2y_in = u2y_in * m2_in; const amrex::ParticleReal p2z_in = u2z_in * m2_in; // Square norm of the total (sum between the two particles) momenta in the lab frame - auto constexpr pow2 = [](double const x) { return x*x; }; - const amrex::ParticleReal p_total_sq = pow2(p1x_in+p2x_in) + - pow2(p1y_in+p2y_in) + - pow2(p1z_in+p2z_in); + const amrex::ParticleReal p_total_sq = powi<2>(p1x_in+p2x_in) + + powi<2>(p1y_in+p2y_in) + + powi<2>(p1z_in+p2z_in); // Total energy of incident macroparticles in the lab frame const amrex::ParticleReal E_lab = (m1_in * g1_in + m2_in * g2_in) * c_sq; @@ -103,7 +104,7 @@ namespace { // due to possible inconsistencies in how the mass is defined in the code, it is // probably more robust to subtract the rest masses and to add the fusion energy to the // total kinetic energy. - const amrex::ParticleReal E_star_f_sq = pow2(std::sqrt(E_star_sq) + const amrex::ParticleReal E_star_f_sq = powi<2>(std::sqrt(E_star_sq) - E_rest_in + E_rest_out + E_fusion); // Square of the norm of the momentum of the products in the center of mass frame @@ -111,8 +112,8 @@ namespace { // 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 ); + const amrex::ParticleReal p_star_f_sq = m1_out*m2_out*c_sq * ( powi<2>(E_ratio) - 1._prt ) + + powi<2>(m1_out - m2_out)*c_sq*0.25_prt * powi<2>( E_ratio - 1._prt/E_ratio ); // Compute momentum of first product in the center of mass frame, assuming isotropic // distribution |