diff options
Diffstat (limited to 'Source/Particles/Collision/BinaryCollision')
3 files changed, 20 insertions, 14 deletions
diff --git a/Source/Particles/Collision/BinaryCollision/NuclearFusion/ProtonBoronFusionCrossSection.H b/Source/Particles/Collision/BinaryCollision/NuclearFusion/ProtonBoronFusionCrossSection.H index ee81a620a..6cd4385f3 100644 --- a/Source/Particles/Collision/BinaryCollision/NuclearFusion/ProtonBoronFusionCrossSection.H +++ b/Source/Particles/Collision/BinaryCollision/NuclearFusion/ProtonBoronFusionCrossSection.H @@ -71,7 +71,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); + astrophysical_factor = D0 + D1*E_norm + D2*E_norm*E_norm + D5*std::pow(E_norm, 5_prt); } else { @@ -115,7 +115,7 @@ amrex::ParticleReal ProtonBoronFusionCrossSectionBuck (const amrex::ParticleReal constexpr amrex::ParticleReal E_start_fit = 3500._prt; // Fit starts at 3.5 MeV // cross section at E = E_start_fit, in barn constexpr amrex::ParticleReal cross_section_start_fit = 0.2168440845211521_prt; - constexpr amrex::ParticleReal slope_fit = -2.661840717596765; + constexpr amrex::ParticleReal slope_fit = -2.661840717596765_prt; // Compute fitted value return cross_section_start_fit*std::pow(E_keV/E_start_fit, slope_fit); diff --git a/Source/Particles/Collision/BinaryCollision/NuclearFusion/ProtonBoronFusionInitializeMomentum.H b/Source/Particles/Collision/BinaryCollision/NuclearFusion/ProtonBoronFusionInitializeMomentum.H index 415555193..2fe2abb92 100644 --- a/Source/Particles/Collision/BinaryCollision/NuclearFusion/ProtonBoronFusionInitializeMomentum.H +++ b/Source/Particles/Collision/BinaryCollision/NuclearFusion/ProtonBoronFusionInitializeMomentum.H @@ -118,9 +118,10 @@ namespace { 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 - const amrex::ParticleReal p_total_sq = std::pow(p1x+p2x, 2) + - std::pow(p1y+p2y, 2) + - std::pow(p1z+p2z, 2); + auto constexpr pow2 = [](double const x) { return x*x; }; + const amrex::ParticleReal p_total_sq = pow2(p1x+p2x) + + pow2(p1y+p2y) + + pow2(p1z+p2z); // Total energy of proton+boron in the lab frame const amrex::ParticleReal E_lab = (m1 * g1 + m2 * g2) * c_sq; @@ -135,14 +136,15 @@ namespace { // mass of the other elements is the atomic mass, that includes the electrons mass), 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 = std::pow(std::sqrt(E_star_sq) - - E_rest_pb + E_rest_abe + E_fusion, 2); + const amrex::ParticleReal E_star_f_sq = pow2(std::sqrt(E_star_sq) + - E_rest_pb + E_rest_abe + E_fusion); // Square of the norm of the momentum of Beryllium or alpha 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 - (ma_sq + mBe_sq)*c_sq*0.5_prt + - std::pow(c_sq, 3)*0.25_prt*std::pow(ma_sq - mBe_sq, 2)/E_star_f_sq; + pow3(c_sq)*0.25_prt * pow2(ma_sq - mBe_sq) / E_star_f_sq; // Compute momentum of first alpha in the center of mass frame, assuming isotropic // distribution diff --git a/Source/Particles/Collision/BinaryCollision/NuclearFusion/SingleNuclearFusionEvent.H b/Source/Particles/Collision/BinaryCollision/NuclearFusion/SingleNuclearFusionEvent.H index c8633b52b..c59ad1a12 100644 --- a/Source/Particles/Collision/BinaryCollision/NuclearFusion/SingleNuclearFusionEvent.H +++ b/Source/Particles/Collision/BinaryCollision/NuclearFusion/SingleNuclearFusionEvent.H @@ -69,6 +69,8 @@ void SingleNuclearFusionEvent (const amrex::ParticleReal& u1x, const amrex::Part // x_sq denotes the square of x // x_star denotes the value of x in the center of mass frame + using namespace amrex::literals; + const amrex::ParticleReal w_min = amrex::min(w1, w2); const amrex::ParticleReal w_max = amrex::max(w1, w2); @@ -93,9 +95,10 @@ 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 - const amrex::ParticleReal p_total_sq = std::pow(p1x+p2x, 2) + - std::pow(p1y+p2y, 2) + - std::pow(p1z+p2z, 2); + auto constexpr pow2 = [](double const x) { return x*x; }; + const amrex::ParticleReal p_total_sq = pow2(p1x + p2x) + + pow2(p1y+p2y) + + pow2(p1z+p2z); // Total energy in the lab frame const amrex::ParticleReal E_lab = (m1 * g1 + m2 * g2) * c_sq; @@ -115,13 +118,14 @@ 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 + - std::pow(c_sq, 3)*inv_four_pr*std::pow(m1_sq - m2_sq, 2)/E_star_sq; + pow3(c_sq) * inv_four_pr * pow2(m1_sq - m2_sq) / E_star_sq; // Lorentz factors in the center of mass frame - const amrex::ParticleReal g1_star = std::sqrt(1 + p_star_sq / (m1_sq*c_sq)); - const amrex::ParticleReal g2_star = std::sqrt(1 + p_star_sq / (m2_sq*c_sq)); + 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)); // relative velocity in the center of mass frame const amrex::ParticleReal v_rel = std::sqrt(p_star_sq) * (one_pr/(m1*g1_star) + |