diff options
| -rw-r--r-- | Regression/Checksum/benchmarks_json/background_mcc_dp_psp.json | 16 | ||||
| -rw-r--r-- | Source/FieldSolver/ElectrostaticSolver.cpp | 5 | ||||
| -rw-r--r-- | Source/Particles/MultiParticleContainer.H | 2 | ||||
| -rw-r--r-- | Source/Particles/PhysicalParticleContainer.cpp | 4 | ||||
| -rw-r--r-- | Source/Particles/Pusher/PushSelector.H | 12 | ||||
| -rw-r--r-- | Source/Particles/Pusher/UpdateMomentumBoris.H | 28 | ||||
| -rw-r--r-- | Source/Particles/Pusher/UpdateMomentumBorisWithRadiationReaction.H | 56 | ||||
| -rw-r--r-- | Source/Particles/Pusher/UpdateMomentumHigueraCary.H | 46 | ||||
| -rw-r--r-- | Source/Particles/Pusher/UpdateMomentumVay.H | 46 | ||||
| -rw-r--r-- | Source/Particles/Pusher/UpdatePosition.H | 4 | ||||
| -rw-r--r-- | Source/Particles/Pusher/UpdatePositionPhoton.H | 4 | ||||
| -rw-r--r-- | Source/Particles/RigidInjectedParticleContainer.H | 10 | ||||
| -rw-r--r-- | Source/Particles/RigidInjectedParticleContainer.cpp | 36 | ||||
| -rw-r--r-- | Source/Particles/WarpXParticleContainer.H | 10 | ||||
| -rw-r--r-- | Source/Particles/WarpXParticleContainer.cpp | 34 |
15 files changed, 157 insertions, 156 deletions
diff --git a/Regression/Checksum/benchmarks_json/background_mcc_dp_psp.json b/Regression/Checksum/benchmarks_json/background_mcc_dp_psp.json index aec36e4c5..fe97ca66f 100644 --- a/Regression/Checksum/benchmarks_json/background_mcc_dp_psp.json +++ b/Regression/Checksum/benchmarks_json/background_mcc_dp_psp.json @@ -2,11 +2,11 @@ "electrons": { "particle_cpu": 187506.0, "particle_id": 58260371529.0, - "particle_momentum_x": 1.0114377612798924e-18, - "particle_momentum_y": 2.8184898200776714e-19, - "particle_momentum_z": 2.810005779758007e-19, - "particle_position_x": 17134.123168296646, - "particle_position_y": 935.6698541620412, + "particle_momentum_x": 1.011437782317433e-18, + "particle_momentum_y": 2.818489864433302e-19, + "particle_momentum_z": 2.810005819942864e-19, + "particle_position_x": 17134.12344905036, + "particle_position_y": 935.6698553080840, "particle_weight": 61112621534.71875 }, "he_ions": { @@ -20,7 +20,7 @@ "particle_weight": 71973184795.40625 }, "lev=0": { - "rho_electrons": 0.03566096918517184, - "rho_he_ions": 0.04192386049381602 + "rho_electrons": 0.03566096818084508, + "rho_he_ions": 0.04192385953761138 } -}
\ No newline at end of file +} diff --git a/Source/FieldSolver/ElectrostaticSolver.cpp b/Source/FieldSolver/ElectrostaticSolver.cpp index e760864b7..164b24d00 100644 --- a/Source/FieldSolver/ElectrostaticSolver.cpp +++ b/Source/FieldSolver/ElectrostaticSolver.cpp @@ -177,8 +177,9 @@ WarpX::AddSpaceChargeField (WarpXParticleContainer& pc) // Get the particle beta vector bool const local_average = false; // Average across all MPI ranks - std::array<Real, 3> beta = pc.meanParticleVelocity(local_average); - for (Real& beta_comp : beta) beta_comp /= PhysConst::c; // Normalize + std::array<ParticleReal, 3> beta_pr = pc.meanParticleVelocity(local_average); + std::array<Real, 3> beta; + for (int i=0 ; i < static_cast<int>(beta.size()) ; i++) beta[i] = beta_pr[i]/PhysConst::c; // Normalize // Compute the potential phi, by solving the Poisson equation computePhi( rho, phi, beta, pc.self_fields_required_precision, 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; |