Starting with the pusher, consistently use ParticleReal (#3259)

* Starting with the pusher, consistently use ParticleReal

* Update benchmarks for background_mcc_dp_psp

13 files changed, 146 insertions, 146 deletions

diff --git a/Source/Particles/MultiParticleContainer.H b/Source/Particles/MultiParticleContainer.H
index 88981c110..bf68dc068 100644
--- a/Source/Particles/MultiParticleContainer.H
+++ b/Source/Particles/MultiParticleContainer.H
@@ -85,7 +85,7 @@ public:
       return  allcontainers[ispecies];
     }
 #endif
-    std::array<amrex::Real, 3> meanParticleVelocity(int ispecies) {
+    std::array<amrex::ParticleReal, 3> meanParticleVelocity(int ispecies) {
         return allcontainers[ispecies]->meanParticleVelocity();
     }
 
diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp
index 8dfca52ef..c075c6d25 100644
--- a/Source/Particles/PhysicalParticleContainer.cpp
+++ b/Source/Particles/PhysicalParticleContainer.cpp
@@ -2646,8 +2646,8 @@ PhysicalParticleContainer::PushPX (WarpXParIter& pti,
     }
 
     // Loop over the particles and update their momentum
-    const amrex::Real q = this->charge;
-    const amrex::Real m = this-> mass;
+    const amrex::ParticleReal q = this->charge;
+    const amrex::ParticleReal m = this-> mass;
 
     const auto pusher_algo = WarpX::particle_pusher_algo;
     const auto do_crr = do_classical_radiation_reaction;
diff --git a/Source/Particles/Pusher/PushSelector.H b/Source/Particles/Pusher/PushSelector.H
index 3022471d9..ed439b4b3 100644
--- a/Source/Particles/Pusher/PushSelector.H
+++ b/Source/Particles/Pusher/PushSelector.H
@@ -55,8 +55,8 @@ void doParticlePush(const GetParticlePosition& GetPosition,
                     const amrex::ParticleReal By,
                     const amrex::ParticleReal Bz,
                     const int ion_lev,
-                    const amrex::Real m,
-                    const amrex::Real q,
+                    const amrex::ParticleReal m,
+                    const amrex::ParticleReal q,
                     const int pusher_algo,
                     const int do_crr,
                     const int do_copy,
@@ -97,7 +97,7 @@ void doParticlePush(const GetParticlePosition& GetPosition,
             SetPosition(i, x, y, z);
         }
 #else
-        amrex::Real qp = q;
+        amrex::ParticleReal qp = q;
         if (ion_lev) { qp *= ion_lev; }
         UpdateMomentumBorisWithRadiationReaction(ux, uy, uz,
                                                  Ex, Ey, Ez, Bx,
@@ -108,7 +108,7 @@ void doParticlePush(const GetParticlePosition& GetPosition,
         SetPosition(i, x, y, z);
 #endif
     } else if (pusher_algo == ParticlePusherAlgo::Boris) {
-        amrex::Real qp = q;
+        amrex::ParticleReal qp = q;
         if (ion_lev) { qp *= ion_lev; }
         UpdateMomentumBoris( ux, uy, uz,
                              Ex, Ey, Ez, Bx,
@@ -118,7 +118,7 @@ void doParticlePush(const GetParticlePosition& GetPosition,
         UpdatePosition(x, y, z, ux, uy, uz, dt );
         SetPosition(i, x, y, z);
     } else if (pusher_algo == ParticlePusherAlgo::Vay) {
-        amrex::Real qp = q;
+        amrex::ParticleReal qp = q;
         if (ion_lev){ qp *= ion_lev; }
         UpdateMomentumVay( ux, uy, uz,
                            Ex, Ey, Ez, Bx,
@@ -128,7 +128,7 @@ void doParticlePush(const GetParticlePosition& GetPosition,
         UpdatePosition(x, y, z, ux, uy, uz, dt );
         SetPosition(i, x, y, z);
     } else if (pusher_algo == ParticlePusherAlgo::HigueraCary) {
-        amrex::Real qp = q;
+        amrex::ParticleReal qp = q;
         if (ion_lev){ qp *= ion_lev; }
         UpdateMomentumHigueraCary( ux, uy, uz,
                                    Ex, Ey, Ez, Bx,
diff --git a/Source/Particles/Pusher/UpdateMomentumBoris.H b/Source/Particles/Pusher/UpdateMomentumBoris.H
index ecb379956..c26a97aa4 100644
--- a/Source/Particles/Pusher/UpdateMomentumBoris.H
+++ b/Source/Particles/Pusher/UpdateMomentumBoris.H
@@ -17,31 +17,31 @@ void UpdateMomentumBoris(
     amrex::ParticleReal& ux, amrex::ParticleReal& uy, amrex::ParticleReal& uz,
     const amrex::ParticleReal Ex, const amrex::ParticleReal Ey, const amrex::ParticleReal Ez,
     const amrex::ParticleReal Bx, const amrex::ParticleReal By, const amrex::ParticleReal Bz,
-    const amrex::Real q, const amrex::Real m, const amrex::Real dt )
+    const amrex::ParticleReal q, const amrex::ParticleReal m, const amrex::Real dt )
 {
     using namespace amrex::literals;
 
-    const amrex::Real econst = 0.5_rt*q*dt/m;
+    const amrex::ParticleReal econst = 0.5_prt*q*dt/m;
 
     // First half-push for E
     ux += econst*Ex;
     uy += econst*Ey;
     uz += econst*Ez;
     // Compute temporary gamma factor
-    constexpr amrex::Real inv_c2 = 1._rt/(PhysConst::c*PhysConst::c);
-    const amrex::Real inv_gamma = 1._rt/std::sqrt(1._rt + (ux*ux + uy*uy + uz*uz)*inv_c2);
+    constexpr amrex::ParticleReal inv_c2 = 1._prt/(PhysConst::c*PhysConst::c);
+    const amrex::ParticleReal inv_gamma = 1._prt/std::sqrt(1._prt + (ux*ux + uy*uy + uz*uz)*inv_c2);
     // Magnetic rotation
     // - Compute temporary variables
-    const amrex::Real tx = econst*inv_gamma*Bx;
-    const amrex::Real ty = econst*inv_gamma*By;
-    const amrex::Real tz = econst*inv_gamma*Bz;
-    const amrex::Real tsqi = 2._rt/(1._rt + tx*tx + ty*ty + tz*tz);
-    const amrex::Real sx = tx*tsqi;
-    const amrex::Real sy = ty*tsqi;
-    const amrex::Real sz = tz*tsqi;
-    const amrex::Real ux_p = ux + uy*tz - uz*ty;
-    const amrex::Real uy_p = uy + uz*tx - ux*tz;
-    const amrex::Real uz_p = uz + ux*ty - uy*tx;
+    const amrex::ParticleReal tx = econst*inv_gamma*Bx;
+    const amrex::ParticleReal ty = econst*inv_gamma*By;
+    const amrex::ParticleReal tz = econst*inv_gamma*Bz;
+    const amrex::ParticleReal tsqi = 2._prt/(1._prt + tx*tx + ty*ty + tz*tz);
+    const amrex::ParticleReal sx = tx*tsqi;
+    const amrex::ParticleReal sy = ty*tsqi;
+    const amrex::ParticleReal sz = tz*tsqi;
+    const amrex::ParticleReal ux_p = ux + uy*tz - uz*ty;
+    const amrex::ParticleReal uy_p = uy + uz*tx - ux*tz;
+    const amrex::ParticleReal uz_p = uz + ux*ty - uy*tx;
     // - Update momentum
     ux += uy_p*sz - uz_p*sy;
     uy += uz_p*sx - ux_p*sz;
diff --git a/Source/Particles/Pusher/UpdateMomentumBorisWithRadiationReaction.H b/Source/Particles/Pusher/UpdateMomentumBorisWithRadiationReaction.H
index 6f0eae476..b637f865a 100644
--- a/Source/Particles/Pusher/UpdateMomentumBorisWithRadiationReaction.H
+++ b/Source/Particles/Pusher/UpdateMomentumBorisWithRadiationReaction.H
@@ -22,17 +22,17 @@ void UpdateMomentumBorisWithRadiationReaction(
     amrex::ParticleReal& ux, amrex::ParticleReal& uy, amrex::ParticleReal& uz,
     const amrex::ParticleReal Ex, const amrex::ParticleReal Ey, const amrex::ParticleReal Ez,
     const amrex::ParticleReal Bx, const amrex::ParticleReal By, const amrex::ParticleReal Bz,
-    const amrex::Real q, const amrex::Real m, const amrex::Real dt )
+    const amrex::ParticleReal q, const amrex::ParticleReal m, const amrex::Real dt )
 {
     using namespace amrex::literals;
 
     //RR algorithm needs to store old value of the normalized momentum
-    const amrex::Real ux_old = ux;
-    const amrex::Real uy_old = uy;
-    const amrex::Real uz_old = uz;
+    const amrex::ParticleReal ux_old = ux;
+    const amrex::ParticleReal uy_old = uy;
+    const amrex::ParticleReal uz_old = uz;
 
     //Useful constant
-    constexpr amrex::Real inv_c2 = 1._rt/(PhysConst::c*PhysConst::c);
+    constexpr amrex::ParticleReal inv_c2 = 1._prt/(PhysConst::c*PhysConst::c);
 
     //Call to regular Boris pusher
     UpdateMomentumBoris(
@@ -42,44 +42,44 @@ void UpdateMomentumBorisWithRadiationReaction(
         q, m, dt );
 
     //Estimation of the normalized momentum at intermediate (integer) time
-    const amrex::Real ux_n = (ux+ux_old)*0.5_rt;
-    const amrex::Real uy_n = (uy+uy_old)*0.5_rt;
-    const amrex::Real uz_n = (uz+uz_old)*0.5_rt;
+    const amrex::ParticleReal ux_n = (ux+ux_old)*0.5_prt;
+    const amrex::ParticleReal uy_n = (uy+uy_old)*0.5_prt;
+    const amrex::ParticleReal uz_n = (uz+uz_old)*0.5_prt;
 
     // Compute Lorentz factor (and inverse) at intermediate (integer) time
-    const amrex::Real gamma_n = std::sqrt( 1._rt +
+    const amrex::ParticleReal gamma_n = std::sqrt( 1._prt +
         (ux_n*ux_n + uy_n*uy_n + uz_n*uz_n)*inv_c2);
-    const amrex::Real inv_gamma_n = 1.0_rt/gamma_n;
+    const amrex::ParticleReal inv_gamma_n = 1.0_prt/gamma_n;
 
     //Estimation of the velocity at intermediate (integer) time
-    const amrex::Real vx_n = ux_n*inv_gamma_n;
-    const amrex::Real vy_n = uy_n*inv_gamma_n;
-    const amrex::Real vz_n = uz_n*inv_gamma_n;
-    const amrex::Real bx_n = vx_n/PhysConst::c;
-    const amrex::Real by_n = vy_n/PhysConst::c;
-    const amrex::Real bz_n = vz_n/PhysConst::c;
+    const amrex::ParticleReal vx_n = ux_n*inv_gamma_n;
+    const amrex::ParticleReal vy_n = uy_n*inv_gamma_n;
+    const amrex::ParticleReal vz_n = uz_n*inv_gamma_n;
+    const amrex::ParticleReal bx_n = vx_n/PhysConst::c;
+    const amrex::ParticleReal by_n = vy_n/PhysConst::c;
+    const amrex::ParticleReal bz_n = vz_n/PhysConst::c;
 
     //Lorentz force over charge
-    const amrex::Real flx_q = (Ex + vy_n*Bz - vz_n*By);
-    const amrex::Real fly_q = (Ey + vz_n*Bx - vx_n*Bz);
-    const amrex::Real flz_q = (Ez + vx_n*By - vy_n*Bx);
-    const amrex::Real fl_q2 = flx_q*flx_q + fly_q*fly_q + flz_q*flz_q;
+    const amrex::ParticleReal flx_q = (Ex + vy_n*Bz - vz_n*By);
+    const amrex::ParticleReal fly_q = (Ey + vz_n*Bx - vx_n*Bz);
+    const amrex::ParticleReal flz_q = (Ez + vx_n*By - vy_n*Bx);
+    const amrex::ParticleReal fl_q2 = flx_q*flx_q + fly_q*fly_q + flz_q*flz_q;
 
     //Calculation of auxiliary quantities
-    const amrex::Real bdotE = (bx_n*Ex + by_n*Ey + bz_n*Ez);
-    const amrex::Real bdotE2 = bdotE*bdotE;
-    const amrex::Real coeff = gamma_n*gamma_n*(fl_q2-bdotE2);
+    const amrex::ParticleReal bdotE = (bx_n*Ex + by_n*Ey + bz_n*Ez);
+    const amrex::ParticleReal bdotE2 = bdotE*bdotE;
+    const amrex::ParticleReal coeff = gamma_n*gamma_n*(fl_q2-bdotE2);
 
     //Radiation reaction constant
-    const amrex::Real q_over_mc = q/(m*PhysConst::c);
-    const amrex::Real RRcoeff = (2.0_rt/3.0_rt)*PhysConst::r_e*q_over_mc*q_over_mc;
+    const amrex::ParticleReal q_over_mc = q/(m*PhysConst::c);
+    const amrex::ParticleReal RRcoeff = (2.0_prt/3.0_prt)*PhysConst::r_e*q_over_mc*q_over_mc;
 
     //Compute the components of the RR force
-    const amrex::Real frx =
+    const amrex::ParticleReal frx =
         RRcoeff*(PhysConst::c*(fly_q*Bz - flz_q*By) + bdotE*Ex - coeff*bx_n);
-    const amrex::Real fry =
+    const amrex::ParticleReal fry =
         RRcoeff*(PhysConst::c*(flz_q*Bx - flx_q*Bz) + bdotE*Ey - coeff*by_n);
-    const amrex::Real frz =
+    const amrex::ParticleReal frz =
         RRcoeff*(PhysConst::c*(flx_q*By - fly_q*Bx) + bdotE*Ez - coeff*bz_n);
 
     //Update momentum using the RR force
diff --git a/Source/Particles/Pusher/UpdateMomentumHigueraCary.H b/Source/Particles/Pusher/UpdateMomentumHigueraCary.H
index 9d42788a4..359fc34d5 100644
--- a/Source/Particles/Pusher/UpdateMomentumHigueraCary.H
+++ b/Source/Particles/Pusher/UpdateMomentumHigueraCary.H
@@ -23,43 +23,43 @@ void UpdateMomentumHigueraCary(
     amrex::ParticleReal& ux, amrex::ParticleReal& uy, amrex::ParticleReal& uz,
     const amrex::ParticleReal Ex, const amrex::ParticleReal Ey, const amrex::ParticleReal Ez,
     const amrex::ParticleReal Bx, const amrex::ParticleReal By, const amrex::ParticleReal Bz,
-    const amrex::Real q, const amrex::Real m, const amrex::Real dt )
+    const amrex::ParticleReal q, const amrex::ParticleReal m, const amrex::Real dt )
 {
     using namespace amrex::literals;
 
     // Constants
-    const amrex::Real qmt = 0.5_rt*q*dt/m;
-    constexpr amrex::Real invclight = 1._rt/PhysConst::c;
-    constexpr amrex::Real invclightsq = 1._rt/(PhysConst::c*PhysConst::c);
+    const amrex::ParticleReal qmt = 0.5_prt*q*dt/m;
+    constexpr amrex::ParticleReal invclight = 1._prt/PhysConst::c;
+    constexpr amrex::ParticleReal invclightsq = 1._prt/(PhysConst::c*PhysConst::c);
     // Compute u_minus
-    const amrex::Real umx = ux + qmt*Ex;
-    const amrex::Real umy = uy + qmt*Ey;
-    const amrex::Real umz = uz + qmt*Ez;
+    const amrex::ParticleReal umx = ux + qmt*Ex;
+    const amrex::ParticleReal umy = uy + qmt*Ey;
+    const amrex::ParticleReal umz = uz + qmt*Ez;
     // Compute gamma squared of u_minus
-    amrex::Real gamma = 1._rt + (umx*umx + umy*umy + umz*umz)*invclightsq;
+    amrex::ParticleReal gamma = 1._prt + (umx*umx + umy*umy + umz*umz)*invclightsq;
     // Compute beta and betam squared
-    const amrex::Real betax = qmt*Bx;
-    const amrex::Real betay = qmt*By;
-    const amrex::Real betaz = qmt*Bz;
-    const amrex::Real betam = betax*betax + betay*betay + betaz*betaz;
+    const amrex::ParticleReal betax = qmt*Bx;
+    const amrex::ParticleReal betay = qmt*By;
+    const amrex::ParticleReal betaz = qmt*Bz;
+    const amrex::ParticleReal betam = betax*betax + betay*betay + betaz*betaz;
     // Compute sigma
-    const amrex::Real sigma = gamma - betam;
+    const amrex::ParticleReal sigma = gamma - betam;
     // Get u*
-    const amrex::Real ust = (umx*betax + umy*betay + umz*betaz)*invclight;
+    const amrex::ParticleReal ust = (umx*betax + umy*betay + umz*betaz)*invclight;
     // Get new gamma inversed
-    gamma = 1._rt/std::sqrt(0.5_rt*(sigma + std::sqrt(sigma*sigma + 4._rt*(betam + ust*ust)) ));
+    gamma = 1._prt/std::sqrt(0.5_prt*(sigma + std::sqrt(sigma*sigma + 4._prt*(betam + ust*ust)) ));
     // Compute t
-    const amrex::Real tx = gamma*betax;
-    const amrex::Real ty = gamma*betay;
-    const amrex::Real tz = gamma*betaz;
+    const amrex::ParticleReal tx = gamma*betax;
+    const amrex::ParticleReal ty = gamma*betay;
+    const amrex::ParticleReal tz = gamma*betaz;
     // Compute s
-    const amrex::Real s = 1._rt/(1._rt+(tx*tx + ty*ty + tz*tz));
+    const amrex::ParticleReal s = 1._prt/(1._prt+(tx*tx + ty*ty + tz*tz));
     // Compute um dot t
-    const amrex::Real umt = umx*tx + umy*ty + umz*tz;
+    const amrex::ParticleReal umt = umx*tx + umy*ty + umz*tz;
     // Compute u_plus
-    const amrex::Real upx = s*( umx + umt*tx + umy*tz - umz*ty );
-    const amrex::Real upy = s*( umy + umt*ty + umz*tx - umx*tz );
-    const amrex::Real upz = s*( umz + umt*tz + umx*ty - umy*tx );
+    const amrex::ParticleReal upx = s*( umx + umt*tx + umy*tz - umz*ty );
+    const amrex::ParticleReal upy = s*( umy + umt*ty + umz*tx - umx*tz );
+    const amrex::ParticleReal upz = s*( umz + umt*tz + umx*ty - umy*tx );
     // Get new u
     ux = upx + qmt*Ex + upy*tz - upz*ty;
     uy = upy + qmt*Ey + upz*tx - upx*tz;
diff --git a/Source/Particles/Pusher/UpdateMomentumVay.H b/Source/Particles/Pusher/UpdateMomentumVay.H
index 610d644d2..c52e6505b 100644
--- a/Source/Particles/Pusher/UpdateMomentumVay.H
+++ b/Source/Particles/Pusher/UpdateMomentumVay.H
@@ -21,40 +21,40 @@ void UpdateMomentumVay(
     amrex::ParticleReal& ux, amrex::ParticleReal& uy, amrex::ParticleReal& uz,
     const amrex::ParticleReal Ex, const amrex::ParticleReal Ey, const amrex::ParticleReal Ez,
     const amrex::ParticleReal Bx, const amrex::ParticleReal By, const amrex::ParticleReal Bz,
-    const amrex::Real q, const amrex::Real m, const amrex::Real dt )
+    const amrex::ParticleReal q, const amrex::ParticleReal m, const amrex::Real dt )
 {
     using namespace amrex::literals;
 
     // Constants
-    const amrex::Real econst = q*dt/m;
-    const amrex::Real bconst = 0.5_rt*q*dt/m;
-    constexpr amrex::Real invclight = 1._rt/PhysConst::c;
-    constexpr amrex::Real invclightsq = 1._rt/(PhysConst::c*PhysConst::c);
+    const amrex::ParticleReal econst = q*dt/m;
+    const amrex::ParticleReal bconst = 0.5_prt*q*dt/m;
+    constexpr amrex::ParticleReal invclight = 1._prt/PhysConst::c;
+    constexpr amrex::ParticleReal invclightsq = 1._prt/(PhysConst::c*PhysConst::c);
     // Compute initial gamma
-    const amrex::Real inv_gamma = 1._rt/std::sqrt(1._rt + (ux*ux + uy*uy + uz*uz)*invclightsq);
+    const amrex::ParticleReal inv_gamma = 1._prt/std::sqrt(1._prt + (ux*ux + uy*uy + uz*uz)*invclightsq);
     // Get tau
-    const amrex::Real taux = bconst*Bx;
-    const amrex::Real tauy = bconst*By;
-    const amrex::Real tauz = bconst*Bz;
-    const amrex::Real tausq = taux*taux+tauy*tauy+tauz*tauz;
+    const amrex::ParticleReal taux = bconst*Bx;
+    const amrex::ParticleReal tauy = bconst*By;
+    const amrex::ParticleReal tauz = bconst*Bz;
+    const amrex::ParticleReal tausq = taux*taux+tauy*tauy+tauz*tauz;
     // Get U', gamma'^2
-    const amrex::Real uxpr = ux + econst*Ex + (uy*tauz-uz*tauy)*inv_gamma;
-    const amrex::Real uypr = uy + econst*Ey + (uz*taux-ux*tauz)*inv_gamma;
-    const amrex::Real uzpr = uz + econst*Ez + (ux*tauy-uy*taux)*inv_gamma;
-    const amrex::Real gprsq = (1._rt + (uxpr*uxpr + uypr*uypr + uzpr*uzpr)*invclightsq);
+    const amrex::ParticleReal uxpr = ux + econst*Ex + (uy*tauz-uz*tauy)*inv_gamma;
+    const amrex::ParticleReal uypr = uy + econst*Ey + (uz*taux-ux*tauz)*inv_gamma;
+    const amrex::ParticleReal uzpr = uz + econst*Ez + (ux*tauy-uy*taux)*inv_gamma;
+    const amrex::ParticleReal gprsq = (1._prt + (uxpr*uxpr + uypr*uypr + uzpr*uzpr)*invclightsq);
     // Get u*
-    const amrex::Real ust = (uxpr*taux + uypr*tauy + uzpr*tauz)*invclight;
+    const amrex::ParticleReal ust = (uxpr*taux + uypr*tauy + uzpr*tauz)*invclight;
     // Get new gamma
-    const amrex::Real sigma = gprsq-tausq;
-    const amrex::Real gisq = 2._rt/(sigma + std::sqrt(sigma*sigma + 4._rt*(tausq + ust*ust)) );
+    const amrex::ParticleReal sigma = gprsq-tausq;
+    const amrex::ParticleReal gisq = 2._prt/(sigma + std::sqrt(sigma*sigma + 4._prt*(tausq + ust*ust)) );
     // Get t, s
-    const amrex::Real bg = bconst*std::sqrt(gisq);
-    const amrex::Real tx = bg*Bx;
-    const amrex::Real ty = bg*By;
-    const amrex::Real tz = bg*Bz;
-    const amrex::Real s = 1._rt/(1._rt+tausq*gisq);
+    const amrex::ParticleReal bg = bconst*std::sqrt(gisq);
+    const amrex::ParticleReal tx = bg*Bx;
+    const amrex::ParticleReal ty = bg*By;
+    const amrex::ParticleReal tz = bg*Bz;
+    const amrex::ParticleReal s = 1._prt/(1._prt+tausq*gisq);
     // Get t.u'
-    const amrex::Real tu = tx*uxpr + ty*uypr + tz*uzpr;
+    const amrex::ParticleReal tu = tx*uxpr + ty*uypr + tz*uzpr;
     // Get new U
     ux = s*(uxpr+tx*tu+uypr*tz-uzpr*ty);
     uy = s*(uypr+ty*tu+uzpr*tx-uxpr*tz);
diff --git a/Source/Particles/Pusher/UpdatePosition.H b/Source/Particles/Pusher/UpdatePosition.H
index d1473eb38..4bce5c102 100644
--- a/Source/Particles/Pusher/UpdatePosition.H
+++ b/Source/Particles/Pusher/UpdatePosition.H
@@ -24,10 +24,10 @@ void UpdatePosition(amrex::ParticleReal& x, amrex::ParticleReal& y, amrex::Parti
 {
     using namespace amrex::literals;
 
-    constexpr amrex::Real inv_c2 = 1._rt/(PhysConst::c*PhysConst::c);
+    constexpr amrex::ParticleReal inv_c2 = 1._prt/(PhysConst::c*PhysConst::c);
 
     // Compute inverse Lorentz factor
-    const amrex::Real inv_gamma = 1._rt/std::sqrt(1._rt + (ux*ux + uy*uy + uz*uz)*inv_c2);
+    const amrex::ParticleReal inv_gamma = 1._prt/std::sqrt(1._prt + (ux*ux + uy*uy + uz*uz)*inv_c2);
     // Update positions over one time step
 #if (AMREX_SPACEDIM >= 2)
     x += ux * inv_gamma * dt;
diff --git a/Source/Particles/Pusher/UpdatePositionPhoton.H b/Source/Particles/Pusher/UpdatePositionPhoton.H
index 5e958c2c1..c570968d4 100644
--- a/Source/Particles/Pusher/UpdatePositionPhoton.H
+++ b/Source/Particles/Pusher/UpdatePositionPhoton.H
@@ -28,8 +28,8 @@ void UpdatePositionPhoton(
 
     // Compute speed of light over inverse of momentum modulus, avoiding a division by zero in the
     // case where the photon has exactly zero momentum
-    const amrex::Real u_norm = std::sqrt(ux*ux + uy*uy + uz*uz);
-    const amrex::Real c_over_umod = (u_norm == 0._rt) ? 0._rt: PhysConst::c/u_norm;
+    const amrex::ParticleReal u_norm = std::sqrt(ux*ux + uy*uy + uz*uz);
+    const amrex::ParticleReal c_over_umod = (u_norm == 0._prt) ? 0._prt: PhysConst::c/u_norm;
 
     // Update positions over one time step
 #if (AMREX_SPACEDIM >= 2)
diff --git a/Source/Particles/RigidInjectedParticleContainer.H b/Source/Particles/RigidInjectedParticleContainer.H
index 81b9448d6..949a6f96f 100644
--- a/Source/Particles/RigidInjectedParticleContainer.H
+++ b/Source/Particles/RigidInjectedParticleContainer.H
@@ -109,16 +109,16 @@ public:
 private:
 
     // User input quantities
-    amrex::Real zinject_plane = 0.;
+    amrex::ParticleReal zinject_plane = 0.;
     bool rigid_advance = true; // When true, particles are advance with vzbar before injection
 
-    amrex::Real vzbeam_ave_boosted;
+    amrex::ParticleReal vzbeam_ave_boosted;
 
-    amrex::Vector<amrex::Real> zinject_plane_levels;
+    amrex::Vector<amrex::ParticleReal> zinject_plane_levels;
 
     // Temporary quantites
-    amrex::Real zinject_plane_lev;
-    amrex::Real zinject_plane_lev_previous;
+    amrex::ParticleReal zinject_plane_lev;
+    amrex::ParticleReal zinject_plane_lev_previous;
     bool done_injecting_lev;
 
 };
diff --git a/Source/Particles/RigidInjectedParticleContainer.cpp b/Source/Particles/RigidInjectedParticleContainer.cpp
index 9ca55067e..489f286b5 100644
--- a/Source/Particles/RigidInjectedParticleContainer.cpp
+++ b/Source/Particles/RigidInjectedParticleContainer.cpp
@@ -92,10 +92,10 @@ RigidInjectedParticleContainer::RemapParticles()
 
         // For now, start with the assumption that this will only happen
         // at the start of the simulation.
-        const Real t_lab = 0._rt;
+        const ParticleReal t_lab = 0._prt;
 
-        const Real uz_boost = WarpX::gamma_boost*WarpX::beta_boost*PhysConst::c;
-        const Real csqi = 1._rt/(PhysConst::c*PhysConst::c);
+        const ParticleReal uz_boost = WarpX::gamma_boost*WarpX::beta_boost*PhysConst::c;
+        const ParticleReal csqi = 1._prt/(PhysConst::c*PhysConst::c);
 
         vzbeam_ave_boosted = meanParticleVelocity(false)[2];
 
@@ -121,21 +121,21 @@ RigidInjectedParticleContainer::RemapParticles()
 
                     // Loop over particles
                     const long np = pti.numParticles();
-                    const Real lvzbeam_ave_boosted = vzbeam_ave_boosted;
-                    const Real gamma_boost = WarpX::gamma_boost;
+                    const ParticleReal lvzbeam_ave_boosted = vzbeam_ave_boosted;
+                    const ParticleReal gamma_boost = WarpX::gamma_boost;
                     amrex::ParallelFor( np, [=] AMREX_GPU_DEVICE (long i)
                     {
                         ParticleReal xp, yp, zp;
                         GetPosition(i, xp, yp, zp);
 
-                        const Real gammapr = std::sqrt(1._rt + (uxp[i]*uxp[i] + uyp[i]*uyp[i] + uzp[i]*uzp[i])*csqi);
-                        const Real vzpr = uzp[i]/gammapr;
+                        const ParticleReal gammapr = std::sqrt(1._prt + (uxp[i]*uxp[i] + uyp[i]*uyp[i] + uzp[i]*uzp[i])*csqi);
+                        const ParticleReal vzpr = uzp[i]/gammapr;
 
                         // Back out the value of z_lab
-                        const Real z_lab = (zp + uz_boost*t_lab + gamma_boost*t_lab*vzpr)/(gamma_boost + uz_boost*vzpr*csqi);
+                        const ParticleReal z_lab = (zp + uz_boost*t_lab + gamma_boost*t_lab*vzpr)/(gamma_boost + uz_boost*vzpr*csqi);
 
                         // Time of the particle in the boosted frame given its position in the lab frame at t=0.
-                        const Real tpr = gamma_boost*t_lab - uz_boost*z_lab*csqi;
+                        const ParticleReal tpr = gamma_boost*t_lab - uz_boost*z_lab*csqi;
 
                         // Adjust the position, taking away its motion from its own velocity and adding
                         // the motion from the average velocity
@@ -253,10 +253,10 @@ RigidInjectedParticleContainer::PushPX (WarpXParIter& pti,
 
         // Undo the push for particles not injected yet.
         // The zp are advanced a fixed amount.
-        const amrex::Real z_plane_lev = zinject_plane_lev;
-        const amrex::Real vz_ave_boosted = vzbeam_ave_boosted;
+        const amrex::ParticleReal z_plane_lev = zinject_plane_lev;
+        const amrex::ParticleReal vz_ave_boosted = vzbeam_ave_boosted;
         const bool rigid = rigid_advance;
-        constexpr amrex::Real inv_csq = 1._rt/(PhysConst::c*PhysConst::c);
+        constexpr amrex::ParticleReal inv_csq = 1._prt/(PhysConst::c*PhysConst::c);
         amrex::ParallelFor( np_to_push,
                             [=] AMREX_GPU_DEVICE (long i) {
                                 amrex::ParticleReal xp, yp, zp;
@@ -271,7 +271,7 @@ RigidInjectedParticleContainer::PushPX (WarpXParIter& pti,
                                         zp = z_save[i] + dt*vz_ave_boosted;
                                     }
                                     else {
-                                        const amrex::Real gi = 1._rt/std::sqrt(1._rt + (ux[i]*ux[i]
+                                        const amrex::ParticleReal gi = 1._prt/std::sqrt(1._prt + (ux[i]*ux[i]
                                                              + uy[i]*uy[i] + uz[i]*uz[i])*inv_csq);
                                         zp = z_save[i] + dt*uz[i]*gi;
                                     }
@@ -400,8 +400,8 @@ RigidInjectedParticleContainer::PushP (int lev, Real dt,
             ParticleReal* const AMREX_RESTRICT uz_save = uzp_save.dataPtr();
 
             // Loop over the particles and update their momentum
-            const amrex::Real q = this->charge;
-            const amrex::Real m = this-> mass;
+            const amrex::ParticleReal q = this->charge;
+            const amrex::ParticleReal m = this-> mass;
 
             const auto pusher_algo = WarpX::particle_pusher_algo;
             const auto do_crr = do_classical_radiation_reaction;
@@ -415,8 +415,8 @@ RigidInjectedParticleContainer::PushP (int lev, Real dt,
                 amrex::ParticleReal xp, yp, zp;
                 getPosition(ip, xp, yp, zp);
 
-                amrex::ParticleReal Exp = 0._rt, Eyp = 0._rt, Ezp = 0._rt;
-                amrex::ParticleReal Bxp = 0._rt, Byp = 0._rt, Bzp = 0._rt;
+                amrex::ParticleReal Exp = 0._prt, Eyp = 0._prt, Ezp = 0._prt;
+                amrex::ParticleReal Bxp = 0._prt, Byp = 0._prt, Bzp = 0._prt;
 
                 // first gather E and B to the particle positions
                 doGatherShapeN(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
@@ -426,7 +426,7 @@ RigidInjectedParticleContainer::PushP (int lev, Real dt,
                                nox, galerkin_interpolation);
                 getExternalEB(ip, Exp, Eyp, Ezp, Bxp, Byp, Bzp);
 
-                amrex::Real qp = q;
+                amrex::ParticleReal qp = q;
                 if (ion_lev) { qp *= ion_lev[ip]; }
 
                 if (do_crr) {
diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H
index e735927e9..f194255bf 100644
--- a/Source/Particles/WarpXParticleContainer.H
+++ b/Source/Particles/WarpXParticleContainer.H
@@ -253,11 +253,11 @@ public:
     /// This returns the total charge for all the particles in this ParticleContainer.
     /// This is needed when solving Poisson's equation with periodic boundary conditions.
     ///
-    amrex::Real sumParticleCharge(bool local = false);
+    amrex::ParticleReal sumParticleCharge(bool local = false);
 
-    std::array<amrex::Real, 3> meanParticleVelocity(bool local = false);
+    std::array<amrex::ParticleReal, 3> meanParticleVelocity(bool local = false);
 
-    amrex::Real maxParticleVelocity(bool local = false);
+    amrex::ParticleReal maxParticleVelocity(bool local = false);
 
     void AddNParticles (int lev,
                         int n, const amrex::ParticleReal* x, const amrex::ParticleReal* y, const amrex::ParticleReal* z,
@@ -349,8 +349,8 @@ public:
 protected:
     int species_id;
 
-    amrex::Real charge;
-    amrex::Real mass;
+    amrex::ParticleReal charge;
+    amrex::ParticleReal mass;
     PhysicalSpecies physical_species;
 
     // Controls boundaries for particles exiting the domain
diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp
index 6af4b5aa7..d3060e8e3 100644
--- a/Source/Particles/WarpXParticleContainer.cpp
+++ b/Source/Particles/WarpXParticleContainer.cpp
@@ -343,7 +343,7 @@ WarpXParticleContainer::DepositCurrent (WarpXParIter& pti,
         "Particles shape does not fit within tile (CPU) or guard cells (GPU) used for current deposition");
 
     const std::array<Real,3>& dx = WarpX::CellSize(std::max(depos_lev,0));
-    Real q = this->charge;
+    amrex::ParticleReal q = this->charge;
 
     WARPX_PROFILE_VAR_NS("WarpXParticleContainer::DepositCurrent::CurrentDeposition", blp_deposit);
     WARPX_PROFILE_VAR_NS("WarpXParticleContainer::DepositCurrent::Accumulate", blp_accumulate);
@@ -767,9 +767,9 @@ WarpXParticleContainer::GetChargeDensity (int lev, bool local)
     return rho;
 }
 
-Real WarpXParticleContainer::sumParticleCharge(bool local) {
+amrex::ParticleReal WarpXParticleContainer::sumParticleCharge(bool local) {
 
-    amrex::Real total_charge = 0.0;
+    amrex::ParticleReal total_charge = 0.0;
 
     const int nLevels = finestLevel();
     for (int lev = 0; lev <= nLevels; ++lev)
@@ -792,15 +792,15 @@ Real WarpXParticleContainer::sumParticleCharge(bool local) {
     return total_charge;
 }
 
-std::array<Real, 3> WarpXParticleContainer::meanParticleVelocity(bool local) {
+std::array<ParticleReal, 3> WarpXParticleContainer::meanParticleVelocity(bool local) {
 
-    amrex::Real vx_total = 0.0_rt;
-    amrex::Real vy_total = 0.0_rt;
-    amrex::Real vz_total = 0.0_rt;
+    amrex::ParticleReal vx_total = 0.0_prt;
+    amrex::ParticleReal vy_total = 0.0_prt;
+    amrex::ParticleReal vz_total = 0.0_prt;
 
     amrex::Long np_total = 0;
 
-    amrex::Real inv_clight_sq = 1.0_rt/PhysConst::c/PhysConst::c;
+    amrex::ParticleReal inv_clight_sq = 1.0_prt/PhysConst::c/PhysConst::c;
 
     const int nLevels = finestLevel();
 
@@ -808,7 +808,7 @@ std::array<Real, 3> WarpXParticleContainer::meanParticleVelocity(bool local) {
     if (Gpu::inLaunchRegion())
     {
         ReduceOps<ReduceOpSum, ReduceOpSum, ReduceOpSum> reduce_op;
-        ReduceData<Real, Real, Real> reduce_data(reduce_op);
+        ReduceData<ParticleReal, ParticleReal, ParticleReal> reduce_data(reduce_op);
         using ReduceTuple = typename decltype(reduce_data)::Type;
         for (int lev = 0; lev <= nLevels; ++lev) {
             for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti)
@@ -823,10 +823,10 @@ std::array<Real, 3> WarpXParticleContainer::meanParticleVelocity(bool local) {
                 reduce_op.eval(np, reduce_data,
                                [=] AMREX_GPU_DEVICE (int i) -> ReduceTuple
                                {
-                                   Real usq = (uxp[i]*uxp[i] +
-                                               uyp[i]*uyp[i] +
-                                               uzp[i]*uzp[i])*inv_clight_sq;
-                                   Real gaminv = 1.0_rt/std::sqrt(1.0_rt + usq);
+                                   amrex::ParticleReal usq = (uxp[i]*uxp[i] +
+                                                              uyp[i]*uyp[i] +
+                                                              uzp[i]*uzp[i])*inv_clight_sq;
+                                   amrex::ParticleReal gaminv = 1.0_prt/std::sqrt(1.0_prt + usq);
                                    return {uxp[i]*gaminv,  uyp[i]*gaminv, uzp[i]*gaminv};
                                });
             }
@@ -853,8 +853,8 @@ std::array<Real, 3> WarpXParticleContainer::meanParticleVelocity(bool local) {
                 np_total += pti.numParticles();
 
                 for (unsigned long i = 0; i < ux.size(); i++) {
-                    Real usq = (ux[i]*ux[i] + uy[i]*uy[i] + uz[i]*uz[i])*inv_clight_sq;
-                    Real gaminv = 1.0_rt/std::sqrt(1.0_rt + usq);
+                    amrex::ParticleReal usq = (ux[i]*ux[i] + uy[i]*uy[i] + uz[i]*uz[i])*inv_clight_sq;
+                    amrex::ParticleReal gaminv = 1.0_prt/std::sqrt(1.0_prt + usq);
                     vx_total += ux[i]*gaminv;
                     vy_total += uy[i]*gaminv;
                     vz_total += uz[i]*gaminv;
@@ -870,7 +870,7 @@ std::array<Real, 3> WarpXParticleContainer::meanParticleVelocity(bool local) {
         ParallelDescriptor::ReduceLongSum(np_total);
     }
 
-    std::array<Real, 3> mean_v;
+    std::array<amrex::ParticleReal, 3> mean_v;
     if (np_total > 0) {
         mean_v[0] = vx_total / np_total;
         mean_v[1] = vy_total / np_total;
@@ -880,7 +880,7 @@ std::array<Real, 3> WarpXParticleContainer::meanParticleVelocity(bool local) {
     return mean_v;
 }
 
-Real WarpXParticleContainer::maxParticleVelocity(bool local) {
+amrex::ParticleReal WarpXParticleContainer::maxParticleVelocity(bool local) {
 
     amrex::ParticleReal max_v = 0.0;