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

* Starting with the pusher, consistently use ParticleReal

* Update benchmarks for background_mcc_dp_psp

7 files changed, 98 insertions, 98 deletions

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)