diff options
6 files changed, 55 insertions, 48 deletions
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmComoving.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmComoving.cpp index 8f576dad7..aecc754d4 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmComoving.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmComoving.cpp @@ -13,6 +13,7 @@ #include <AMReX_GpuComplex.H> #include <AMReX_GpuLaunch.H> #include <AMReX_GpuQualifiers.H> +#include <AMReX_Math.H> #include <AMReX_MFIter.H> #include <AMReX_PODVector.H> @@ -198,21 +199,21 @@ void PsatdAlgorithmComoving::InitializeSpectralCoefficients (const SpectralKSpac { // Calculate norm of finite-order k vector const amrex::Real knorm_mod = std::sqrt( - std::pow(kx_mod[i], 2) + + amrex::Math::powi<2>(kx_mod[i]) + #if defined(WARPX_DIM_3D) - std::pow(ky_mod[j], 2) + - std::pow(kz_mod[k], 2)); + amrex::Math::powi<2>(ky_mod[j]) + + amrex::Math::powi<2>(kz_mod[k])); #else - std::pow(kz_mod[j], 2)); + amrex::Math::powi<2>(kz_mod[j])); #endif // Calculate norm of infinite-order k vector const amrex::Real knorm = std::sqrt( - std::pow(kx[i], 2) + + amrex::Math::powi<2>(kx[i]) + #if defined(WARPX_DIM_3D) - std::pow(ky[j], 2) + - std::pow(kz[k], 2)); + amrex::Math::powi<2>(ky[j]) + + amrex::Math::powi<2>(kz[k])); #else - std::pow(kz[j], 2)); + amrex::Math::powi<2>(kz[j])); #endif // Physical constants c, c**2, and epsilon_0, and imaginary unit constexpr amrex::Real c = PhysConst::c; diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmJConstantInTime.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmJConstantInTime.cpp index 8c28c42b2..ae1cf65eb 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmJConstantInTime.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmJConstantInTime.cpp @@ -19,6 +19,7 @@ #include <AMReX_GpuLaunch.H> #include <AMReX_GpuQualifiers.H> #include <AMReX_IntVect.H> +#include <AMReX_Math.H> #include <AMReX_MFIter.H> #include <AMReX_PODVector.H> @@ -405,19 +406,19 @@ void PsatdAlgorithmJConstantInTime::InitializeSpectralCoefficients ( { // Calculate norm of k vector const amrex::Real knorm_s = std::sqrt( - std::pow(kx_s[i], 2) + + amrex::Math::powi<2>(kx_s[i]) + #if defined(WARPX_DIM_3D) - std::pow(ky_s[j], 2) + std::pow(kz_s[k], 2)); + amrex::Math::powi<2>(ky_s[j]) + amrex::Math::powi<2>(kz_s[k])); #else - std::pow(kz_s[j], 2)); + amrex::Math::powi<2>(kz_s[j])); #endif // Physical constants and imaginary unit constexpr amrex::Real c = PhysConst::c; constexpr amrex::Real ep0 = PhysConst::ep0; constexpr Complex I = Complex{0._rt, 1._rt}; - const amrex::Real c2 = std::pow(c, 2); - const amrex::Real dt2 = std::pow(dt, 2); + const amrex::Real c2 = amrex::Math::powi<2>(c); + const amrex::Real dt2 = amrex::Math::powi<2>(dt); // Calculate the dot product of the k vector with the Galilean velocity. // This has to be computed always with the centered (collocated) finite-order @@ -429,10 +430,10 @@ void PsatdAlgorithmJConstantInTime::InitializeSpectralCoefficients ( #else kz_c[j]*vg_z; #endif - const amrex::Real w2_c = std::pow(w_c, 2); + const amrex::Real w2_c = amrex::Math::powi<2>(w_c); const amrex::Real om_s = c * knorm_s; - const amrex::Real om2_s = std::pow(om_s, 2); + const amrex::Real om2_s = amrex::Math::powi<2>(om_s); const Complex theta_c = amrex::exp( I * w_c * dt * 0.5_rt); const Complex theta2_c = amrex::exp( I * w_c * dt); @@ -566,19 +567,19 @@ void PsatdAlgorithmJConstantInTime::InitializeSpectralCoefficientsAveraging ( { // Calculate norm of k vector const amrex::Real knorm_s = std::sqrt( - std::pow(kx_s[i], 2) + + amrex::Math::powi<2>(kx_s[i]) + #if defined(WARPX_DIM_3D) - std::pow(ky_s[j], 2) + std::pow(kz_s[k], 2)); + amrex::Math::powi<2>(ky_s[j]) + amrex::Math::powi<2>(kz_s[k])); #else - std::pow(kz_s[j], 2)); + amrex::Math::powi<2>(kz_s[j])); #endif // Physical constants and imaginary unit constexpr amrex::Real c = PhysConst::c; constexpr amrex::Real ep0 = PhysConst::ep0; constexpr Complex I = Complex{0._rt, 1._rt}; - const amrex::Real c2 = std::pow(c, 2); - const amrex::Real dt2 = std::pow(dt, 2); + const amrex::Real c2 = amrex::Math::powi<2>(c); + const amrex::Real dt2 = amrex::Math::powi<2>(dt); // Calculate the dot product of the k vector with the Galilean velocity. // This has to be computed always with the centered (collocated) finite-order @@ -590,12 +591,12 @@ void PsatdAlgorithmJConstantInTime::InitializeSpectralCoefficientsAveraging ( #else kz_c[j]*vg_z; #endif - const amrex::Real w2_c = std::pow(w_c, 2); - const amrex::Real w3_c = std::pow(w_c, 3); + const amrex::Real w2_c = amrex::Math::powi<2>(w_c); + const amrex::Real w3_c = amrex::Math::powi<3>(w_c); const amrex::Real om_s = c * knorm_s; - const amrex::Real om2_s = std::pow(om_s, 2); - const amrex::Real om4_s = std::pow(om_s, 4); + const amrex::Real om2_s = amrex::Math::powi<2>(om_s); + const amrex::Real om4_s = amrex::Math::powi<4>(om_s); const Complex theta_c = amrex::exp(I * w_c * dt * 0.5_rt); const Complex theta2_c = amrex::exp(I * w_c * dt); diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmJLinearInTime.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmJLinearInTime.cpp index b0580422b..ecdd6169c 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmJLinearInTime.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmJLinearInTime.cpp @@ -18,6 +18,7 @@ #include <AMReX_GpuLaunch.H> #include <AMReX_GpuQualifiers.H> #include <AMReX_IntVect.H> +#include <AMReX_Math.H> #include <AMReX_MFIter.H> #include <AMReX_PODVector.H> @@ -291,21 +292,21 @@ void PsatdAlgorithmJLinearInTime::InitializeSpectralCoefficients ( { // Calculate norm of k vector const amrex::Real knorm_s = std::sqrt( - std::pow(kx_s[i], 2) + + amrex::Math::powi<2>(kx_s[i]) + #if defined(WARPX_DIM_3D) - std::pow(ky_s[j], 2) + std::pow(kz_s[k], 2)); + amrex::Math::powi<2>(ky_s[j]) + amrex::Math::powi<2>(kz_s[k])); #else - std::pow(kz_s[j], 2)); + amrex::Math::powi<2>(kz_s[j])); #endif // Physical constants and imaginary unit constexpr amrex::Real c = PhysConst::c; constexpr amrex::Real ep0 = PhysConst::ep0; - const amrex::Real c2 = std::pow(c, 2); - const amrex::Real dt2 = std::pow(dt, 2); + const amrex::Real c2 = amrex::Math::powi<2>(c); + const amrex::Real dt2 = amrex::Math::powi<2>(dt); const amrex::Real om_s = c * knorm_s; - const amrex::Real om2_s = std::pow(om_s, 2); + const amrex::Real om2_s = amrex::Math::powi<2>(om_s); // C C(i,j,k) = std::cos(om_s * dt); @@ -383,11 +384,11 @@ void PsatdAlgorithmJLinearInTime::InitializeSpectralCoefficientsAveraging ( { // Calculate norm of k vector const amrex::Real knorm_s = std::sqrt( - std::pow(kx_s[i], 2) + + amrex::Math::powi<2>(kx_s[i]) + #if defined(WARPX_DIM_3D) - std::pow(ky_s[j], 2) + std::pow(kz_s[k], 2)); + amrex::Math::powi<2>(ky_s[j]) + amrex::Math::powi<2>(kz_s[k])); #else - std::pow(kz_s[j], 2)); + amrex::Math::powi<2>(kz_s[j])); #endif // Physical constants and imaginary unit constexpr amrex::Real c = PhysConst::c; 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 |