[WIP] Replace sin(theta), cos(theta) with amrex::Math::sincos(theta) when appropriate

Source/Diagnostics/ReducedDiags/BeamRelevant.cpp

@@ -13,6 +13,7 @@
 #include "WarpX.H"
 
 #include <AMReX_GpuQualifiers.H>
+#include <AMReX_Math.H>
 #include <AMReX_PODVector.H>
 #include <AMReX_ParallelDescriptor.H>
 #include <AMReX_ParmParse.H>
@@ -220,8 +221,9 @@ void BeamRelevant::ComputeDiags (int step)
                 const ParticleReal p_y_mean = p_pos1*p_w;
 #elif defined(WARPX_DIM_RZ)
                 const ParticleReal p_theta = p.rdata(PIdx::theta);
-                const ParticleReal p_x_mean = p_pos0*std::cos(p_theta)*p_w;
-                const ParticleReal p_y_mean = p_pos0*std::sin(p_theta)*p_w;
+                auto const [sin_p_theta, cos_p_theta] = amrex::Math::sincos(p_theta);
+                const ParticleReal p_x_mean = p_pos0*cos_p_theta*p_w;
+                const ParticleReal p_y_mean = p_pos0*sin_p_theta*p_w;
 #elif defined(WARPX_DIM_XZ)
                 const ParticleReal p_x_mean = p_pos0*p_w;
                 const ParticleReal p_y_mean = 0;

Source/Diagnostics/WarpXOpenPMD.cpp

@@ -31,6 +31,7 @@
 #include <AMReX_FabArray.H>
 #include <AMReX_GpuQualifiers.H>
 #include <AMReX_IntVect.H>
+#include <AMReX_Math.H>
 #include <AMReX_MFIter.H>
 #include <AMReX_MultiFab.H>
 #include <AMReX_PODVector.H>
@@ -768,8 +769,9 @@ WarpXOpenPMDPlot::DumpToFile (ParticleContainer* pc,
                     );
                     for (auto i=0; i<numParticleOnTile; i++) {
                         auto const r = aos[i].pos(0);  // {0: "r", 1: "z"}
-                        x.get()[i] = r * std::cos(theta[i]);
-                        y.get()[i] = r * std::sin(theta[i]);
+                        auto const [sin_theta_i, cos_theta_i] = amrex::Math::sincos(theta[i]);
+                        x.get()[i] = r * cos_theta_i;
+                        y.get()[i] = r * sin_theta_i;
                     }
                     currSpecies["position"]["x"].storeChunk(x, {offset}, {numParticleOnTile64});
                     currSpecies["position"]["y"].storeChunk(y, {offset}, {numParticleOnTile64});

Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmComoving.cpp

@@ -238,13 +238,15 @@ void PsatdAlgorithmComoving::InitializeSpectralCoefficients (const SpectralKSpac
                 const amrex::Real om2_mod = om_mod * om_mod;
                 const amrex::Real om  = c * knorm;
                 const amrex::Real om2 = om * om;
-                const Complex tmp1 = amrex::exp(  I * om_mod * dt);
-                const Complex tmp2 = amrex::exp(- I * om_mod * dt);
-                const Complex tmp1_sqrt = amrex::exp(  I * om_mod * dt * 0.5_rt);
-                const Complex tmp2_sqrt = amrex::exp(- I * om_mod * dt * 0.5_rt);
+                const auto phi = om_mod * dt;
+                const Complex tmp1 = amrex::exp(  I * phi);
+                const Complex tmp2 = amrex::exp(- I * phi);
+                const Complex tmp1_sqrt = amrex::exp(  I * phi * 0.5_rt);
+                const Complex tmp2_sqrt = amrex::exp(- I * phi * 0.5_rt);
 
-                C   (i,j,k) = std::cos(om_mod * dt);
-                S_ck(i,j,k) = std::sin(om_mod * dt) / om_mod;
+                auto const [sin_phi, cos_phi] = amrex::Math::sincos(phi);
+                C   (i,j,k) = cos_phi;
+                S_ck(i,j,k) = sin_phi / om_mod;
 
                 const amrex::Real nu = - kv / om;
                 const Complex theta      = amrex::exp(  I * nu * om * dt * 0.5_rt);
@@ -327,8 +329,10 @@ void PsatdAlgorithmComoving::InitializeSpectralCoefficients (const SpectralKSpac
                 const amrex::Real om_mod  = c * knorm_mod;
                 const amrex::Real om2_mod = om_mod * om_mod;
 
-                C   (i,j,k) = std::cos(om_mod * dt);
-                S_ck(i,j,k) = std::sin(om_mod * dt) / om_mod;
+                const auto phi = om_mod * dt;
+                auto const [sin_phi, cos_phi] = amrex::Math::sincos(phi);
+                C   (i,j,k) = cos_phi;
+                S_ck(i,j,k) = sin_phi / om_mod;
                 T2(i,j,k) = 1._rt;
 
                 // X1 multiplies i*(k \times J) in the update equation for B

Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmFirstOrder.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>
 
@@ -151,8 +152,7 @@ PsatdAlgorithmFirstOrder::pushSpectralFields (SpectralFieldData& f) const
             const amrex::Real inv_knorm2 = 1._rt/knorm2;
             const amrex::Real inv_knorm4 = 1._rt/knorm4;
 
-            const amrex::Real C = std::cos(c*knorm*dt);
-            const amrex::Real S = std::sin(c*knorm*dt);
+            auto const [S, C] = amrex::Math::sincos(c*knorm*dt);
 
             // Update equations
 

Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmGalileanRZ.cpp

@@ -11,6 +11,8 @@
 #include "Utils/WarpXProfilerWrapper.H"
 #include "WarpX.H"
 
+#include <AMReX_Math.H>
+
 #include <cmath>
 
 using namespace amrex::literals;
@@ -221,16 +223,18 @@ void PsatdAlgorithmGalileanRZ::InitializeSpectralCoefficients (SpectralFieldData
             // Calculate coefficients
             if (k_norm != 0._rt){
 
-                C(i,j,k,mode) = std::cos(c*k_norm*dt);
-                S_ck(i,j,k,mode) = std::sin(c*k_norm*dt)/(c*k_norm);
+                auto const phi = c*k_norm*dt;
+                auto const [sin_phi, cos_phi] = amrex::Math::sincos(phi);
+                C(i,j,k,mode) = cos_phi;
+                S_ck(i,j,k,mode) = sin_phi/(c*k_norm);
 
                 // Calculate dot product with galilean velocity
                 amrex::Real const kv = kz*vz;
 
                 amrex::Real const nu = kv/(k_norm*c);
                 Complex const theta = amrex::exp( 0.5_rt*I*kv*dt );
                 Complex const theta_star = amrex::exp( -0.5_rt*I*kv*dt );
-                Complex const e_theta = amrex::exp( I*c*k_norm*dt );
+                Complex const e_theta = amrex::exp( I*phi );
 
                 Theta2(i,j,k,mode) = theta*theta;
 
@@ -267,10 +271,10 @@ void PsatdAlgorithmGalileanRZ::InitializeSpectralCoefficients (SpectralFieldData
                     X4(i,j,k,mode) = -S_ck(i,j,k,mode)/ep0;
 
                 } else if ( nu == 1._rt) {
-                    X1(i,j,k,mode) = (1._rt - e_theta*e_theta + 2._rt*I*c*k_norm*dt) / (4._rt*c*c*ep0*k_norm*k_norm);
-                    X2(i,j,k,mode) = (3._rt - 4._rt*e_theta + e_theta*e_theta + 2._rt*I*c*k_norm*dt) / (4._rt*ep0*k_norm*k_norm*(1._rt - e_theta));
-                    X3(i,j,k,mode) = (3._rt - 2._rt/e_theta - 2._rt*e_theta + e_theta*e_theta - 2._rt*I*c*k_norm*dt) / (4._rt*ep0*(e_theta - 1._rt)*k_norm*k_norm);
-                    X4(i,j,k,mode) = I*(-1._rt + e_theta*e_theta + 2._rt*I*c*k_norm*dt) / (4._rt*ep0*c*k_norm);
+                    X1(i,j,k,mode) = (1._rt - e_theta*e_theta + 2._rt*I*phi) / (4._rt*c*c*ep0*k_norm*k_norm);
+                    X2(i,j,k,mode) = (3._rt - 4._rt*e_theta + e_theta*e_theta + 2._rt*I*phi) / (4._rt*ep0*k_norm*k_norm*(1._rt - e_theta));
+                    X3(i,j,k,mode) = (3._rt - 2._rt/e_theta - 2._rt*e_theta + e_theta*e_theta - 2._rt*I*phi) / (4._rt*ep0*(e_theta - 1._rt)*k_norm*k_norm);
+                    X4(i,j,k,mode) = I*(-1._rt + e_theta*e_theta + 2._rt*I*phi) / (4._rt*ep0*c*k_norm);
                 }
 
             } else { // Handle k_norm = 0, by using the analytical limit

Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmPml.cpp

@@ -21,6 +21,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>
 
@@ -443,8 +444,10 @@ void PsatdAlgorithmPml::InitializeSpectralCoefficients (
             // Coefficients for knorm = 0 do not need to be set
             if (knorm != 0._rt)
             {
-                C(i,j,k) = std::cos(c*knorm*dt);
-                S_ck(i,j,k) = std::sin(c*knorm*dt) / (c*knorm);
+                auto const phi = c*knorm*dt;
+                auto const [sin_phi, cos_phi] = amrex::Math::sincos(phi);
+                C(i,j,k) = cos_phi;
+                S_ck(i,j,k) = sin_phi / (c*knorm);
                 inv_k2(i,j,k) = 1._rt / (knorm*knorm);
 
                 if (is_galilean)

Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmPmlRZ.cpp

@@ -11,6 +11,8 @@
 #include "Utils/WarpXProfilerWrapper.H"
 #include "WarpX.H"
 
+#include <AMReX_Math.H>
+
 #include <cmath>
 
 using namespace amrex::literals;
@@ -146,8 +148,10 @@ void PsatdAlgorithmPmlRZ::InitializeSpectralCoefficients (SpectralFieldDataRZ co
             // Calculate coefficients
             constexpr amrex::Real c = PhysConst::c;
             if (k_norm != 0._rt){
-                C(i,j,k,mode) = std::cos(c*k_norm*dt);
-                S_ck(i,j,k,mode) = std::sin(c*k_norm*dt)/(c*k_norm);
+                auto const phi = c*k_norm*dt;
+                auto const [sin_phi, cos_phi] = amrex::Math::sincos(phi);
+                C(i,j,k,mode) = cos_phi;
+                S_ck(i,j,k,mode) = sin_phi/(c*k_norm);
             } else { // Handle k_norm = 0, by using the analytical limit
                 C(i,j,k,mode) = 1._rt;
                 S_ck(i,j,k,mode) = dt;

Source/Initialization/InjectorMomentum.H

@@ -19,6 +19,7 @@
 #include <AMReX_Config.H>
 #include <AMReX_Dim3.H>
 #include <AMReX_GpuQualifiers.H>
+#include <AMReX_Math.H>
 #include <AMReX_REAL.H>
 #include <AMReX_Parser.H>
 #include <AMReX_Random.H>
@@ -404,8 +405,9 @@ struct InjectorMomentumJuttner
         x2 = amrex::Random(engine);
         // Direction dir is an input parameter that sets the boost direction:
         // 'x' -> d = 0, 'y' -> d = 1, 'z' -> d = 2.
-        u[(dir+1)%3] = 2._rt*u[dir]*std::sqrt(x1*(1._rt-x1))*std::sin(2._rt*MathConst::pi*x2);
-        u[(dir+2)%3] = 2._rt*u[dir]*std::sqrt(x1*(1._rt-x1))*std::cos(2._rt*MathConst::pi*x2);
+        auto const [sin_phi, cos_phi] = amrex::Math::sincospi(2._rt*x2);
+        u[(dir+1)%3] = 2._rt*u[dir]*std::sqrt(x1*(1._rt-x1))*sin_phi;
+        u[(dir+2)%3] = 2._rt*u[dir]*std::sqrt(x1*(1._rt-x1))*cos_phi;
         // The value of dir is the boost direction to be transformed.
         u[dir] = u[dir]*(2._rt*x1-1._rt);
         x1 = amrex::Random(engine);

Source/Laser/LaserProfilesImpl/LaserProfileFromFile.cpp

@@ -22,6 +22,7 @@
 #include <AMReX_GpuDevice.H>
 #include <AMReX_GpuLaunch.H>
 #include <AMReX_GpuQualifiers.H>
+#include <AMReX_Math.H>
 #include <AMReX_PODVector.H>
 #include <AMReX_ParallelDescriptor.H>
 #include <AMReX_ParmParse.H>
@@ -433,8 +434,8 @@ WarpXLaserProfiles::FromFileLaserProfile::internal_fill_amplitude_uniform_cartes
 {
     // Copy member variables to tmp copies
     // and get pointers to underlying data for GPU.
-    const amrex::Real omega_t = 2._rt*MathConst::pi*PhysConst::c*t/m_common_params.wavelength;
-    const Complex exp_omega_t = Complex{ std::cos(-omega_t), std::sin(-omega_t) };
+    const auto [omega_im, omega_real] = amrex::Math::sincospi(-2._rt*PhysConst::c*t/m_common_params.wavelength);
+    const Complex exp_omega_t = Complex{ omega_real, omega_im };
     const auto tmp_x_min = m_params.x_min;
     const auto tmp_x_max = m_params.x_max;
     const auto tmp_y_min = m_params.y_min;
@@ -518,8 +519,8 @@ WarpXLaserProfiles::FromFileLaserProfile::internal_fill_amplitude_uniform_cylind
 {
     // Copy member variables to tmp copies
     // and get pointers to underlying data for GPU.
-    const amrex::Real omega_t = 2._rt*MathConst::pi*PhysConst::c*t/m_common_params.wavelength;
-    const Complex exp_omega_t = Complex{ std::cos(-omega_t), std::sin(-omega_t) };
+    const auto [omega_im, omega_real] = amrex::Math::sincospi(-2._rt*PhysConst::c*t/m_common_params.wavelength);
+    const Complex exp_omega_t = Complex{ omega_real, omega_im };
     const auto tmp_r_min = m_params.r_min;
     const auto tmp_r_max = m_params.r_max;
     const auto tmp_nr = m_params.nr;

Source/Laser/LaserProfilesImpl/LaserProfileGaussian.cpp

@@ -18,6 +18,7 @@
 #include <AMReX_GpuComplex.H>
 #include <AMReX_GpuLaunch.H>
 #include <AMReX_GpuQualifiers.H>
+#include <AMReX_Math.H>
 #include <AMReX_ParmParse.H>
 #include <AMReX_REAL.H>
 #include <AMReX_Vector.H>
@@ -146,10 +147,11 @@ WarpXLaserProfiles::GaussianLaserProfile::fill_amplitude (
     amrex::ParallelFor(
         np,
         [=] AMREX_GPU_DEVICE (int i) {
+            auto const [sin_tmp_theta_stc, cos_tmp_theta_stc] = amrex::Math::sincos(tmp_theta_stc);
             const Complex stc_exponent = 1._rt / stretch_factor * inv_tau2 *
                 amrex::pow((t - tmp_profile_t_peak -
-                    tmp_beta*k0*(Xp[i]*std::cos(tmp_theta_stc) + Yp[i]*std::sin(tmp_theta_stc)) -
-                    2._rt *I*(Xp[i]*std::cos(tmp_theta_stc) + Yp[i]*std::sin(tmp_theta_stc))
+                    tmp_beta*k0*(Xp[i]*cos_tmp_theta_stc + Yp[i]*sin_tmp_theta_stc) -
+                    2._rt *I*(Xp[i]*cos_tmp_theta_stc + Yp[i]*sin_tmp_theta_stc)
                     *( tmp_zeta - tmp_beta*tmp_profile_focal_distance ) * inv_complex_waist_2),2);
             // stcfactor = everything but complex transverse envelope
             const Complex stcfactor = prefactor * amrex::exp( - stc_exponent );

Source/Particles/Collision/BinaryCollision/Coulomb/ElasticCollisionPerez.H

@@ -12,6 +12,7 @@
 #include "Particles/WarpXParticleContainer.H"
 #include "Utils/WarpXConst.H"
 
+#include <AMReX_Math.H>
 #include <AMReX_Random.H>
 
 
@@ -141,8 +142,9 @@ void ElasticCollisionPerez (
            * there is a corresponding assert statement at initialization.) */
           T_PR const theta = theta2[I2[i2]]-theta1[I1[i1]];
           T_PR const u1xbuf = u1x[I1[i1]];
-          u1x[I1[i1]] = u1xbuf*std::cos(theta) - u1y[I1[i1]]*std::sin(theta);
-          u1y[I1[i1]] = u1xbuf*std::sin(theta) + u1y[I1[i1]]*std::cos(theta);
+          const auto [sin_theta, cos_theta] = amrex::Math::sincos(theta);
+          u1x[I1[i1]] = u1xbuf*cos_theta - u1y[I1[i1]]*sin_theta;
+          u1y[I1[i1]] = u1xbuf*sin_theta + u1y[I1[i1]]*cos_theta;
 #endif
 
           UpdateMomentumPerezElastic(
@@ -155,8 +157,8 @@ void ElasticCollisionPerez (
 
 #if (defined WARPX_DIM_RZ)
           T_PR const u1xbuf_new = u1x[I1[i1]];
-          u1x[I1[i1]] = u1xbuf_new*std::cos(-theta) - u1y[I1[i1]]*std::sin(-theta);
-          u1y[I1[i1]] = u1xbuf_new*std::sin(-theta) + u1y[I1[i1]]*std::cos(-theta);
+          u1x[I1[i1]] =  u1xbuf_new*cos_theta + u1y[I1[i1]]*sin_theta;
+          u1y[I1[i1]] = -u1xbuf_new*sin_theta + u1y[I1[i1]]*cos_theta;
 #endif
 
           ++i1; if ( i1 == I1e ) { i1 = I1s; }

Source/Particles/Collision/BinaryCollision/Coulomb/UpdateMomentumPerezElastic.H

@@ -192,9 +192,7 @@ void UpdateMomentumPerezElastic (
         sinXs = std::sqrt(T_PR(1.0) - cosXs*cosXs);
 
         // Get random azimuthal angle
-        T_PR const phis = amrex::Random(engine) * T_PR(2.0) * MathConst::pi;
-        T_PR const cosphis = std::cos(phis);
-        T_PR const sinphis = std::sin(phis);
+        auto const [sinphis, cosphis] = amrex::Math::sincospi(amrex::Random(engine) * T_PR(2.0));
 
         // Compute post-collision momenta pfs in COM
         T_PR p1fsx;

Source/Particles/Collision/BinaryCollision/DSMC/CollisionFilterFunc.H

@@ -12,6 +12,7 @@
 #include "Particles/Collision/BinaryCollision/BinaryCollisionUtils.H"
 #include "Particles/Collision/ScatteringProcess.H"
 
+#include <AMReX_Math.H>
 #include <AMReX_Random.H>
 
 /**
@@ -194,8 +195,9 @@ void CollisionFilter (
             * there is a corresponding assert statement at initialization.) */
         amrex::ParticleReal const theta = theta2[I2[i2]]-theta1[I1[i1]];
         amrex::ParticleReal const u1xbuf = u1x[I1[i1]];
-        u1x[I1[i1]] = u1xbuf*std::cos(theta) - u1y[I1[i1]]*std::sin(theta);
-        u1y[I1[i1]] = u1xbuf*std::sin(theta) + u1y[I1[i1]]*std::cos(theta);
+        const auto [sin_theta, cos_theta] = amrex::Math::sincos(theta);
+        u1x[I1[i1]] = u1xbuf*cos_theta - u1y[I1[i1]]*sin_theta;
+        u1y[I1[i1]] = u1xbuf*sin_theta + u1y[I1[i1]]*cos_theta;
 #endif
 
         CollisionPairFilter(
@@ -207,8 +209,8 @@ void CollisionFilter (
 
 #if (defined WARPX_DIM_RZ)
         amrex::ParticleReal const u1xbuf_new = u1x[I1[i1]];
-        u1x[I1[i1]] = u1xbuf_new*std::cos(-theta) - u1y[I1[i1]]*std::sin(-theta);
-        u1y[I1[i1]] = u1xbuf_new*std::sin(-theta) + u1y[I1[i1]]*std::cos(-theta);
+        u1x[I1[i1]] =  u1xbuf_new*cos_theta + u1y[I1[i1]]*sin_theta;
+        u1y[I1[i1]] = -u1xbuf_new*sin_theta + u1y[I1[i1]]*cos_theta;
 #endif
 
         p_pair_indices_1[pair_index] = I1[i1];

Source/Particles/Collision/BinaryCollision/NuclearFusion/NuclearFusionFunc.H

@@ -20,6 +20,7 @@
 
 #include <AMReX_Algorithm.H>
 #include <AMReX_DenseBins.H>
+#include <AMReX_Math.H>
 #include <AMReX_ParmParse.H>
 #include <AMReX_Random.H>
 #include <AMReX_REAL.H>
@@ -184,8 +185,9 @@ public:
                  * there is a corresponding assert statement at initialization.) */
                 amrex::ParticleReal const theta = theta2[I2[i2]]-theta1[I1[i1]];
                 amrex::ParticleReal const u1xbuf = u1x[I1[i1]];
-                u1x[I1[i1]] = u1xbuf*std::cos(theta) - u1y[I1[i1]]*std::sin(theta);
-                u1y[I1[i1]] = u1xbuf*std::sin(theta) + u1y[I1[i1]]*std::cos(theta);
+                auto const [sin_theta, cos_theta] = amrex::Math::sincos(theta);
+                u1x[I1[i1]] = u1xbuf*cos_theta - u1y[I1[i1]]*sin_theta;
+                u1y[I1[i1]] = u1xbuf*sin_theta + u1y[I1[i1]]*cos_theta;
 #endif
 
                 SingleNuclearFusionEvent(
@@ -200,8 +202,8 @@ public:
 
 #if (defined WARPX_DIM_RZ)
                 amrex::ParticleReal const u1xbuf_new = u1x[I1[i1]];
-                u1x[I1[i1]] = u1xbuf_new*std::cos(-theta) - u1y[I1[i1]]*std::sin(-theta);
-                u1y[I1[i1]] = u1xbuf_new*std::sin(-theta) + u1y[I1[i1]]*std::cos(-theta);
+                u1x[I1[i1]] =  u1xbuf_new*cos_theta + u1y[I1[i1]]*sin_theta;
+                u1y[I1[i1]] = -u1xbuf_new*sin_theta + u1y[I1[i1]]*cos_theta;
 #endif
 
                 p_pair_indices_1[pair_index] = I1[i1];

Source/Particles/LaserParticleContainer.cpp

@@ -522,10 +522,10 @@ LaserParticleContainer::InitData (int lev)
                     const auto n_spokes =
                         static_cast<int>((WarpX::n_rz_azimuthal_modes - 1)*m_min_particles_per_mode);
                     for (int spoke = 0 ; spoke < n_spokes ; spoke++) {
-                        const Real phase = 2._rt*MathConst::pi*spoke/n_spokes;
                         for (int k = 0; k<2; ++k) {
-                            particle_x.push_back(pos[0]*std::cos(phase));
-                            particle_y.push_back(pos[0]*std::sin(phase));
+                            auto const [sin_phase, cos_phase] = amrex::Math::sincospi(2._rt*spoke/n_spokes);
+                            particle_x.push_back(pos[0]*cos_phase);
+                            particle_y.push_back(pos[0]*sin_phase);
                             particle_z.push_back(pos[2]);
                         }
                         const Real r_weight = m_weight*2._rt*MathConst::pi*pos[0]/n_spokes;

Source/Particles/ParticleBoundaries_K.H

@@ -10,6 +10,7 @@
 #include "ParticleBoundaries.H"
 
 #include <AMReX_AmrCore.H>
+#include <AMReX_Math.H>
 
 namespace ApplyParticleBoundaries {
 
@@ -133,11 +134,12 @@ namespace ApplyParticleBoundaries {
         } else if (change_sign_ux) {
             // Reflect only ur
             // Note that y is theta
-            amrex::Real ur = ux*std::cos(y) + uy*std::sin(y);
-            const amrex::Real ut = -ux*std::sin(y) + uy*std::cos(y);
+            auto const [siny, cosy] = amrex::Math::sincos(y);
+            amrex::Real ur = ux*cosy + uy*siny;
+            const amrex::Real ut = -ux*siny + uy*cosy;
             ur = -ur;
-            ux = ur*std::cos(y) - ut*std::sin(y);
-            uy = ur*std::sin(y) + ut*std::cos(y);
+            ux = ur*cosy - ut*siny;
+            uy = ur*siny + ut*cosy;
         }
 #else
         if (change_sign_ux) { ux = -ux; }