Use amrex::Math::powi<N>(x) in WarpX (#3846)

* use_powi

* fix bug

* fix bug

* fix bug

* fix bug

3 files changed, 36 insertions, 33 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;