diff options
Diffstat (limited to 'Source')
45 files changed, 823 insertions, 671 deletions
diff --git a/Source/BoundaryConditions/PML.cpp b/Source/BoundaryConditions/PML.cpp index edf8c8358..572030f73 100644 --- a/Source/BoundaryConditions/PML.cpp +++ b/Source/BoundaryConditions/PML.cpp @@ -23,7 +23,6 @@ namespace int olo = overlap.smallEnd(idim); int ohi = overlap.bigEnd(idim); int slo = sigma.m_lo; - int shi = sigma.m_hi; int sslo = sigma_star.m_lo; for (int i = olo; i <= ohi+1; ++i) @@ -51,7 +50,6 @@ namespace int olo = overlap.smallEnd(idim); int ohi = overlap.bigEnd(idim); int slo = sigma.m_lo; - int shi = sigma.m_hi; int sslo = sigma_star.m_lo; for (int i = olo; i <= ohi+1; ++i) { diff --git a/Source/Diagnostics/BoostedFrameDiagnostic.cpp b/Source/Diagnostics/BoostedFrameDiagnostic.cpp index abd9e99fe..45343a0cb 100644 --- a/Source/Diagnostics/BoostedFrameDiagnostic.cpp +++ b/Source/Diagnostics/BoostedFrameDiagnostic.cpp @@ -876,37 +876,37 @@ writeParticleData(const WarpXParticleContainer::DiagnosticParticleData& pdata, field_name = name + Concatenate("w_", i_lab, 5) + "_" + std::to_string(MyProc); ofs.open(field_name.c_str(), std::ios::out|std::ios::binary); - writeRealData(pdata.GetRealData(DiagIdx::w).data(), np, ofs); + writeData(pdata.GetRealData(DiagIdx::w).data(), np, ofs); ofs.close(); field_name = name + Concatenate("x_", i_lab, 5) + "_" + std::to_string(MyProc); ofs.open(field_name.c_str(), std::ios::out|std::ios::binary); - writeRealData(pdata.GetRealData(DiagIdx::x).data(), np, ofs); + writeData(pdata.GetRealData(DiagIdx::x).data(), np, ofs); ofs.close(); field_name = name + Concatenate("y_", i_lab, 5) + "_" + std::to_string(MyProc); ofs.open(field_name.c_str(), std::ios::out|std::ios::binary); - writeRealData(pdata.GetRealData(DiagIdx::y).data(), np, ofs); + writeData(pdata.GetRealData(DiagIdx::y).data(), np, ofs); ofs.close(); field_name = name + Concatenate("z_", i_lab, 5) + "_" + std::to_string(MyProc); ofs.open(field_name.c_str(), std::ios::out|std::ios::binary); - writeRealData(pdata.GetRealData(DiagIdx::z).data(), np, ofs); + writeData(pdata.GetRealData(DiagIdx::z).data(), np, ofs); ofs.close(); field_name = name + Concatenate("ux_", i_lab, 5) + "_" + std::to_string(MyProc); ofs.open(field_name.c_str(), std::ios::out|std::ios::binary); - writeRealData(pdata.GetRealData(DiagIdx::ux).data(), np, ofs); + writeData(pdata.GetRealData(DiagIdx::ux).data(), np, ofs); ofs.close(); field_name = name + Concatenate("uy_", i_lab, 5) + "_" + std::to_string(MyProc); ofs.open(field_name.c_str(), std::ios::out|std::ios::binary); - writeRealData(pdata.GetRealData(DiagIdx::uy).data(), np, ofs); + writeData(pdata.GetRealData(DiagIdx::uy).data(), np, ofs); ofs.close(); field_name = name + Concatenate("uz_", i_lab, 5) + "_" + std::to_string(MyProc); ofs.open(field_name.c_str(), std::ios::out|std::ios::binary); - writeRealData(pdata.GetRealData(DiagIdx::uz).data(), np, ofs); + writeData(pdata.GetRealData(DiagIdx::uz).data(), np, ofs); ofs.close(); } diff --git a/Source/Diagnostics/ParticleIO.cpp b/Source/Diagnostics/ParticleIO.cpp index 5cf3f3047..2a9c16aa8 100644 --- a/Source/Diagnostics/ParticleIO.cpp +++ b/Source/Diagnostics/ParticleIO.cpp @@ -174,9 +174,9 @@ PhysicalParticleContainer::ConvertUnits(ConvertDirection convert_direction) { // - momenta are stored as a struct of array, in `attribs` auto& attribs = pti.GetAttribs(); - Real* AMREX_RESTRICT ux = attribs[PIdx::ux].dataPtr(); - Real* AMREX_RESTRICT uy = attribs[PIdx::uy].dataPtr(); - Real* AMREX_RESTRICT uz = attribs[PIdx::uz].dataPtr(); + ParticleReal* AMREX_RESTRICT ux = attribs[PIdx::ux].dataPtr(); + ParticleReal* AMREX_RESTRICT uy = attribs[PIdx::uy].dataPtr(); + ParticleReal* AMREX_RESTRICT uz = attribs[PIdx::uz].dataPtr(); // Loop over the particles and convert momentum const long np = pti.numParticles(); ParallelFor( np, diff --git a/Source/FieldSolver/WarpX_FDTD.H b/Source/FieldSolver/WarpX_FDTD.H index 8ead945e6..b93c0f40a 100644 --- a/Source/FieldSolver/WarpX_FDTD.H +++ b/Source/FieldSolver/WarpX_FDTD.H @@ -16,7 +16,7 @@ void warpx_push_bx_yee(int i, int j, int k, + dtsdz * (Ey(i,j ,k+1) - Ey(i,j,k)); #elif (defined WARPX_DIM_XZ) Bx(i,j,0) += + dtsdz * (Ey(i,j+1,0) - Ey(i,j,0)); -#if (defined WARPX_DIM_RZ) +#elif (defined WARPX_DIM_RZ) if (i != 0 || rmin != 0.) { Bx(i,j,0,0) += + dtsdz * (Ey(i,j+1,0,0) - Ey(i,j,0,0)); } else { @@ -29,9 +29,9 @@ void warpx_push_bx_yee(int i, int j, int k, // (due to the 1/r terms). The following expressions regularize // these divergences by assuming, on axis : // Ez/r = 0/r + dEz/dr - Bx(i,j,0,2*imode-1) = Bx(i,j,0,2*imode-1) - imode*dtsdx*Ez(i+1,j,0,2*imode) & + Bx(i,j,0,2*imode-1) = Bx(i,j,0,2*imode-1) - imode*dtsdx*Ez(i+1,j,0,2*imode) + dtsdz*(Ey(i,j+1,0,2*imode-1) - Ey(i,j,0,2*imode-1)); - Bx(i,j,0,2*imode) = Bx(i,j,0,2*imode) + imode*dtsdx*Ez(i+1,j,0,2*imode-1) & + Bx(i,j,0,2*imode) = Bx(i,j,0,2*imode) + imode*dtsdx*Ez(i+1,j,0,2*imode-1) + dtsdz*(Ey(i,j+1,0,2*imode) - Ey(i,j,0,2*imode)); } else { Bx(i,j,0,2*imode-1) = 0.; @@ -40,14 +40,13 @@ void warpx_push_bx_yee(int i, int j, int k, } else { // Br(i,j,m) = Br(i,j,m) + I*m*dt*Ez(i,j,m)/r + dtsdz*(Et(i,j+1,m) - Et(i,j,m)) const amrex::Real r = rmin*dxinv + i; - Bx(i,j,0,2*imode-1) = Bx(i,j,0,2*imode-1) - imode*dtsdx*Ez(i,j,0,2*imode)/r & + Bx(i,j,0,2*imode-1) = Bx(i,j,0,2*imode-1) - imode*dtsdx*Ez(i,j,0,2*imode)/r + dtsdz*(Ey(i,j+1,0,2*imode-1) - Ey(i,j,0,2*imode-1)); - Bx(i,j,0,2*imode) = Bx(i,j,0,2*imode) + imode*dtsdx*Ez(i,j,0,2*imode-1)/r & + Bx(i,j,0,2*imode) = Bx(i,j,0,2*imode) + imode*dtsdx*Ez(i,j,0,2*imode-1)/r + dtsdz*(Ey(i,j+1,0,2*imode) - Ey(i,j,0,2*imode)); } } #endif -#endif } AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE diff --git a/Source/FortranInterface/WarpX_f.H b/Source/FortranInterface/WarpX_f.H index 9a0aeb819..2d5d346e2 100644 --- a/Source/FortranInterface/WarpX_f.H +++ b/Source/FortranInterface/WarpX_f.H @@ -84,14 +84,14 @@ extern "C" #endif const amrex::Real* dx); - void WRPX_DEPOSIT_CIC(const amrex::Real* particles, int ns, int np, - const amrex::Real* weights, + void WRPX_DEPOSIT_CIC(const amrex::ParticleReal* particles, int ns, int np, + const amrex::ParticleReal* weights, const amrex::Real* charge, amrex::Real* rho, const int* lo, const int* hi, const amrex::Real* plo, const amrex::Real* dx, const int* ng); - void WRPX_INTERPOLATE_CIC_TWO_LEVELS(const amrex::Real* particles, int ns, int np, + void WRPX_INTERPOLATE_CIC_TWO_LEVELS(const amrex::ParticleReal* particles, int ns, int np, amrex::Real* Ex_p, amrex::Real* Ey_p, #if (AMREX_SPACEDIM == 3) amrex::Real* Ez_p, @@ -109,7 +109,7 @@ extern "C" const int* clo, const int* chi, const amrex::Real* cdx, const amrex::Real* plo, const int* ng, const int* lev); - void WRPX_INTERPOLATE_CIC(const amrex::Real* particles, int ns, int np, + void WRPX_INTERPOLATE_CIC(const amrex::ParticleReal* particles, int ns, int np, amrex::Real* Ex_p, amrex::Real* Ey_p, #if (AMREX_SPACEDIM == 3) amrex::Real* Ez_p, @@ -122,10 +122,10 @@ extern "C" const amrex::Real* plo, const amrex::Real* dx, const int* ng); - void WRPX_PUSH_LEAPFROG(amrex::Real* particles, int ns, int np, - amrex::Real* vx_p, amrex::Real* vy_p, + void WRPX_PUSH_LEAPFROG(amrex::ParticleReal* particles, int ns, int np, + amrex::ParticleReal* vx_p, amrex::ParticleReal* vy_p, #if (AMREX_SPACEDIM == 3) - amrex::Real* vz_p, + amrex::ParticleReal* vz_p, #endif const amrex::Real* Ex_p, const amrex::Real* Ey_p, #if (AMREX_SPACEDIM == 3) @@ -134,10 +134,10 @@ extern "C" const amrex::Real* charge, const amrex::Real* mass, const amrex::Real* dt, const amrex::Real* prob_lo, const amrex::Real* prob_hi); - void WRPX_PUSH_LEAPFROG_POSITIONS(amrex::Real* particles, int ns, int np, - amrex::Real* vx_p, amrex::Real* vy_p, + void WRPX_PUSH_LEAPFROG_POSITIONS(amrex::ParticleReal* particles, int ns, int np, + amrex::ParticleReal* vx_p, amrex::ParticleReal* vy_p, #if (AMREX_SPACEDIM == 3) - amrex::Real* vz_p, + amrex::ParticleReal* vz_p, #endif const amrex::Real* dt, const amrex::Real* prob_lo, const amrex::Real* prob_hi); diff --git a/Source/Initialization/InjectorPosition.H b/Source/Initialization/InjectorPosition.H index f8f16746c..6ecae93e0 100644 --- a/Source/Initialization/InjectorPosition.H +++ b/Source/Initialization/InjectorPosition.H @@ -41,7 +41,9 @@ struct InjectorPositionRegular int ix_part = i_part/(ny*nz); // written this way backward compatibility int iz_part = (i_part-ix_part*(ny*nz)) / ny; int iy_part = (i_part-ix_part*(ny*nz)) - ny*iz_part; - return amrex::XDim3{(0.5+ix_part)/nx, (0.5+iy_part)/ny, (0.5+iz_part) / nz}; + return amrex::XDim3{(amrex::Real(0.5)+ix_part)/nx, + (amrex::Real(0.5)+iy_part)/ny, + (amrex::Real(0.5)+iz_part) / nz}; } private: amrex::Dim3 ppc; diff --git a/Source/Initialization/PlasmaInjector.H b/Source/Initialization/PlasmaInjector.H index a944165d6..56b32c827 100644 --- a/Source/Initialization/PlasmaInjector.H +++ b/Source/Initialization/PlasmaInjector.H @@ -43,9 +43,9 @@ public: bool doInjection () const noexcept { return inj_pos != NULL;} bool add_single_particle = false; - amrex::Vector<amrex::Real> single_particle_pos; - amrex::Vector<amrex::Real> single_particle_vel; - amrex::Real single_particle_weight; + amrex::Vector<amrex::ParticleReal> single_particle_pos; + amrex::Vector<amrex::ParticleReal> single_particle_vel; + amrex::ParticleReal single_particle_weight; bool gaussian_beam = false; amrex::Real x_m; diff --git a/Source/Laser/LaserParticleContainer.H b/Source/Laser/LaserParticleContainer.H index 31fe7cee4..e2a0743bc 100644 --- a/Source/Laser/LaserParticleContainer.H +++ b/Source/Laser/LaserParticleContainer.H @@ -56,10 +56,10 @@ public: amrex::Real * AMREX_RESTRICT const pplane_Xp, amrex::Real * AMREX_RESTRICT const pplane_Yp); - void update_laser_particle (const int np, amrex::Real * AMREX_RESTRICT const puxp, - amrex::Real * AMREX_RESTRICT const puyp, - amrex::Real * AMREX_RESTRICT const puzp, - amrex::Real const * AMREX_RESTRICT const pwp, + void update_laser_particle (const int np, amrex::ParticleReal * AMREX_RESTRICT const puxp, + amrex::ParticleReal * AMREX_RESTRICT const puyp, + amrex::ParticleReal * AMREX_RESTRICT const puzp, + amrex::ParticleReal const * AMREX_RESTRICT const pwp, amrex::Real const * AMREX_RESTRICT const amplitude, const amrex::Real dt, const int thread_num); @@ -81,6 +81,8 @@ private: amrex::Real wavelength = std::numeric_limits<amrex::Real>::quiet_NaN(); amrex::Real Z0_lab = 0; // Position of the antenna in the lab frame + long min_particles_per_mode = 4; + // computed using runtime parameters amrex::Vector<amrex::Real> p_Y; amrex::Vector<amrex::Real> u_X; diff --git a/Source/Laser/LaserParticleContainer.cpp b/Source/Laser/LaserParticleContainer.cpp index dca73de50..8571c74ad 100644 --- a/Source/Laser/LaserParticleContainer.cpp +++ b/Source/Laser/LaserParticleContainer.cpp @@ -40,115 +40,116 @@ LaserParticleContainer::LaserParticleContainer (AmrCore* amr_core, int ispecies, profile = laser_t::Harris; } else if(laser_type_s == "parse_field_function") { profile = laser_t::parse_field_function; - } else { - amrex::Abort("Unknown laser type"); - } + } else { + amrex::Abort("Unknown laser type"); + } - // Parse the properties of the antenna - pp.getarr("position", position); - pp.getarr("direction", nvec); - pp.getarr("polarization", p_X); - pp.query("pusher_algo", pusher_algo); - pp.get("wavelength", wavelength); - pp.get("e_max", e_max); - pp.query("do_continuous_injection", do_continuous_injection); - - if ( profile == laser_t::Gaussian ) { - // Parse the properties of the Gaussian profile - pp.get("profile_waist", profile_waist); - pp.get("profile_duration", profile_duration); - pp.get("profile_t_peak", profile_t_peak); - pp.get("profile_focal_distance", profile_focal_distance); - stc_direction = p_X; - pp.queryarr("stc_direction", stc_direction); - pp.query("zeta", zeta); - pp.query("beta", beta); - pp.query("phi2", phi2); - } + // Parse the properties of the antenna + pp.getarr("position", position); + pp.getarr("direction", nvec); + pp.getarr("polarization", p_X); + pp.query("pusher_algo", pusher_algo); + pp.get("wavelength", wavelength); + pp.get("e_max", e_max); + pp.query("do_continuous_injection", do_continuous_injection); + pp.query("min_particles_per_mode", min_particles_per_mode); + + if ( profile == laser_t::Gaussian ) { + // Parse the properties of the Gaussian profile + pp.get("profile_waist", profile_waist); + pp.get("profile_duration", profile_duration); + pp.get("profile_t_peak", profile_t_peak); + pp.get("profile_focal_distance", profile_focal_distance); + stc_direction = p_X; + pp.queryarr("stc_direction", stc_direction); + pp.query("zeta", zeta); + pp.query("beta", beta); + pp.query("phi2", phi2); + } - if ( profile == laser_t::Harris ) { - // Parse the properties of the Harris profile - pp.get("profile_waist", profile_waist); - pp.get("profile_duration", profile_duration); - pp.get("profile_focal_distance", profile_focal_distance); - } + if ( profile == laser_t::Harris ) { + // Parse the properties of the Harris profile + pp.get("profile_waist", profile_waist); + pp.get("profile_duration", profile_duration); + pp.get("profile_focal_distance", profile_focal_distance); + } - if ( profile == laser_t::parse_field_function ) { - // Parse the properties of the parse_field_function profile - pp.get("field_function(X,Y,t)", field_function); - parser.define(field_function); - parser.registerVariables({"X","Y","t"}); - - ParmParse ppc("my_constants"); - std::set<std::string> symbols = parser.symbols(); - symbols.erase("X"); - symbols.erase("Y"); - symbols.erase("t"); // after removing variables, we are left with constants - for (auto it = symbols.begin(); it != symbols.end(); ) { - Real v; - if (ppc.query(it->c_str(), v)) { - parser.setConstant(*it, v); - it = symbols.erase(it); - } else { - ++it; - } - } - for (auto const& s : symbols) { // make sure there no unknown symbols - amrex::Abort("Laser Profile: Unknown symbol "+s); + if ( profile == laser_t::parse_field_function ) { + // Parse the properties of the parse_field_function profile + pp.get("field_function(X,Y,t)", field_function); + parser.define(field_function); + parser.registerVariables({"X","Y","t"}); + + ParmParse ppc("my_constants"); + std::set<std::string> symbols = parser.symbols(); + symbols.erase("X"); + symbols.erase("Y"); + symbols.erase("t"); // after removing variables, we are left with constants + for (auto it = symbols.begin(); it != symbols.end(); ) { + Real v; + if (ppc.query(it->c_str(), v)) { + parser.setConstant(*it, v); + it = symbols.erase(it); + } else { + ++it; } } - - // Plane normal - Real s = 1.0/std::sqrt(nvec[0]*nvec[0] + nvec[1]*nvec[1] + nvec[2]*nvec[2]); - nvec = { nvec[0]*s, nvec[1]*s, nvec[2]*s }; - - if (WarpX::gamma_boost > 1.) { - // Check that the laser direction is equal to the boost direction - AMREX_ALWAYS_ASSERT_WITH_MESSAGE( nvec[0]*WarpX::boost_direction[0] - + nvec[1]*WarpX::boost_direction[1] - + nvec[2]*WarpX::boost_direction[2] - 1. < 1.e-12, - "The Lorentz boost should be in the same direction as the laser propagation"); - // Get the position of the plane, along the boost direction, in the lab frame - // and convert the position of the antenna to the boosted frame - Z0_lab = nvec[0]*position[0] + nvec[1]*position[1] + nvec[2]*position[2]; - Real Z0_boost = Z0_lab/WarpX::gamma_boost; - position[0] += (Z0_boost-Z0_lab)*nvec[0]; - position[1] += (Z0_boost-Z0_lab)*nvec[1]; - position[2] += (Z0_boost-Z0_lab)*nvec[2]; + for (auto const& s : symbols) { // make sure there no unknown symbols + amrex::Abort("Laser Profile: Unknown symbol "+s); } + } - // The first polarization vector - s = 1.0/std::sqrt(p_X[0]*p_X[0] + p_X[1]*p_X[1] + p_X[2]*p_X[2]); - p_X = { p_X[0]*s, p_X[1]*s, p_X[2]*s }; + // Plane normal + Real s = 1.0/std::sqrt(nvec[0]*nvec[0] + nvec[1]*nvec[1] + nvec[2]*nvec[2]); + nvec = { nvec[0]*s, nvec[1]*s, nvec[2]*s }; + + if (WarpX::gamma_boost > 1.) { + // Check that the laser direction is equal to the boost direction + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( nvec[0]*WarpX::boost_direction[0] + + nvec[1]*WarpX::boost_direction[1] + + nvec[2]*WarpX::boost_direction[2] - 1. < 1.e-12, + "The Lorentz boost should be in the same direction as the laser propagation"); + // Get the position of the plane, along the boost direction, in the lab frame + // and convert the position of the antenna to the boosted frame + Z0_lab = nvec[0]*position[0] + nvec[1]*position[1] + nvec[2]*position[2]; + Real Z0_boost = Z0_lab/WarpX::gamma_boost; + position[0] += (Z0_boost-Z0_lab)*nvec[0]; + position[1] += (Z0_boost-Z0_lab)*nvec[1]; + position[2] += (Z0_boost-Z0_lab)*nvec[2]; + } + + // The first polarization vector + s = 1.0/std::sqrt(p_X[0]*p_X[0] + p_X[1]*p_X[1] + p_X[2]*p_X[2]); + p_X = { p_X[0]*s, p_X[1]*s, p_X[2]*s }; - Real dp = std::inner_product(nvec.begin(), nvec.end(), p_X.begin(), 0.0); - AMREX_ALWAYS_ASSERT_WITH_MESSAGE(std::abs(dp) < 1.0e-14, - "Laser plane vector is not perpendicular to the main polarization vector"); + Real dp = std::inner_product(nvec.begin(), nvec.end(), p_X.begin(), 0.0); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(std::abs(dp) < 1.0e-14, + "Laser plane vector is not perpendicular to the main polarization vector"); - p_Y = CrossProduct(nvec, p_X); // The second polarization vector + p_Y = CrossProduct(nvec, p_X); // The second polarization vector - s = 1.0/std::sqrt(stc_direction[0]*stc_direction[0] + stc_direction[1]*stc_direction[1] + stc_direction[2]*stc_direction[2]); - stc_direction = { stc_direction[0]*s, stc_direction[1]*s, stc_direction[2]*s }; - dp = std::inner_product(nvec.begin(), nvec.end(), stc_direction.begin(), 0.0); - AMREX_ALWAYS_ASSERT_WITH_MESSAGE(std::abs(dp) < 1.0e-14, - "stc_direction is not perpendicular to the laser plane vector"); + s = 1.0/std::sqrt(stc_direction[0]*stc_direction[0] + stc_direction[1]*stc_direction[1] + stc_direction[2]*stc_direction[2]); + stc_direction = { stc_direction[0]*s, stc_direction[1]*s, stc_direction[2]*s }; + dp = std::inner_product(nvec.begin(), nvec.end(), stc_direction.begin(), 0.0); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(std::abs(dp) < 1.0e-14, + "stc_direction is not perpendicular to the laser plane vector"); - // Get angle between p_X and stc_direction - // in 2d, stcs are in the simulation plane + // Get angle between p_X and stc_direction + // in 2d, stcs are in the simulation plane #if AMREX_SPACEDIM == 3 - theta_stc = acos(stc_direction[0]*p_X[0] + + theta_stc = acos(stc_direction[0]*p_X[0] + stc_direction[1]*p_X[1] + stc_direction[2]*p_X[2]); #else - theta_stc = 0.; + theta_stc = 0.; #endif -#if AMREX_SPACEDIM == 3 - u_X = p_X; - u_Y = p_Y; +#if (defined WARPX_DIM_3D) || (defined WARPX_DIM_RZ) + u_X = p_X; + u_Y = p_Y; #else - u_X = CrossProduct({0., 1., 0.}, nvec); - u_Y = {0., 1., 0.}; + u_X = CrossProduct({0., 1., 0.}, nvec); + u_Y = {0., 1., 0.}; #endif laser_injection_box= Geom(0).ProbDomain(); @@ -271,9 +272,15 @@ LaserParticleContainer::InitData (int lev) position[1] + (S_X*(i+0.5))*u_X[1] + (S_Y*(j+0.5))*u_Y[1], position[2] + (S_X*(i+0.5))*u_X[2] + (S_Y*(j+0.5))*u_Y[2] }; #else +# if (defined WARPX_DIM_RZ) + return { position[0] + (S_X*(i+0.5)), + 0.0, + position[2]}; +# else return { position[0] + (S_X*(i+0.5))*u_X[0], - 0.0, - position[2] + (S_X*(i+0.5))*u_X[2] }; + 0.0, + position[2] + (S_X*(i+0.5))*u_X[2] }; +# endif #endif }; @@ -283,7 +290,11 @@ LaserParticleContainer::InitData (int lev) return {u_X[0]*(pos[0]-position[0])+u_X[1]*(pos[1]-position[1])+u_X[2]*(pos[2]-position[2]), u_Y[0]*(pos[0]-position[0])+u_Y[1]*(pos[1]-position[1])+u_Y[2]*(pos[2]-position[2])}; #else +# if (defined WARPX_DIM_RZ) + return {pos[0]-position[0], 0.0}; +# else return {u_X[0]*(pos[0]-position[0])+u_X[2]*(pos[2]-position[2]), 0.0}; +# endif #endif }; @@ -364,6 +375,7 @@ LaserParticleContainer::InitData (int lev) #endif if (laser_injection_box.contains(x)) { +#ifndef WARPX_DIM_RZ for (int k = 0; k<2; ++k) { particle_x.push_back(pos[0]); particle_y.push_back(pos[1]); @@ -371,6 +383,21 @@ LaserParticleContainer::InitData (int lev) } particle_w.push_back( weight); particle_w.push_back(-weight); +#else + // Particles are laid out in radial spokes + const int n_spokes = (WarpX::n_rz_azimuthal_modes - 1)*min_particles_per_mode; + for (int spoke = 0 ; spoke < n_spokes ; spoke++) { + const Real phase = 2.*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)); + particle_z.push_back(pos[2]); + } + const Real r_weight = weight*2.*MathConst::pi*pos[0]/n_spokes; + particle_w.push_back( r_weight); + particle_w.push_back(-r_weight); + } +#endif } } } @@ -564,13 +591,17 @@ LaserParticleContainer::ComputeSpacing (int lev, Real& Sx, Real& Sy) const #if (AMREX_SPACEDIM == 3) Sx = std::min(std::min(dx[0]/(std::abs(u_X[0])+eps), dx[1]/(std::abs(u_X[1])+eps)), - dx[2]/(std::abs(u_X[2])+eps)); + dx[2]/(std::abs(u_X[2])+eps)); Sy = std::min(std::min(dx[0]/(std::abs(u_Y[0])+eps), dx[1]/(std::abs(u_Y[1])+eps)), - dx[2]/(std::abs(u_Y[2])+eps)); + dx[2]/(std::abs(u_Y[2])+eps)); #else +# if (defined WARPX_DIM_RZ) + Sx = dx[0]; +# else Sx = std::min(dx[0]/(std::abs(u_X[0])+eps), dx[2]/(std::abs(u_X[2])+eps)); +# endif Sy = 1.0; #endif } @@ -579,7 +610,7 @@ void LaserParticleContainer::ComputeWeightMobility (Real Sx, Real Sy) { constexpr Real eps = 0.01; - constexpr Real fac = 1.0/(2.0*3.1415926535897932*PhysConst::mu0*PhysConst::c*PhysConst::c*eps); + constexpr Real fac = 1.0/(2.0*MathConst::pi*PhysConst::mu0*PhysConst::c*PhysConst::c*eps); weight = fac * wavelength * Sx * Sy / std::min(Sx,Sy) * e_max; // The mobility is the constant of proportionality between the field to @@ -612,14 +643,14 @@ LaserParticleContainer::calculate_laser_plane_coordinates ( Real * AMREX_RESTRICT const pplane_Xp, Real * AMREX_RESTRICT const pplane_Yp) { - Real const * const AMREX_RESTRICT xp = m_xp[thread_num].dataPtr(); - Real const * const AMREX_RESTRICT yp = m_yp[thread_num].dataPtr(); - Real const * const AMREX_RESTRICT zp = m_zp[thread_num].dataPtr(); + ParticleReal const * const AMREX_RESTRICT xp = m_xp[thread_num].dataPtr(); + ParticleReal const * const AMREX_RESTRICT yp = m_yp[thread_num].dataPtr(); + ParticleReal const * const AMREX_RESTRICT zp = m_zp[thread_num].dataPtr(); Real tmp_u_X_0 = u_X[0]; Real tmp_u_X_2 = u_X[2]; Real tmp_position_0 = position[0]; Real tmp_position_2 = position[2]; -#if (AMREX_SPACEDIM == 3) +#if (defined WARPX_DIM_3D) || (defined WARPX_DIM_RZ) Real tmp_u_X_1 = u_X[1]; Real tmp_u_Y_0 = u_Y[0]; Real tmp_u_Y_1 = u_Y[1]; @@ -630,7 +661,7 @@ LaserParticleContainer::calculate_laser_plane_coordinates ( amrex::ParallelFor( np, [=] AMREX_GPU_DEVICE (int i) { -#if (AMREX_SPACEDIM == 3) +#if (defined WARPX_DIM_3D) || (defined WARPX_DIM_RZ) pplane_Xp[i] = tmp_u_X_0 * (xp[i] - tmp_position_0) + tmp_u_X_1 * (yp[i] - tmp_position_1) + @@ -660,14 +691,14 @@ LaserParticleContainer::calculate_laser_plane_coordinates ( */ void LaserParticleContainer::update_laser_particle( - const int np, Real * AMREX_RESTRICT const puxp, Real * AMREX_RESTRICT const puyp, - Real * AMREX_RESTRICT const puzp, Real const * AMREX_RESTRICT const pwp, + const int np, ParticleReal * AMREX_RESTRICT const puxp, ParticleReal * AMREX_RESTRICT const puyp, + ParticleReal * AMREX_RESTRICT const puzp, ParticleReal const * AMREX_RESTRICT const pwp, Real const * AMREX_RESTRICT const amplitude, const Real dt, const int thread_num) { - Real * const AMREX_RESTRICT xp = m_xp[thread_num].dataPtr(); - Real * const AMREX_RESTRICT yp = m_yp[thread_num].dataPtr(); - Real * const AMREX_RESTRICT zp = m_zp[thread_num].dataPtr(); + ParticleReal * const AMREX_RESTRICT xp = m_xp[thread_num].dataPtr(); + ParticleReal * const AMREX_RESTRICT yp = m_yp[thread_num].dataPtr(); + ParticleReal * const AMREX_RESTRICT zp = m_zp[thread_num].dataPtr(); Real tmp_p_X_0 = p_X[0]; Real tmp_p_X_1 = p_X[1]; Real tmp_p_X_2 = p_X[2]; @@ -702,7 +733,7 @@ LaserParticleContainer::update_laser_particle( puzp[i] = gamma * vz; // Push the the particle positions xp[i] += vx * dt; -#if (AMREX_SPACEDIM == 3) +#if (defined WARPX_DIM_3D) || (defined WARPX_DIM_RZ) yp[i] += vy * dt; #endif zp[i] += vz * dt; diff --git a/Source/Laser/LaserProfiles.cpp b/Source/Laser/LaserProfiles.cpp index 69804b17c..281ab2101 100644 --- a/Source/Laser/LaserProfiles.cpp +++ b/Source/Laser/LaserProfiles.cpp @@ -28,16 +28,16 @@ LaserParticleContainer::gaussian_laser_profile ( const Real oscillation_phase = k0 * PhysConst::c * ( t - profile_t_peak ); // The coefficients below contain info about Gouy phase, // laser diffraction, and phase front curvature - const Complex diffract_factor = 1. + I * profile_focal_distance - * 2./( k0 * profile_waist * profile_waist ); - const Complex inv_complex_waist_2 = 1./( profile_waist*profile_waist * diffract_factor ); + const Complex diffract_factor = Real(1.) + I * profile_focal_distance + * Real(2.)/( k0 * profile_waist * profile_waist ); + const Complex inv_complex_waist_2 = Real(1.)/( profile_waist*profile_waist * diffract_factor ); // Time stretching due to STCs and phi2 complex envelope // (1 if zeta=0, beta=0, phi2=0) - const Complex stretch_factor = 1. + 4. * + const Complex stretch_factor = Real(1.) + Real(4.) * (zeta+beta*profile_focal_distance) * (zeta+beta*profile_focal_distance) * (inv_tau2*inv_complex_waist_2) + - 2.*I*(phi2 - beta*beta*k0*profile_focal_distance) * inv_tau2; + Real(2.)*I*(phi2 - beta*beta*k0*profile_focal_distance) * inv_tau2; // Amplitude and monochromatic oscillations Complex prefactor = e_max * MathFunc::exp( I * oscillation_phase ); @@ -61,10 +61,10 @@ LaserParticleContainer::gaussian_laser_profile ( amrex::ParallelFor( np, [=] AMREX_GPU_DEVICE (int i) { - const Complex stc_exponent = 1./stretch_factor * inv_tau2 * + const Complex stc_exponent = Real(1.)/stretch_factor * inv_tau2 * MathFunc::pow((t - tmp_profile_t_peak - tmp_beta*k0*(Xp[i]*std::cos(tmp_theta_stc) + Yp[i]*std::sin(tmp_theta_stc)) - - 2.*I*(Xp[i]*std::cos(tmp_theta_stc) + Yp[i]*std::sin(tmp_theta_stc)) + Real(2.)*I*(Xp[i]*std::cos(tmp_theta_stc) + Yp[i]*std::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 * MathFunc::exp( - stc_exponent ); diff --git a/Source/Make.WarpX b/Source/Make.WarpX index 1dfffbbb6..13d4f4554 100644 --- a/Source/Make.WarpX +++ b/Source/Make.WarpX @@ -175,11 +175,7 @@ ifeq ($(USE_HDF5),TRUE) LIBRARY_LOCATIONS += $(HDF5_HOME)/lib endif DEFINES += -DWARPX_USE_HDF5 -<<<<<<< HEAD - LIBRARIES += -lhdf5 -lz -======= libraries += -lhdf5 -lz ->>>>>>> dev endif # job_info support diff --git a/Source/Parser/GpuParser.H b/Source/Parser/GpuParser.H index e49671e06..c158ee314 100644 --- a/Source/Parser/GpuParser.H +++ b/Source/Parser/GpuParser.H @@ -16,16 +16,16 @@ public: void clear (); AMREX_GPU_HOST_DEVICE - double - operator() (double x, double y, double z) const noexcept + amrex::Real + operator() (amrex::Real x, amrex::Real y, amrex::Real z) const noexcept { #ifdef AMREX_USE_GPU #ifdef AMREX_DEVICE_COMPILE // WarpX compiled for GPU, function compiled for __device__ // the 3D position of each particle is stored in shared memory. - amrex::Gpu::SharedMemory<double> gsm; - double* p = gsm.dataPtr(); + amrex::Gpu::SharedMemory<amrex::Real> gsm; + amrex::Real* p = gsm.dataPtr(); int tid = threadIdx.x + threadIdx.y*blockDim.x + threadIdx.z*(blockDim.x*blockDim.y); p[tid*3] = x; p[tid*3+1] = y; diff --git a/Source/Parser/WarpXParser.H b/Source/Parser/WarpXParser.H index ffa61e457..8c1d854d8 100644 --- a/Source/Parser/WarpXParser.H +++ b/Source/Parser/WarpXParser.H @@ -6,6 +6,8 @@ #include <string> #include <set> +#include <AMReX_REAL.H> + #include "wp_parser_c.h" #include "wp_parser_y.h" @@ -23,15 +25,15 @@ public: ~WarpXParser (); void define (std::string const& func_body); - void setConstant (std::string const& name, double c); + void setConstant (std::string const& name, amrex::Real c); // // Option 1: Register every variable to an address provided. // Assign values to external variables. // Call eval(). - void registerVariable (std::string const& name, double& var); + void registerVariable (std::string const& name, amrex::Real& var); // - inline double eval () const noexcept; + inline amrex::Real eval () const noexcept; // // Option 2: Register all variables at once. Parser will create @@ -40,7 +42,7 @@ public: void registerVariables (std::vector<std::string> const& names); // template <typename T, typename... Ts> inline - double eval (T x, Ts... yz) const noexcept; + amrex::Real eval (T x, Ts... yz) const noexcept; void print () const; @@ -54,23 +56,23 @@ private: void clear (); template <typename T> inline - void unpack (double* p, T x) const noexcept; + void unpack (amrex::Real* p, T x) const noexcept; template <typename T, typename... Ts> inline - void unpack (double* p, T x, Ts... yz) const noexcept; + void unpack (amrex::Real* p, T x, Ts... yz) const noexcept; std::string m_expression; #ifdef _OPENMP std::vector<struct wp_parser*> m_parser; - mutable std::vector<std::array<double,16> > m_variables; + mutable std::vector<std::array<amrex::Real,16> > m_variables; #else struct wp_parser* m_parser = nullptr; - mutable std::array<double,16> m_variables; + mutable std::array<amrex::Real,16> m_variables; #endif }; inline -double +amrex::Real WarpXParser::eval () const noexcept { #ifdef _OPENMP @@ -82,7 +84,7 @@ WarpXParser::eval () const noexcept template <typename T, typename... Ts> inline -double +amrex::Real WarpXParser::eval (T x, Ts... yz) const noexcept { #ifdef _OPENMP @@ -96,7 +98,7 @@ WarpXParser::eval (T x, Ts... yz) const noexcept template <typename T> inline void -WarpXParser::unpack (double* p, T x) const noexcept +WarpXParser::unpack (amrex::Real* p, T x) const noexcept { *p = x; } @@ -104,7 +106,7 @@ WarpXParser::unpack (double* p, T x) const noexcept template <typename T, typename... Ts> inline void -WarpXParser::unpack (double* p, T x, Ts... yz) const noexcept +WarpXParser::unpack (amrex::Real* p, T x, Ts... yz) const noexcept { *p++ = x; unpack(p, yz...); diff --git a/Source/Parser/WarpXParser.cpp b/Source/Parser/WarpXParser.cpp index 3237086f2..ced536327 100644 --- a/Source/Parser/WarpXParser.cpp +++ b/Source/Parser/WarpXParser.cpp @@ -69,7 +69,7 @@ WarpXParser::clear () } void -WarpXParser::registerVariable (std::string const& name, double& var) +WarpXParser::registerVariable (std::string const& name, amrex::Real& var) { // We assume this is called inside OMP parallel region #ifdef _OPENMP @@ -105,7 +105,7 @@ WarpXParser::registerVariables (std::vector<std::string> const& names) } void -WarpXParser::setConstant (std::string const& name, double c) +WarpXParser::setConstant (std::string const& name, amrex::Real c) { #ifdef _OPENMP diff --git a/Source/Parser/wp_parser.tab.c b/Source/Parser/wp_parser.tab.c index 3981894a5..0f7c2403d 100644 --- a/Source/Parser/wp_parser.tab.c +++ b/Source/Parser/wp_parser.tab.c @@ -138,7 +138,7 @@ union YYSTYPE #line 19 "wp_parser.y" /* yacc.c:352 */ struct wp_node* n; - double d; + amrex_real d; struct wp_symbol* s; enum wp_f1_t f1; enum wp_f2_t f2; diff --git a/Source/Parser/wp_parser.tab.h b/Source/Parser/wp_parser.tab.h index b50516808..0c859fc03 100644 --- a/Source/Parser/wp_parser.tab.h +++ b/Source/Parser/wp_parser.tab.h @@ -75,7 +75,7 @@ union YYSTYPE #line 19 "wp_parser.y" /* yacc.c:1921 */ struct wp_node* n; - double d; + amrex_real d; struct wp_symbol* s; enum wp_f1_t f1; enum wp_f2_t f2; diff --git a/Source/Parser/wp_parser.y b/Source/Parser/wp_parser.y index 453eda1cd..809dbfa5e 100644 --- a/Source/Parser/wp_parser.y +++ b/Source/Parser/wp_parser.y @@ -18,7 +18,7 @@ */ %union { struct wp_node* n; - double d; + amrex_real d; struct wp_symbol* s; enum wp_f1_t f1; enum wp_f2_t f2; diff --git a/Source/Parser/wp_parser_c.h b/Source/Parser/wp_parser_c.h index 3aafdec65..970d6b355 100644 --- a/Source/Parser/wp_parser_c.h +++ b/Source/Parser/wp_parser_c.h @@ -4,6 +4,7 @@ #include "wp_parser_y.h" #include <AMReX_GpuQualifiers.H> #include <AMReX_Extension.H> +#include <AMReX_REAL.H> #ifdef __cplusplus extern "C" { @@ -21,15 +22,15 @@ extern "C" { #include <string> AMREX_GPU_HOST_DEVICE -inline double +inline amrex_real wp_ast_eval (struct wp_node* node) { - double result; + amrex_real result; #ifdef AMREX_DEVICE_COMPILE - extern __shared__ double extern_xyz[]; + extern __shared__ amrex_real extern_xyz[]; int tid = threadIdx.x + threadIdx.y*blockDim.x + threadIdx.z*(blockDim.x*blockDim.y); - double* x = extern_xyz + tid*3; + amrex_real* x = extern_xyz + tid*3; #endif switch (node->type) diff --git a/Source/Parser/wp_parser_y.c b/Source/Parser/wp_parser_y.c index 259f9368b..b71b42638 100644 --- a/Source/Parser/wp_parser_y.c +++ b/Source/Parser/wp_parser_y.c @@ -6,8 +6,6 @@ #include "wp_parser_y.h" #include "wp_parser.tab.h" -#include <AMReX_GpuQualifiers.H> - static struct wp_node* wp_root = NULL; /* This is called by a bison rule to store the original AST in a @@ -21,7 +19,7 @@ wp_defexpr (struct wp_node* body) } struct wp_node* -wp_newnumber (double d) +wp_newnumber (amrex_real d) { struct wp_number* r = (struct wp_number*) malloc(sizeof(struct wp_number)); r->type = WP_NUMBER; @@ -154,8 +152,8 @@ wp_parser_dup (struct wp_parser* source) } AMREX_GPU_HOST_DEVICE -double -wp_call_f1 (enum wp_f1_t type, double a) +amrex_real +wp_call_f1 (enum wp_f1_t type, amrex_real a) { switch (type) { case WP_SQRT: return sqrt(a); @@ -185,8 +183,8 @@ wp_call_f1 (enum wp_f1_t type, double a) } AMREX_GPU_HOST_DEVICE -double -wp_call_f2 (enum wp_f2_t type, double a, double b) +amrex_real +wp_call_f2 (enum wp_f2_t type, amrex_real a, amrex_real b) { switch (type) { case WP_POW: @@ -356,13 +354,13 @@ wp_parser_ast_dup (struct wp_parser* my_parser, struct wp_node* node, int move) #define WP_MOVEUP_R(node, v) \ struct wp_node* n = node->r->r; \ - double* p = node->r->rip.p; \ + amrex_real* p = node->r->rip.p; \ node->r = n; \ node->lvp.v = v; \ node->rip.p = p; #define WP_MOVEUP_L(node, v) \ struct wp_node* n = node->l->r; \ - double* p = node->l->rip.p; \ + amrex_real* p = node->l->rip.p; \ node->r = n; \ node->lvp.v = v; \ node->rip.p = p; @@ -392,7 +390,7 @@ wp_ast_optimize (struct wp_node* node) if (node->l->type == WP_NUMBER && node->r->type == WP_NUMBER) { - double v = ((struct wp_number*)(node->l))->value + amrex_real v = ((struct wp_number*)(node->l))->value + ((struct wp_number*)(node->r))->value; ((struct wp_number*)node)->type = WP_NUMBER; ((struct wp_number*)node)->value = v; @@ -422,28 +420,28 @@ wp_ast_optimize (struct wp_node* node) else if (node->l->type == WP_NUMBER && node->r->type == WP_ADD_VP) { - double v = ((struct wp_number*)(node->l))->value + WP_EVAL_R(node); + amrex_real v = ((struct wp_number*)(node->l))->value + WP_EVAL_R(node); WP_MOVEUP_R(node, v); node->type = WP_ADD_VP; } else if (node->l->type == WP_NUMBER && node->r->type == WP_SUB_VP) { - double v = ((struct wp_number*)(node->l))->value + WP_EVAL_R(node); + amrex_real v = ((struct wp_number*)(node->l))->value + WP_EVAL_R(node); WP_MOVEUP_R(node, v); node->type = WP_SUB_VP; } else if (node->l->type == WP_ADD_VP && node->r->type == WP_NUMBER) { - double v = WP_EVAL_L(node) + ((struct wp_number*)(node->r))->value; + amrex_real v = WP_EVAL_L(node) + ((struct wp_number*)(node->r))->value; WP_MOVEUP_L(node, v); node->type = WP_ADD_VP; } else if (node->l->type == WP_SUB_VP && node->r->type == WP_NUMBER) { - double v = WP_EVAL_L(node) + ((struct wp_number*)(node->r))->value; + amrex_real v = WP_EVAL_L(node) + ((struct wp_number*)(node->r))->value; WP_MOVEUP_L(node, v); node->type = WP_SUB_VP; } @@ -455,7 +453,7 @@ wp_ast_optimize (struct wp_node* node) if (node->l->type == WP_NUMBER && node->r->type == WP_NUMBER) { - double v = ((struct wp_number*)(node->l))->value + amrex_real v = ((struct wp_number*)(node->l))->value - ((struct wp_number*)(node->r))->value; ((struct wp_number*)node)->type = WP_NUMBER; ((struct wp_number*)node)->value = v; @@ -485,28 +483,28 @@ wp_ast_optimize (struct wp_node* node) else if (node->l->type == WP_NUMBER && node->r->type == WP_ADD_VP) { - double v = ((struct wp_number*)(node->l))->value - WP_EVAL_R(node); + amrex_real v = ((struct wp_number*)(node->l))->value - WP_EVAL_R(node); WP_MOVEUP_R(node, v); node->type = WP_SUB_VP; } else if (node->l->type == WP_NUMBER && node->r->type == WP_SUB_VP) { - double v = ((struct wp_number*)(node->l))->value - WP_EVAL_R(node); + amrex_real v = ((struct wp_number*)(node->l))->value - WP_EVAL_R(node); WP_MOVEUP_R(node, v); node->type = WP_ADD_VP; } else if (node->l->type == WP_ADD_VP && node->r->type == WP_NUMBER) { - double v = WP_EVAL_L(node) - ((struct wp_number*)(node->r))->value; + amrex_real v = WP_EVAL_L(node) - ((struct wp_number*)(node->r))->value; WP_MOVEUP_L(node, v); node->type = WP_ADD_VP; } else if (node->l->type == WP_SUB_VP && node->r->type == WP_NUMBER) { - double v = WP_EVAL_L(node) - ((struct wp_number*)(node->r))->value; + amrex_real v = WP_EVAL_L(node) - ((struct wp_number*)(node->r))->value; WP_MOVEUP_L(node, v); node->type = WP_SUB_VP; } @@ -518,7 +516,7 @@ wp_ast_optimize (struct wp_node* node) if (node->l->type == WP_NUMBER && node->r->type == WP_NUMBER) { - double v = ((struct wp_number*)(node->l))->value + amrex_real v = ((struct wp_number*)(node->l))->value * ((struct wp_number*)(node->r))->value; ((struct wp_number*)node)->type = WP_NUMBER; ((struct wp_number*)node)->value = v; @@ -548,28 +546,28 @@ wp_ast_optimize (struct wp_node* node) else if (node->l->type == WP_NUMBER && node->r->type == WP_MUL_VP) { - double v = ((struct wp_number*)(node->l))->value * WP_EVAL_R(node); + amrex_real v = ((struct wp_number*)(node->l))->value * WP_EVAL_R(node); WP_MOVEUP_R(node, v); node->type = WP_MUL_VP; } else if (node->l->type == WP_NUMBER && node->r->type == WP_DIV_VP) { - double v = ((struct wp_number*)(node->l))->value * WP_EVAL_R(node); + amrex_real v = ((struct wp_number*)(node->l))->value * WP_EVAL_R(node); WP_MOVEUP_R(node, v); node->type = WP_DIV_VP; } else if (node->l->type == WP_MUL_VP && node->r->type == WP_NUMBER) { - double v = WP_EVAL_L(node) * ((struct wp_number*)(node->r))->value; + amrex_real v = WP_EVAL_L(node) * ((struct wp_number*)(node->r))->value; WP_MOVEUP_L(node, v); node->type = WP_MUL_VP; } else if (node->l->type == WP_DIV_VP && node->r->type == WP_NUMBER) { - double v = WP_EVAL_L(node) * ((struct wp_number*)(node->r))->value; + amrex_real v = WP_EVAL_L(node) * ((struct wp_number*)(node->r))->value; WP_MOVEUP_L(node, v); node->type = WP_DIV_VP; } @@ -581,7 +579,7 @@ wp_ast_optimize (struct wp_node* node) if (node->l->type == WP_NUMBER && node->r->type == WP_NUMBER) { - double v = ((struct wp_number*)(node->l))->value + amrex_real v = ((struct wp_number*)(node->l))->value / ((struct wp_number*)(node->r))->value; ((struct wp_number*)node)->type = WP_NUMBER; ((struct wp_number*)node)->value = v; @@ -611,28 +609,28 @@ wp_ast_optimize (struct wp_node* node) else if (node->l->type == WP_NUMBER && node->r->type == WP_MUL_VP) { - double v = ((struct wp_number*)(node->l))->value / WP_EVAL_R(node); + amrex_real v = ((struct wp_number*)(node->l))->value / WP_EVAL_R(node); WP_MOVEUP_R(node, v); node->type = WP_DIV_VP; } else if (node->l->type == WP_NUMBER && node->r->type == WP_DIV_VP) { - double v = ((struct wp_number*)(node->l))->value / WP_EVAL_R(node); + amrex_real v = ((struct wp_number*)(node->l))->value / WP_EVAL_R(node); WP_MOVEUP_R(node, v); node->type = WP_MUL_VP; } else if (node->l->type == WP_MUL_VP && node->r->type == WP_NUMBER) { - double v = WP_EVAL_L(node) / ((struct wp_number*)(node->r))->value; + amrex_real v = WP_EVAL_L(node) / ((struct wp_number*)(node->r))->value; WP_MOVEUP_L(node, v); node->type = WP_MUL_VP; } else if (node->l->type == WP_DIV_VP && node->r->type == WP_NUMBER) { - double v = WP_EVAL_L(node) / ((struct wp_number*)(node->r))->value; + amrex_real v = WP_EVAL_L(node) / ((struct wp_number*)(node->r))->value; WP_MOVEUP_L(node, v); node->type = WP_DIV_VP; } @@ -641,7 +639,7 @@ wp_ast_optimize (struct wp_node* node) wp_ast_optimize(node->l); if (node->l->type == WP_NUMBER) { - double v = -((struct wp_number*)(node->l))->value; + amrex_real v = -((struct wp_number*)(node->l))->value; ((struct wp_number*)node)->type = WP_NUMBER; ((struct wp_number*)node)->value = v; } @@ -675,7 +673,7 @@ wp_ast_optimize (struct wp_node* node) wp_ast_optimize(node->l); if (node->l->type == WP_NUMBER) { - double v = wp_call_f1 + amrex_real v = wp_call_f1 (((struct wp_f1*)node)->ftype, ((struct wp_number*)(((struct wp_f1*)node)->l))->value); ((struct wp_number*)node)->type = WP_NUMBER; @@ -688,7 +686,7 @@ wp_ast_optimize (struct wp_node* node) if (node->l->type == WP_NUMBER && node->r->type == WP_NUMBER) { - double v = wp_call_f2 + amrex_real v = wp_call_f2 (((struct wp_f2*)node)->ftype, ((struct wp_number*)(((struct wp_f2*)node)->l))->value, ((struct wp_number*)(((struct wp_f2*)node)->r))->value); @@ -698,7 +696,7 @@ wp_ast_optimize (struct wp_node* node) else if (node->r->type == WP_NUMBER && ((struct wp_f2*)node)->ftype == WP_POW) { struct wp_node* n = node->l; - double v = ((struct wp_number*)(node->r))->value; + amrex_real v = ((struct wp_number*)(node->r))->value; if (-3.0 == v) { ((struct wp_f1*)node)->type = WP_F1; ((struct wp_f1*)node)->l = n; @@ -733,7 +731,7 @@ wp_ast_optimize (struct wp_node* node) wp_ast_optimize(node->r); if (node->r->type == WP_NUMBER) { - double v = node->lvp.v + ((struct wp_number*)(node->r))->value; + amrex_real v = node->lvp.v + ((struct wp_number*)(node->r))->value; ((struct wp_number*)node)->type = WP_NUMBER; ((struct wp_number*)node)->value = v; } @@ -742,7 +740,7 @@ wp_ast_optimize (struct wp_node* node) wp_ast_optimize(node->r); if (node->r->type == WP_NUMBER) { - double v = node->lvp.v - ((struct wp_number*)(node->r))->value; + amrex_real v = node->lvp.v - ((struct wp_number*)(node->r))->value; ((struct wp_number*)node)->type = WP_NUMBER; ((struct wp_number*)node)->value = v; } @@ -751,7 +749,7 @@ wp_ast_optimize (struct wp_node* node) wp_ast_optimize(node->r); if (node->r->type == WP_NUMBER) { - double v = node->lvp.v * ((struct wp_number*)(node->r))->value; + amrex_real v = node->lvp.v * ((struct wp_number*)(node->r))->value; ((struct wp_number*)node)->type = WP_NUMBER; ((struct wp_number*)node)->value = v; } @@ -760,7 +758,7 @@ wp_ast_optimize (struct wp_node* node) wp_ast_optimize(node->r); if (node->r->type == WP_NUMBER) { - double v = node->lvp.v / ((struct wp_number*)(node->r))->value; + amrex_real v = node->lvp.v / ((struct wp_number*)(node->r))->value; ((struct wp_number*)node)->type = WP_NUMBER; ((struct wp_number*)node)->value = v; } @@ -769,7 +767,7 @@ wp_ast_optimize (struct wp_node* node) wp_ast_optimize(node->l); if (node->l->type == WP_NUMBER) { - double v = -((struct wp_number*)(node->l))->value; + amrex_real v = -((struct wp_number*)(node->l))->value; ((struct wp_number*)node)->type = WP_NUMBER; ((struct wp_number*)node)->value = v; } @@ -938,7 +936,7 @@ wp_ast_print (struct wp_node* node) } void -wp_ast_regvar (struct wp_node* node, char const* name, double* p) +wp_ast_regvar (struct wp_node* node, char const* name, amrex_real* p) { switch (node->type) { @@ -1047,7 +1045,7 @@ wp_ast_regvar_gpu (struct wp_node* node, char const* name, int i) } } -void wp_ast_setconst (struct wp_node* node, char const* name, double c) +void wp_ast_setconst (struct wp_node* node, char const* name, amrex_real c) { switch (node->type) { @@ -1099,7 +1097,7 @@ void wp_ast_setconst (struct wp_node* node, char const* name, double c) } void -wp_parser_regvar (struct wp_parser* parser, char const* name, double* p) +wp_parser_regvar (struct wp_parser* parser, char const* name, amrex_real* p) { wp_ast_regvar(parser->ast, name, p); } @@ -1111,7 +1109,7 @@ wp_parser_regvar_gpu (struct wp_parser* parser, char const* name, int i) } void -wp_parser_setconst (struct wp_parser* parser, char const* name, double c) +wp_parser_setconst (struct wp_parser* parser, char const* name, amrex_real c) { wp_ast_setconst(parser->ast, name, c); wp_ast_optimize(parser->ast); diff --git a/Source/Parser/wp_parser_y.h b/Source/Parser/wp_parser_y.h index 8c9f8e4e4..d83815090 100644 --- a/Source/Parser/wp_parser_y.h +++ b/Source/Parser/wp_parser_y.h @@ -2,6 +2,7 @@ #define WP_PARSER_Y_H_ #include <AMReX_GpuQualifiers.H> +#include <AMReX_REAL.H> #ifdef __cplusplus #include <cstdlib> @@ -77,11 +78,11 @@ enum wp_node_t { union wp_ip { int i; - double* p; + amrex_real* p; }; union wp_vp { - double v; + amrex_real v; union wp_ip ip; }; @@ -95,7 +96,7 @@ struct wp_node { struct wp_number { enum wp_node_t type; - double value; + amrex_real value; }; struct wp_symbol { @@ -122,7 +123,7 @@ struct wp_f2 { /* Builtin functions with two arguments */ /* These functions are used in bison rules to generate the original * AST. */ void wp_defexpr (struct wp_node* body); -struct wp_node* wp_newnumber (double d); +struct wp_node* wp_newnumber (amrex_real d); struct wp_symbol* wp_makesymbol (char* name); struct wp_node* wp_newsymbol (struct wp_symbol* sym); struct wp_node* wp_newnode (enum wp_node_t type, struct wp_node* l, @@ -153,20 +154,20 @@ void wp_parser_delete (struct wp_parser* parser); struct wp_parser* wp_parser_dup (struct wp_parser* source); struct wp_node* wp_parser_ast_dup (struct wp_parser* parser, struct wp_node* src, int move); -void wp_parser_regvar (struct wp_parser* parser, char const* name, double* p); +void wp_parser_regvar (struct wp_parser* parser, char const* name, amrex_real* p); void wp_parser_regvar_gpu (struct wp_parser* parser, char const* name, int i); -void wp_parser_setconst (struct wp_parser* parser, char const* name, double c); +void wp_parser_setconst (struct wp_parser* parser, char const* name, amrex_real c); /* We need to walk the tree in these functions */ void wp_ast_optimize (struct wp_node* node); size_t wp_ast_size (struct wp_node* node); void wp_ast_print (struct wp_node* node); -void wp_ast_regvar (struct wp_node* node, char const* name, double* p); +void wp_ast_regvar (struct wp_node* node, char const* name, amrex_real* p); void wp_ast_regvar_gpu (struct wp_node* node, char const* name, int i); -void wp_ast_setconst (struct wp_node* node, char const* name, double c); +void wp_ast_setconst (struct wp_node* node, char const* name, amrex_real c); -AMREX_GPU_HOST_DEVICE double wp_call_f1 (enum wp_f1_t type, double a); -AMREX_GPU_HOST_DEVICE double wp_call_f2 (enum wp_f2_t type, double a, double b); +AMREX_GPU_HOST_DEVICE amrex_real wp_call_f1 (enum wp_f1_t type, amrex_real a); +AMREX_GPU_HOST_DEVICE amrex_real wp_call_f2 (enum wp_f2_t type, amrex_real a, amrex_real b); #ifdef __cplusplus } diff --git a/Source/Particles/Deposition/ChargeDeposition.H b/Source/Particles/Deposition/ChargeDeposition.H index b9210e67c..eec407a2b 100755 --- a/Source/Particles/Deposition/ChargeDeposition.H +++ b/Source/Particles/Deposition/ChargeDeposition.H @@ -18,10 +18,10 @@ * /param q : species charge. */ template <int depos_order> -void doChargeDepositionShapeN(const amrex::Real * const xp, - const amrex::Real * const yp, - const amrex::Real * const zp, - const amrex::Real * const wp, +void doChargeDepositionShapeN(const amrex::ParticleReal * const xp, + const amrex::ParticleReal * const yp, + const amrex::ParticleReal * const zp, + const amrex::ParticleReal * const wp, const int * const ion_lev, const amrex::Array4<amrex::Real>& rho_arr, const long np_to_depose, diff --git a/Source/Particles/Deposition/CurrentDeposition.H b/Source/Particles/Deposition/CurrentDeposition.H index 7a96dab9a..6da0f1155 100644 --- a/Source/Particles/Deposition/CurrentDeposition.H +++ b/Source/Particles/Deposition/CurrentDeposition.H @@ -24,13 +24,13 @@ * /param q : species charge. */ template <int depos_order> -void doDepositionShapeN(const amrex::Real * const xp, - const amrex::Real * const yp, - const amrex::Real * const zp, - const amrex::Real * const wp, - const amrex::Real * const uxp, - const amrex::Real * const uyp, - const amrex::Real * const uzp, +void doDepositionShapeN(const amrex::ParticleReal * const xp, + const amrex::ParticleReal * const yp, + const amrex::ParticleReal * const zp, + const amrex::ParticleReal * const wp, + const amrex::ParticleReal * const uxp, + const amrex::ParticleReal * const uyp, + const amrex::ParticleReal * const uzp, const int * const ion_lev, const amrex::Array4<amrex::Real>& jx_arr, const amrex::Array4<amrex::Real>& jy_arr, @@ -189,13 +189,13 @@ void doDepositionShapeN(const amrex::Real * const xp, * \param n_rz_azimuthal_modes: Number of azimuthal modes when using RZ geometry */ template <int depos_order> -void doEsirkepovDepositionShapeN (const amrex::Real * const xp, - const amrex::Real * const yp, - const amrex::Real * const zp, - const amrex::Real * const wp, - const amrex::Real * const uxp, - const amrex::Real * const uyp, - const amrex::Real * const uzp, +void doEsirkepovDepositionShapeN (const amrex::ParticleReal * const xp, + const amrex::ParticleReal * const yp, + const amrex::ParticleReal * const zp, + const amrex::ParticleReal * const wp, + const amrex::ParticleReal * const uxp, + const amrex::ParticleReal * const uyp, + const amrex::ParticleReal * const uzp, const int * ion_lev, const amrex::Array4<amrex::Real>& Jx_arr, const amrex::Array4<amrex::Real>& Jy_arr, @@ -397,7 +397,7 @@ void doEsirkepovDepositionShapeN (const amrex::Real * const xp, Complex xy_mid = xy_mid0; // Throughout the following loop, xy_mid takes the value e^{i m theta} for (int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) { // The factor 2 comes from the normalization of the modes - const Complex djr_cmplx = 2.*sdxi*xy_mid; + const Complex djr_cmplx = amrex::Real(2.)*sdxi*xy_mid; amrex::Gpu::Atomic::Add( &Jx_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode-1), djr_cmplx.real()); amrex::Gpu::Atomic::Add( &Jx_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode), djr_cmplx.imag()); xy_mid = xy_mid*xy_mid0; @@ -418,7 +418,7 @@ void doEsirkepovDepositionShapeN (const amrex::Real * const xp, for (int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) { // The factor 2 comes from the normalization of the modes // The minus sign comes from the different convention with respect to Davidson et al. - const Complex djt_cmplx = -2.*I*(i_new-1 + i + xmin*dxi)*wq*invdtdx/(double)imode* + const Complex djt_cmplx = -amrex::Real(2.)*I*(i_new-1 + i + xmin*dxi)*wq*invdtdx/(amrex::Real)imode* (sx_new[i]*sz_new[k]*(xy_new - xy_mid) + sx_old[i]*sz_old[k]*(xy_mid - xy_old)); amrex::Gpu::Atomic::Add( &Jy_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode-1), djt_cmplx.real()); amrex::Gpu::Atomic::Add( &Jy_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode), djt_cmplx.imag()); @@ -438,7 +438,7 @@ void doEsirkepovDepositionShapeN (const amrex::Real * const xp, Complex xy_mid = xy_mid0; // Throughout the following loop, xy_mid takes the value e^{i m theta} for (int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) { // The factor 2 comes from the normalization of the modes - const Complex djz_cmplx = 2.*sdzk*xy_mid; + const Complex djz_cmplx = amrex::Real(2.)*sdzk*xy_mid; amrex::Gpu::Atomic::Add( &Jz_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode-1), djz_cmplx.real()); amrex::Gpu::Atomic::Add( &Jz_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode), djz_cmplx.imag()); xy_mid = xy_mid*xy_mid0; diff --git a/Source/Particles/Gather/FieldGather.H b/Source/Particles/Gather/FieldGather.H index 6727b0aa9..c5ec6fb5b 100644 --- a/Source/Particles/Gather/FieldGather.H +++ b/Source/Particles/Gather/FieldGather.H @@ -19,12 +19,12 @@ * \param n_rz_azimuthal_modes: Number of azimuthal modes when using RZ geometry */ template <int depos_order, int lower_in_v> -void doGatherShapeN(const amrex::Real * const xp, - const amrex::Real * const yp, - const amrex::Real * const zp, - amrex::Real * const Exp, amrex::Real * const Eyp, - amrex::Real * const Ezp, amrex::Real * const Bxp, - amrex::Real * const Byp, amrex::Real * const Bzp, +void doGatherShapeN(const amrex::ParticleReal * const xp, + const amrex::ParticleReal * const yp, + const amrex::ParticleReal * const zp, + amrex::ParticleReal * const Exp, amrex::ParticleReal * const Eyp, + amrex::ParticleReal * const Ezp, amrex::ParticleReal * const Bxp, + amrex::ParticleReal * const Byp, amrex::ParticleReal * const Bzp, const amrex::Array4<const amrex::Real>& ex_arr, const amrex::Array4<const amrex::Real>& ey_arr, const amrex::Array4<const amrex::Real>& ez_arr, diff --git a/Source/Particles/Make.package b/Source/Particles/Make.package index 95f36cf65..c9d520873 100644 --- a/Source/Particles/Make.package +++ b/Source/Particles/Make.package @@ -18,6 +18,7 @@ CEXE_headers += ShapeFactors.H include $(WARPX_HOME)/Source/Particles/Pusher/Make.package include $(WARPX_HOME)/Source/Particles/Deposition/Make.package include $(WARPX_HOME)/Source/Particles/Gather/Make.package +include $(WARPX_HOME)/Source/Particles/Sorting/Make.package INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Particles VPATH_LOCATIONS += $(WARPX_HOME)/Source/Particles diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index bb795465e..715c97b99 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -545,12 +545,12 @@ namespace WarpXParticleContainer::ParticleType* particles_source = ptile_source.GetArrayOfStructs()().data(); // --- source SoA particle data auto& soa_source = ptile_source.GetStructOfArrays(); - GpuArray<Real*,PIdx::nattribs> attribs_source; + GpuArray<ParticleReal*,PIdx::nattribs> attribs_source; for (int ia = 0; ia < PIdx::nattribs; ++ia) { attribs_source[ia] = soa_source.GetRealData(ia).data(); } // --- source runtime attribs - GpuArray<Real*,3> runtime_uold_source; + GpuArray<ParticleReal*,3> runtime_uold_source; // Prepare arrays for boosted frame diagnostics. runtime_uold_source[0] = soa_source.GetRealData(PIdx::ux).data(); runtime_uold_source[1] = soa_source.GetRealData(PIdx::uy).data(); @@ -590,13 +590,13 @@ namespace WarpXParticleContainer::ParticleType* particles_product = ptile_product.GetArrayOfStructs()().data() + np_product_old; // --- product SoA particle data auto& soa_product = ptile_product.GetStructOfArrays(); - GpuArray<Real*,PIdx::nattribs> attribs_product; + GpuArray<ParticleReal*,PIdx::nattribs> attribs_product; for (int ia = 0; ia < PIdx::nattribs; ++ia) { // First element is the first newly-created product particle attribs_product[ia] = soa_product.GetRealData(ia).data() + np_product_old; } // --- product runtime attribs - GpuArray<Real*,6> runtime_attribs_product; + GpuArray<ParticleReal*,6> runtime_attribs_product; bool do_boosted_product = WarpX::do_boosted_frame_diagnostic && pc_product->DoBoostedFrameDiags(); if (do_boosted_product) { diff --git a/Source/Particles/PhotonParticleContainer.H b/Source/Particles/PhotonParticleContainer.H index 3ac304bdc..a6ffd1d76 100644 --- a/Source/Particles/PhotonParticleContainer.H +++ b/Source/Particles/PhotonParticleContainer.H @@ -41,9 +41,9 @@ public: DtType a_dt_type=DtType::Full) override; virtual void PushPX(WarpXParIter& pti, - amrex::Cuda::ManagedDeviceVector<amrex::Real>& xp, - amrex::Cuda::ManagedDeviceVector<amrex::Real>& yp, - amrex::Cuda::ManagedDeviceVector<amrex::Real>& zp, + amrex::Cuda::ManagedDeviceVector<amrex::ParticleReal>& xp, + amrex::Cuda::ManagedDeviceVector<amrex::ParticleReal>& yp, + amrex::Cuda::ManagedDeviceVector<amrex::ParticleReal>& zp, amrex::Real dt, DtType a_dt_type=DtType::Full) override; diff --git a/Source/Particles/PhotonParticleContainer.cpp b/Source/Particles/PhotonParticleContainer.cpp index 9de441e5c..4a75ec9f3 100644 --- a/Source/Particles/PhotonParticleContainer.cpp +++ b/Source/Particles/PhotonParticleContainer.cpp @@ -34,27 +34,27 @@ void PhotonParticleContainer::InitData() void PhotonParticleContainer::PushPX(WarpXParIter& pti, - Cuda::ManagedDeviceVector<Real>& xp, - Cuda::ManagedDeviceVector<Real>& yp, - Cuda::ManagedDeviceVector<Real>& zp, + Cuda::ManagedDeviceVector<ParticleReal>& xp, + Cuda::ManagedDeviceVector<ParticleReal>& yp, + Cuda::ManagedDeviceVector<ParticleReal>& zp, Real dt, DtType a_dt_type) { // This wraps the momentum and position advance so that inheritors can modify the call. auto& attribs = pti.GetAttribs(); // Extract pointers to the different particle quantities - Real* const AMREX_RESTRICT x = xp.dataPtr(); - Real* const AMREX_RESTRICT y = yp.dataPtr(); - Real* const AMREX_RESTRICT z = zp.dataPtr(); - Real* const AMREX_RESTRICT ux = attribs[PIdx::ux].dataPtr(); - Real* const AMREX_RESTRICT uy = attribs[PIdx::uy].dataPtr(); - Real* const AMREX_RESTRICT uz = attribs[PIdx::uz].dataPtr(); - const Real* const AMREX_RESTRICT Ex = attribs[PIdx::Ex].dataPtr(); - const Real* const AMREX_RESTRICT Ey = attribs[PIdx::Ey].dataPtr(); - const Real* const AMREX_RESTRICT Ez = attribs[PIdx::Ez].dataPtr(); - const Real* const AMREX_RESTRICT Bx = attribs[PIdx::Bx].dataPtr(); - const Real* const AMREX_RESTRICT By = attribs[PIdx::By].dataPtr(); - const Real* const AMREX_RESTRICT Bz = attribs[PIdx::Bz].dataPtr(); + ParticleReal* const AMREX_RESTRICT x = xp.dataPtr(); + ParticleReal* const AMREX_RESTRICT y = yp.dataPtr(); + ParticleReal* const AMREX_RESTRICT z = zp.dataPtr(); + ParticleReal* const AMREX_RESTRICT ux = attribs[PIdx::ux].dataPtr(); + ParticleReal* const AMREX_RESTRICT uy = attribs[PIdx::uy].dataPtr(); + ParticleReal* const AMREX_RESTRICT uz = attribs[PIdx::uz].dataPtr(); + const ParticleReal* const AMREX_RESTRICT Ex = attribs[PIdx::Ex].dataPtr(); + const ParticleReal* const AMREX_RESTRICT Ey = attribs[PIdx::Ey].dataPtr(); + const ParticleReal* const AMREX_RESTRICT Ez = attribs[PIdx::Ez].dataPtr(); + const ParticleReal* const AMREX_RESTRICT Bx = attribs[PIdx::Bx].dataPtr(); + const ParticleReal* const AMREX_RESTRICT By = attribs[PIdx::By].dataPtr(); + const ParticleReal* const AMREX_RESTRICT Bz = attribs[PIdx::Bz].dataPtr(); if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index ace1ec7f8..b558323a3 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -1,12 +1,13 @@ #ifndef WARPX_PhysicalParticleContainer_H_ #define WARPX_PhysicalParticleContainer_H_ -#include <map> +#include <PlasmaInjector.H> +#include <WarpXParticleContainer.H> +#include <NCIGodfreyFilter.H> #include <AMReX_IArrayBox.H> -#include <PlasmaInjector.H> -#include <WarpXParticleContainer.H> +#include <map> class PhysicalParticleContainer : public WarpXParticleContainer @@ -87,9 +88,9 @@ public: DtType a_dt_type=DtType::Full) override; virtual void PushPX(WarpXParIter& pti, - amrex::Cuda::ManagedDeviceVector<amrex::Real>& xp, - amrex::Cuda::ManagedDeviceVector<amrex::Real>& yp, - amrex::Cuda::ManagedDeviceVector<amrex::Real>& zp, + amrex::Cuda::ManagedDeviceVector<amrex::ParticleReal>& xp, + amrex::Cuda::ManagedDeviceVector<amrex::ParticleReal>& yp, + amrex::Cuda::ManagedDeviceVector<amrex::ParticleReal>& zp, amrex::Real dt, DtType a_dt_type=DtType::Full); virtual void PushP (int lev, amrex::Real dt, @@ -100,8 +101,21 @@ public: const amrex::MultiFab& By, const amrex::MultiFab& Bz) override; - void copy_attribs(WarpXParIter& pti,const amrex::Real* xp, - const amrex::Real* yp, const amrex::Real* zp); + void PartitionParticlesInBuffers( + long& nfine_current, + long& nfine_gather, + long const np, + WarpXParIter& pti, + int const lev, + amrex::iMultiFab const* current_masks, + amrex::iMultiFab const* gather_masks, + RealVector& uxp, + RealVector& uyp, + RealVector& uzp, + RealVector& wp ); + + void copy_attribs(WarpXParIter& pti,const amrex::ParticleReal* xp, + const amrex::ParticleReal* yp, const amrex::ParticleReal* zp); virtual void PostRestart () final {} @@ -131,6 +145,33 @@ public: virtual void ConvertUnits (ConvertDirection convert_dir) override; +/** + * \brief Apply NCI Godfrey filter to all components of E and B before gather + * \param lev MR level + * \param box box onto which the filter is applied + * \param exeli safeguard to avoid destructing arrays between ParIter iterations on GPU + * \param filtered_Ex Array containing filtered value + * \param Ex Field array before filtering (not modified) + * \param ex_ptr pointer to the array used for field gather. + * + * The NCI Godfrey filter is applied on Ex, the result is stored in filtered_Ex + * and the pointer is modified (before this function is called, it points to Ex + * and after this function is called, it points to Ex_filtered) + */ + void applyNCIFilter ( + int lev, const amrex::Box& box, + amrex::Elixir& exeli, amrex::Elixir& eyeli, amrex::Elixir& ezeli, + amrex::Elixir& bxeli, amrex::Elixir& byeli, amrex::Elixir& bzeli, + amrex::FArrayBox& filtered_Ex, amrex::FArrayBox& filtered_Ey, + amrex::FArrayBox& filtered_Ez, amrex::FArrayBox& filtered_Bx, + amrex::FArrayBox& filtered_By, amrex::FArrayBox& filtered_Bz, + const amrex::FArrayBox& Ex, const amrex::FArrayBox& Ey, + const amrex::FArrayBox& Ez, const amrex::FArrayBox& Bx, + const amrex::FArrayBox& By, const amrex::FArrayBox& Bz, + amrex::FArrayBox const * & exfab, amrex::FArrayBox const * & eyfab, + amrex::FArrayBox const * & ezfab, amrex::FArrayBox const * & bxfab, + amrex::FArrayBox const * & byfab, amrex::FArrayBox const * & bzfab); + protected: std::string species_name; diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index b1e83d652..02dee1967 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -31,8 +31,11 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp pp.query("boost_adjust_transverse_positions", boost_adjust_transverse_positions); pp.query("do_backward_propagation", do_backward_propagation); + + // Initialize splitting pp.query("do_splitting", do_splitting); pp.query("split_type", split_type); + pp.query("do_continuous_injection", do_continuous_injection); // Whether to plot back-transformed (lab-frame) diagnostics // for this species. @@ -196,7 +199,7 @@ PhysicalParticleContainer::CheckAndAddParticle(Real x, Real y, Real z, std::array<Real, 3> u, Real weight) { - std::array<Real,PIdx::nattribs> attribs; + std::array<ParticleReal,PIdx::nattribs> attribs; attribs.fill(0.0); // update attribs with input arguments @@ -361,13 +364,13 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) for (int dir=0; dir<AMREX_SPACEDIM; dir++) { if ( tile_realbox.lo(dir) <= part_realbox.hi(dir) ) { Real ncells_adjust = std::floor( (tile_realbox.lo(dir) - part_realbox.lo(dir))/dx[dir] ); - overlap_realbox.setLo( dir, part_realbox.lo(dir) + std::max(ncells_adjust, 0.) * dx[dir]); + overlap_realbox.setLo( dir, part_realbox.lo(dir) + std::max(ncells_adjust, Real(0.)) * dx[dir]); } else { no_overlap = true; break; } if ( tile_realbox.hi(dir) >= part_realbox.lo(dir) ) { Real ncells_adjust = std::floor( (part_realbox.hi(dir) - tile_realbox.hi(dir))/dx[dir] ); - overlap_realbox.setHi( dir, part_realbox.hi(dir) - std::max(ncells_adjust, 0.) * dx[dir]); + overlap_realbox.setHi( dir, part_realbox.hi(dir) - std::max(ncells_adjust, Real(0.)) * dx[dir]); } else { no_overlap = true; break; } @@ -437,7 +440,7 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) ParticleType* pp = particle_tile.GetArrayOfStructs()().data() + old_size; auto& soa = particle_tile.GetStructOfArrays(); - GpuArray<Real*,PIdx::nattribs> pa; + GpuArray<ParticleReal*,PIdx::nattribs> pa; for (int ia = 0; ia < PIdx::nattribs; ++ia) { pa[ia] = soa.GetRealData(ia).data() + old_size; } @@ -948,7 +951,6 @@ PhysicalParticleContainer::Evolve (int lev, BL_PROFILE_VAR_NS("PPC::Evolve::Copy", blp_copy); BL_PROFILE_VAR_NS("PPC::FieldGather", blp_fg); BL_PROFILE_VAR_NS("PPC::ParticlePush", blp_ppc_pp); - BL_PROFILE_VAR_NS("PPC::Evolve::partition", blp_partition); const std::array<Real,3>& dx = WarpX::CellSize(lev); const std::array<Real,3>& cdx = WarpX::CellSize(std::max(lev-1,0)); @@ -960,7 +962,6 @@ PhysicalParticleContainer::Evolve (int lev, BL_ASSERT(OnSameGrids(lev,jx)); MultiFab* cost = WarpX::getCosts(lev); - const iMultiFab* current_masks = WarpX::CurrentBufferMasks(lev); const iMultiFab* gather_masks = WarpX::GatherBufferMasks(lev); @@ -991,10 +992,6 @@ PhysicalParticleContainer::Evolve (int lev, FArrayBox filtered_Ex, filtered_Ey, filtered_Ez; FArrayBox filtered_Bx, filtered_By, filtered_Bz; - std::vector<bool> inexflag; - Vector<long> pid; - RealVector tmp; - ParticleVector particle_tmp; for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) { @@ -1026,56 +1023,18 @@ PhysicalParticleContainer::Evolve (int lev, FArrayBox const* bzfab = &(Bz[pti]); Elixir exeli, eyeli, ezeli, bxeli, byeli, bzeli; + if (WarpX::use_fdtd_nci_corr) { -#if (AMREX_SPACEDIM == 2) - const Box& tbox = amrex::grow(pti.tilebox(),{static_cast<int>(WarpX::nox), - static_cast<int>(WarpX::noz)}); -#else - const Box& tbox = amrex::grow(pti.tilebox(),{static_cast<int>(WarpX::nox), - static_cast<int>(WarpX::noy), - static_cast<int>(WarpX::noz)}); -#endif - - // Filter Ex (Both 2D and 3D) - filtered_Ex.resize(amrex::convert(tbox,WarpX::Ex_nodal_flag)); - // Safeguard for GPU - exeli = filtered_Ex.elixir(); - // Apply filter on Ex, result stored in filtered_Ex - nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Ex, Ex[pti], filtered_Ex.box()); - // Update exfab reference - exfab = &filtered_Ex; - - // Filter Ez - filtered_Ez.resize(amrex::convert(tbox,WarpX::Ez_nodal_flag)); - ezeli = filtered_Ez.elixir(); - nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_Ez, Ez[pti], filtered_Ez.box()); - ezfab = &filtered_Ez; - - // Filter By - filtered_By.resize(amrex::convert(tbox,WarpX::By_nodal_flag)); - byeli = filtered_By.elixir(); - nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_By, By[pti], filtered_By.box()); - byfab = &filtered_By; -#if (AMREX_SPACEDIM == 3) - // Filter Ey - filtered_Ey.resize(amrex::convert(tbox,WarpX::Ey_nodal_flag)); - eyeli = filtered_Ey.elixir(); - nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Ey, Ey[pti], filtered_Ey.box()); - eyfab = &filtered_Ey; - - // Filter Bx - filtered_Bx.resize(amrex::convert(tbox,WarpX::Bx_nodal_flag)); - bxeli = filtered_Bx.elixir(); - nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_Bx, Bx[pti], filtered_Bx.box()); - bxfab = &filtered_Bx; - - // Filter Bz - filtered_Bz.resize(amrex::convert(tbox,WarpX::Bz_nodal_flag)); - bzeli = filtered_Bz.elixir(); - nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Bz, Bz[pti], filtered_Bz.box()); - bzfab = &filtered_Bz; -#endif + // Filter arrays Ex[pti], store the result in + // filtered_Ex and update pointer exfab so that it + // points to filtered_Ex (and do the same for all + // components of E and B). + applyNCIFilter(lev, pti.tilebox(), exeli, eyeli, ezeli, bxeli, byeli, bzeli, + filtered_Ex, filtered_Ey, filtered_Ez, + filtered_Bx, filtered_By, filtered_Bz, + Ex[pti], Ey[pti], Ez[pti], Bx[pti], By[pti], Bz[pti], + exfab, eyfab, ezfab, bxfab, byfab, bzfab); } Exp.assign(np,0.0); @@ -1085,99 +1044,21 @@ PhysicalParticleContainer::Evolve (int lev, Byp.assign(np,WarpX::B_external[1]); Bzp.assign(np,WarpX::B_external[2]); - long nfine_current = np; //! number of particles depositing to fine grid - long nfine_gather = np; //! number of particles gathering from fine grid - if (has_buffer && !do_not_push) - { - BL_PROFILE_VAR_START(blp_partition); - inexflag.resize(np); - auto& aos = pti.GetArrayOfStructs(); - // We need to partition the large buffer first - iMultiFab const* bmasks = (WarpX::n_field_gather_buffer >= WarpX::n_current_deposition_buffer) ? - gather_masks : current_masks; - int i = 0; - const auto& msk = (*bmasks)[pti]; - for (const auto& p : aos) { - const IntVect& iv = Index(p, lev); - inexflag[i++] = msk(iv); - } - - pid.resize(np); - std::iota(pid.begin(), pid.end(), 0L); - - auto sep = std::stable_partition(pid.begin(), pid.end(), - [&inexflag](long id) { return inexflag[id]; }); - - if (WarpX::n_current_deposition_buffer == WarpX::n_field_gather_buffer) { - nfine_current = nfine_gather = std::distance(pid.begin(), sep); - } else if (sep != pid.end()) { - int n_buf; - if (bmasks == gather_masks) { - nfine_gather = std::distance(pid.begin(), sep); - bmasks = current_masks; - n_buf = WarpX::n_current_deposition_buffer; - } else { - nfine_current = std::distance(pid.begin(), sep); - bmasks = gather_masks; - n_buf = WarpX::n_field_gather_buffer; - } - if (n_buf > 0) - { - const auto& msk2 = (*bmasks)[pti]; - for (auto it = sep; it != pid.end(); ++it) { - const long id = *it; - const IntVect& iv = Index(aos[id], lev); - inexflag[id] = msk2(iv); - } - - auto sep2 = std::stable_partition(sep, pid.end(), - [&inexflag](long id) {return inexflag[id];}); - if (bmasks == gather_masks) { - nfine_gather = std::distance(pid.begin(), sep2); - } else { - nfine_current = std::distance(pid.begin(), sep2); - } - } - } - - // only deposit / gather to coarsest grid - if (m_deposit_on_main_grid && lev > 0) { - nfine_current = 0; - } - if (m_gather_from_main_grid && lev > 0) { - nfine_gather = 0; - } - - if (nfine_current != np || nfine_gather != np) - { - particle_tmp.resize(np); - for (long ip = 0; ip < np; ++ip) { - particle_tmp[ip] = aos[pid[ip]]; - } - std::swap(aos(), particle_tmp); - - tmp.resize(np); - for (long ip = 0; ip < np; ++ip) { - tmp[ip] = wp[pid[ip]]; - } - std::swap(wp, tmp); - - for (long ip = 0; ip < np; ++ip) { - tmp[ip] = uxp[pid[ip]]; - } - std::swap(uxp, tmp); - - for (long ip = 0; ip < np; ++ip) { - tmp[ip] = uyp[pid[ip]]; - } - std::swap(uyp, tmp); - - for (long ip = 0; ip < np; ++ip) { - tmp[ip] = uzp[pid[ip]]; - } - std::swap(uzp, tmp); - } - BL_PROFILE_VAR_STOP(blp_partition); + // Determine which particles deposit/gather in the buffer, and + // which particles deposit/gather in the fine patch + long nfine_current = np; + long nfine_gather = np; + if (has_buffer && !do_not_push) { + // - Modify `nfine_current` and `nfine_gather` (in place) + // so that they correspond to the number of particles + // that deposit/gather in the fine patch respectively. + // - Reorder the particle arrays, + // so that the `nfine_current`/`nfine_gather` first particles + // deposit/gather in the fine patch + // and (thus) the `np-nfine_current`/`np-nfine_gather` last particles + // deposit/gather in the buffer + PartitionParticlesInBuffers( nfine_current, nfine_gather, np, + pti, lev, current_masks, gather_masks, uxp, uyp, uzp, wp ); } const long np_current = (cjx) ? nfine_current : np; @@ -1235,54 +1116,16 @@ PhysicalParticleContainer::Evolve (int lev, if (WarpX::use_fdtd_nci_corr) { -#if (AMREX_SPACEDIM == 2) - const Box& tbox = amrex::grow(cbox,{static_cast<int>(WarpX::nox), - static_cast<int>(WarpX::noz)}); -#else - const Box& tbox = amrex::grow(cbox,{static_cast<int>(WarpX::nox), - static_cast<int>(WarpX::noy), - static_cast<int>(WarpX::noz)}); -#endif - - // Filter Ex (both 2D and 3D) - filtered_Ex.resize(amrex::convert(tbox,WarpX::Ex_nodal_flag)); - // Safeguard for GPU - exeli = filtered_Ex.elixir(); - // Apply filter on Ex, result stored in filtered_Ex - nci_godfrey_filter_exeybz[lev-1]->ApplyStencil(filtered_Ex, (*cEx)[pti], filtered_Ex.box()); - // Update exfab reference - cexfab = &filtered_Ex; - - // Filter Ez - filtered_Ez.resize(amrex::convert(tbox,WarpX::Ez_nodal_flag)); - ezeli = filtered_Ez.elixir(); - nci_godfrey_filter_bxbyez[lev-1]->ApplyStencil(filtered_Ez, (*cEz)[pti], filtered_Ez.box()); - cezfab = &filtered_Ez; - - // Filter By - filtered_By.resize(amrex::convert(tbox,WarpX::By_nodal_flag)); - byeli = filtered_By.elixir(); - nci_godfrey_filter_bxbyez[lev-1]->ApplyStencil(filtered_By, (*cBy)[pti], filtered_By.box()); - cbyfab = &filtered_By; -#if (AMREX_SPACEDIM == 3) - // Filter Ey - filtered_Ey.resize(amrex::convert(tbox,WarpX::Ey_nodal_flag)); - eyeli = filtered_Ey.elixir(); - nci_godfrey_filter_exeybz[lev-1]->ApplyStencil(filtered_Ey, (*cEy)[pti], filtered_Ey.box()); - ceyfab = &filtered_Ey; - - // Filter Bx - filtered_Bx.resize(amrex::convert(tbox,WarpX::Bx_nodal_flag)); - bxeli = filtered_Bx.elixir(); - nci_godfrey_filter_bxbyez[lev-1]->ApplyStencil(filtered_Bx, (*cBx)[pti], filtered_Bx.box()); - cbxfab = &filtered_Bx; - - // Filter Bz - filtered_Bz.resize(amrex::convert(tbox,WarpX::Bz_nodal_flag)); - bzeli = filtered_Bz.elixir(); - nci_godfrey_filter_exeybz[lev-1]->ApplyStencil(filtered_Bz, (*cBz)[pti], filtered_Bz.box()); - cbzfab = &filtered_Bz; -#endif + // Filter arrays (*cEx)[pti], store the result in + // filtered_Ex and update pointer cexfab so that it + // points to filtered_Ex (and do the same for all + // components of E and B) + applyNCIFilter(lev-1, cbox, exeli, eyeli, ezeli, bxeli, byeli, bzeli, + filtered_Ex, filtered_Ey, filtered_Ez, + filtered_Bx, filtered_By, filtered_Bz, + (*cEx)[pti], (*cEy)[pti], (*cEz)[pti], + (*cBx)[pti], (*cBy)[pti], (*cBz)[pti], + cexfab, ceyfab, cezfab, cbxfab, cbyfab, cbzfab); } // Field gather for particles in gather buffers @@ -1375,6 +1218,74 @@ PhysicalParticleContainer::Evolve (int lev, } } +void +PhysicalParticleContainer::applyNCIFilter ( + int lev, const Box& box, + Elixir& exeli, Elixir& eyeli, Elixir& ezeli, + Elixir& bxeli, Elixir& byeli, Elixir& bzeli, + FArrayBox& filtered_Ex, FArrayBox& filtered_Ey, FArrayBox& filtered_Ez, + FArrayBox& filtered_Bx, FArrayBox& filtered_By, FArrayBox& filtered_Bz, + const FArrayBox& Ex, const FArrayBox& Ey, const FArrayBox& Ez, + const FArrayBox& Bx, const FArrayBox& By, const FArrayBox& Bz, + FArrayBox const * & ex_ptr, FArrayBox const * & ey_ptr, + FArrayBox const * & ez_ptr, FArrayBox const * & bx_ptr, + FArrayBox const * & by_ptr, FArrayBox const * & bz_ptr) +{ + + // Get instances of NCI Godfrey filters + const auto& nci_godfrey_filter_exeybz = WarpX::GetInstance().nci_godfrey_filter_exeybz; + const auto& nci_godfrey_filter_bxbyez = WarpX::GetInstance().nci_godfrey_filter_bxbyez; + +#if (AMREX_SPACEDIM == 2) + const Box& tbox = amrex::grow(box,{static_cast<int>(WarpX::nox), + static_cast<int>(WarpX::noz)}); +#else + const Box& tbox = amrex::grow(box,{static_cast<int>(WarpX::nox), + static_cast<int>(WarpX::noy), + static_cast<int>(WarpX::noz)}); +#endif + + // Filter Ex (Both 2D and 3D) + filtered_Ex.resize(amrex::convert(tbox,WarpX::Ex_nodal_flag)); + // Safeguard for GPU + exeli = filtered_Ex.elixir(); + // Apply filter on Ex, result stored in filtered_Ex + nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Ex, Ex, filtered_Ex.box()); + // Update ex_ptr reference + ex_ptr = &filtered_Ex; + + // Filter Ez + filtered_Ez.resize(amrex::convert(tbox,WarpX::Ez_nodal_flag)); + ezeli = filtered_Ez.elixir(); + nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_Ez, Ez, filtered_Ez.box()); + ez_ptr = &filtered_Ez; + + // Filter By + filtered_By.resize(amrex::convert(tbox,WarpX::By_nodal_flag)); + byeli = filtered_By.elixir(); + nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_By, By, filtered_By.box()); + by_ptr = &filtered_By; +#if (AMREX_SPACEDIM == 3) + // Filter Ey + filtered_Ey.resize(amrex::convert(tbox,WarpX::Ey_nodal_flag)); + eyeli = filtered_Ey.elixir(); + nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Ey, Ey, filtered_Ey.box()); + ey_ptr = &filtered_Ey; + + // Filter Bx + filtered_Bx.resize(amrex::convert(tbox,WarpX::Bx_nodal_flag)); + bxeli = filtered_Bx.elixir(); + nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_Bx, Bx, filtered_Bx.box()); + bx_ptr = &filtered_Bx; + + // Filter Bz + filtered_Bz.resize(amrex::convert(tbox,WarpX::Bz_nodal_flag)); + bzeli = filtered_Bz.elixir(); + nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Bz, Bz, filtered_Bz.box()); + bz_ptr = &filtered_Bz; +#endif +} + // Loop over all particles in the particle container and // split particles tagged with p.id()=DoSplitParticleID void @@ -1382,7 +1293,7 @@ PhysicalParticleContainer::SplitParticles(int lev) { auto& mypc = WarpX::GetInstance().GetPartContainer(); auto& pctmp_split = mypc.GetPCtmp(); - Cuda::ManagedDeviceVector<Real> xp, yp, zp; + Cuda::ManagedDeviceVector<ParticleReal> xp, yp, zp; RealVector psplit_x, psplit_y, psplit_z, psplit_w; RealVector psplit_ux, psplit_uy, psplit_uz; long np_split_to_add = 0; @@ -1398,7 +1309,16 @@ PhysicalParticleContainer::SplitParticles(int lev) for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) { pti.GetPosition(xp, yp, zp); + + // offset for split particles is computed as a function of cell size + // and number of particles per cell, so that a uniform distribution + // before splitting results in a uniform distribution after splitting + const amrex::Vector<int> ppc_nd = plasma_injector->num_particles_per_cell_each_dim; const std::array<Real,3>& dx = WarpX::CellSize(lev); + amrex::Vector<amrex::Real> split_offset = {dx[0]/2./ppc_nd[0], + dx[1]/2./ppc_nd[1], + dx[2]/2./ppc_nd[2]}; + // particle Array Of Structs data auto& particles = pti.GetArrayOfStructs(); // particle Struct Of Arrays data @@ -1421,9 +1341,9 @@ PhysicalParticleContainer::SplitParticles(int lev) for (int ishift = -1; ishift < 2; ishift +=2 ){ for (int kshift = -1; kshift < 2; kshift +=2 ){ // Add one particle with offset in x and z - psplit_x.push_back( xp[i] + ishift*dx[0]/2 ); + psplit_x.push_back( xp[i] + ishift*split_offset[0] ); psplit_y.push_back( yp[i] ); - psplit_z.push_back( zp[i] + kshift*dx[2]/2 ); + psplit_z.push_back( zp[i] + kshift*split_offset[2] ); psplit_ux.push_back( uxp[i] ); psplit_uy.push_back( uyp[i] ); psplit_uz.push_back( uzp[i] ); @@ -1435,7 +1355,7 @@ PhysicalParticleContainer::SplitParticles(int lev) // 4 particles in 2d for (int ishift = -1; ishift < 2; ishift +=2 ){ // Add one particle with offset in x - psplit_x.push_back( xp[i] + ishift*dx[0]/2 ); + psplit_x.push_back( xp[i] + ishift*split_offset[0] ); psplit_y.push_back( yp[i] ); psplit_z.push_back( zp[i] ); psplit_ux.push_back( uxp[i] ); @@ -1445,7 +1365,7 @@ PhysicalParticleContainer::SplitParticles(int lev) // Add one particle with offset in z psplit_x.push_back( xp[i] ); psplit_y.push_back( yp[i] ); - psplit_z.push_back( zp[i] + ishift*dx[2]/2 ); + psplit_z.push_back( zp[i] + ishift*split_offset[2] ); psplit_ux.push_back( uxp[i] ); psplit_uy.push_back( uyp[i] ); psplit_uz.push_back( uzp[i] ); @@ -1460,9 +1380,9 @@ PhysicalParticleContainer::SplitParticles(int lev) for (int jshift = -1; jshift < 2; jshift +=2 ){ for (int kshift = -1; kshift < 2; kshift +=2 ){ // Add one particle with offset in x, y and z - psplit_x.push_back( xp[i] + ishift*dx[0]/2 ); - psplit_y.push_back( yp[i] + jshift*dx[1]/2 ); - psplit_z.push_back( zp[i] + kshift*dx[2]/2 ); + psplit_x.push_back( xp[i] + ishift*split_offset[0] ); + psplit_y.push_back( yp[i] + jshift*split_offset[1] ); + psplit_z.push_back( zp[i] + kshift*split_offset[2] ); psplit_ux.push_back( uxp[i] ); psplit_uy.push_back( uyp[i] ); psplit_uz.push_back( uzp[i] ); @@ -1475,7 +1395,7 @@ PhysicalParticleContainer::SplitParticles(int lev) // 6 particles in 3d for (int ishift = -1; ishift < 2; ishift +=2 ){ // Add one particle with offset in x - psplit_x.push_back( xp[i] + ishift*dx[0]/2 ); + psplit_x.push_back( xp[i] + ishift*split_offset[0] ); psplit_y.push_back( yp[i] ); psplit_z.push_back( zp[i] ); psplit_ux.push_back( uxp[i] ); @@ -1484,7 +1404,7 @@ PhysicalParticleContainer::SplitParticles(int lev) psplit_w.push_back( wp[i]/np_split ); // Add one particle with offset in y psplit_x.push_back( xp[i] ); - psplit_y.push_back( yp[i] + ishift*dx[1]/2 ); + psplit_y.push_back( yp[i] + ishift*split_offset[1] ); psplit_z.push_back( zp[i] ); psplit_ux.push_back( uxp[i] ); psplit_uy.push_back( uyp[i] ); @@ -1493,7 +1413,7 @@ PhysicalParticleContainer::SplitParticles(int lev) // Add one particle with offset in z psplit_x.push_back( xp[i] ); psplit_y.push_back( yp[i] ); - psplit_z.push_back( zp[i] + ishift*dx[2]/2 ); + psplit_z.push_back( zp[i] + ishift*split_offset[2] ); psplit_ux.push_back( uxp[i] ); psplit_uy.push_back( uyp[i] ); psplit_uz.push_back( uzp[i] ); @@ -1531,27 +1451,27 @@ PhysicalParticleContainer::SplitParticles(int lev) void PhysicalParticleContainer::PushPX(WarpXParIter& pti, - Cuda::ManagedDeviceVector<Real>& xp, - Cuda::ManagedDeviceVector<Real>& yp, - Cuda::ManagedDeviceVector<Real>& zp, + Cuda::ManagedDeviceVector<ParticleReal>& xp, + Cuda::ManagedDeviceVector<ParticleReal>& yp, + Cuda::ManagedDeviceVector<ParticleReal>& zp, Real dt, DtType a_dt_type) { // This wraps the momentum and position advance so that inheritors can modify the call. auto& attribs = pti.GetAttribs(); // Extract pointers to the different particle quantities - Real* const AMREX_RESTRICT x = xp.dataPtr(); - Real* const AMREX_RESTRICT y = yp.dataPtr(); - Real* const AMREX_RESTRICT z = zp.dataPtr(); - Real* const AMREX_RESTRICT ux = attribs[PIdx::ux].dataPtr(); - Real* const AMREX_RESTRICT uy = attribs[PIdx::uy].dataPtr(); - Real* const AMREX_RESTRICT uz = attribs[PIdx::uz].dataPtr(); - const Real* const AMREX_RESTRICT Ex = attribs[PIdx::Ex].dataPtr(); - const Real* const AMREX_RESTRICT Ey = attribs[PIdx::Ey].dataPtr(); - const Real* const AMREX_RESTRICT Ez = attribs[PIdx::Ez].dataPtr(); - const Real* const AMREX_RESTRICT Bx = attribs[PIdx::Bx].dataPtr(); - const Real* const AMREX_RESTRICT By = attribs[PIdx::By].dataPtr(); - const Real* const AMREX_RESTRICT Bz = attribs[PIdx::Bz].dataPtr(); + ParticleReal* const AMREX_RESTRICT x = xp.dataPtr(); + ParticleReal* const AMREX_RESTRICT y = yp.dataPtr(); + ParticleReal* const AMREX_RESTRICT z = zp.dataPtr(); + ParticleReal* const AMREX_RESTRICT ux = attribs[PIdx::ux].dataPtr(); + ParticleReal* const AMREX_RESTRICT uy = attribs[PIdx::uy].dataPtr(); + ParticleReal* const AMREX_RESTRICT uz = attribs[PIdx::uz].dataPtr(); + const ParticleReal* const AMREX_RESTRICT Ex = attribs[PIdx::Ex].dataPtr(); + const ParticleReal* const AMREX_RESTRICT Ey = attribs[PIdx::Ey].dataPtr(); + const ParticleReal* const AMREX_RESTRICT Ez = attribs[PIdx::Ez].dataPtr(); + const ParticleReal* const AMREX_RESTRICT Bx = attribs[PIdx::Bx].dataPtr(); + const ParticleReal* const AMREX_RESTRICT By = attribs[PIdx::By].dataPtr(); + const ParticleReal* const AMREX_RESTRICT Bz = attribs[PIdx::Bz].dataPtr(); if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags && (a_dt_type!=DtType::SecondHalf)) { @@ -1660,15 +1580,15 @@ PhysicalParticleContainer::PushP (int lev, Real dt, // This wraps the momentum advance so that inheritors can modify the call. // Extract pointers to the different particle quantities - Real* const AMREX_RESTRICT ux = attribs[PIdx::ux].dataPtr(); - Real* const AMREX_RESTRICT uy = attribs[PIdx::uy].dataPtr(); - Real* const AMREX_RESTRICT uz = attribs[PIdx::uz].dataPtr(); - const Real* const AMREX_RESTRICT Expp = Exp.dataPtr(); - const Real* const AMREX_RESTRICT Eypp = Eyp.dataPtr(); - const Real* const AMREX_RESTRICT Ezpp = Ezp.dataPtr(); - const Real* const AMREX_RESTRICT Bxpp = Bxp.dataPtr(); - const Real* const AMREX_RESTRICT Bypp = Byp.dataPtr(); - const Real* const AMREX_RESTRICT Bzpp = Bzp.dataPtr(); + ParticleReal* const AMREX_RESTRICT ux = attribs[PIdx::ux].dataPtr(); + ParticleReal* const AMREX_RESTRICT uy = attribs[PIdx::uy].dataPtr(); + ParticleReal* const AMREX_RESTRICT uz = attribs[PIdx::uz].dataPtr(); + const ParticleReal* const AMREX_RESTRICT Expp = Exp.dataPtr(); + const ParticleReal* const AMREX_RESTRICT Eypp = Eyp.dataPtr(); + const ParticleReal* const AMREX_RESTRICT Ezpp = Ezp.dataPtr(); + const ParticleReal* const AMREX_RESTRICT Bxpp = Bxp.dataPtr(); + const ParticleReal* const AMREX_RESTRICT Bypp = Byp.dataPtr(); + const ParticleReal* const AMREX_RESTRICT Bzpp = Bzp.dataPtr(); // Loop over the particles and update their momentum const Real q = this->charge; @@ -1694,23 +1614,23 @@ PhysicalParticleContainer::PushP (int lev, Real dt, } } -void PhysicalParticleContainer::copy_attribs(WarpXParIter& pti,const Real* xp, - const Real* yp, const Real* zp) +void PhysicalParticleContainer::copy_attribs(WarpXParIter& pti,const ParticleReal* xp, + const ParticleReal* yp, const ParticleReal* zp) { auto& attribs = pti.GetAttribs(); - Real* AMREX_RESTRICT uxp = attribs[PIdx::ux].dataPtr(); - Real* AMREX_RESTRICT uyp = attribs[PIdx::uy].dataPtr(); - Real* AMREX_RESTRICT uzp = attribs[PIdx::uz].dataPtr(); + ParticleReal* AMREX_RESTRICT uxp = attribs[PIdx::ux].dataPtr(); + ParticleReal* AMREX_RESTRICT uyp = attribs[PIdx::uy].dataPtr(); + ParticleReal* AMREX_RESTRICT uzp = attribs[PIdx::uz].dataPtr(); const auto np = pti.numParticles(); const auto lev = pti.GetLevel(); const auto index = pti.GetPairIndex(); - Real* AMREX_RESTRICT xpold = tmp_particle_data[lev][index][TmpIdx::xold ].dataPtr(); - Real* AMREX_RESTRICT ypold = tmp_particle_data[lev][index][TmpIdx::yold ].dataPtr(); - Real* AMREX_RESTRICT zpold = tmp_particle_data[lev][index][TmpIdx::zold ].dataPtr(); - Real* AMREX_RESTRICT uxpold = tmp_particle_data[lev][index][TmpIdx::uxold].dataPtr(); - Real* AMREX_RESTRICT uypold = tmp_particle_data[lev][index][TmpIdx::uyold].dataPtr(); - Real* AMREX_RESTRICT uzpold = tmp_particle_data[lev][index][TmpIdx::uzold].dataPtr(); + ParticleReal* AMREX_RESTRICT xpold = tmp_particle_data[lev][index][TmpIdx::xold ].dataPtr(); + ParticleReal* AMREX_RESTRICT ypold = tmp_particle_data[lev][index][TmpIdx::yold ].dataPtr(); + ParticleReal* AMREX_RESTRICT zpold = tmp_particle_data[lev][index][TmpIdx::zold ].dataPtr(); + ParticleReal* AMREX_RESTRICT uxpold = tmp_particle_data[lev][index][TmpIdx::uxold].dataPtr(); + ParticleReal* AMREX_RESTRICT uypold = tmp_particle_data[lev][index][TmpIdx::uyold].dataPtr(); + ParticleReal* AMREX_RESTRICT uzpold = tmp_particle_data[lev][index][TmpIdx::uzold].dataPtr(); ParallelFor( np, [=] AMREX_GPU_DEVICE (long i) { @@ -1929,9 +1849,9 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti, const Array4<const Real>& by_arr = byfab->array(); const Array4<const Real>& bz_arr = bzfab->array(); - const Real * const AMREX_RESTRICT xp = m_xp[thread_num].dataPtr() + offset; - const Real * const AMREX_RESTRICT zp = m_zp[thread_num].dataPtr() + offset; - const Real * const AMREX_RESTRICT yp = m_yp[thread_num].dataPtr() + offset; + const ParticleReal * const AMREX_RESTRICT xp = m_xp[thread_num].dataPtr() + offset; + const ParticleReal * const AMREX_RESTRICT zp = m_zp[thread_num].dataPtr() + offset; + const ParticleReal * const AMREX_RESTRICT yp = m_yp[thread_num].dataPtr() + offset; // Lower corner of tile box physical domain const std::array<Real, 3>& xyzmin = WarpX::LowerCorner(box, gather_lev); @@ -2078,15 +1998,15 @@ PhysicalParticleContainer::buildIonizationMask (const amrex::MFIter& mfi, const // Otherwise, resize ionization_mask, and get poiters to attribs arrays. ionization_mask.resize(np); int * const AMREX_RESTRICT p_ionization_mask = ionization_mask.data(); - const Real * const AMREX_RESTRICT ux = soa.GetRealData(PIdx::ux).data(); - const Real * const AMREX_RESTRICT uy = soa.GetRealData(PIdx::uy).data(); - const Real * const AMREX_RESTRICT uz = soa.GetRealData(PIdx::uz).data(); - const Real * const AMREX_RESTRICT ex = soa.GetRealData(PIdx::Ex).data(); - const Real * const AMREX_RESTRICT ey = soa.GetRealData(PIdx::Ey).data(); - const Real * const AMREX_RESTRICT ez = soa.GetRealData(PIdx::Ez).data(); - const Real * const AMREX_RESTRICT bx = soa.GetRealData(PIdx::Bx).data(); - const Real * const AMREX_RESTRICT by = soa.GetRealData(PIdx::By).data(); - const Real * const AMREX_RESTRICT bz = soa.GetRealData(PIdx::Bz).data(); + const ParticleReal * const AMREX_RESTRICT ux = soa.GetRealData(PIdx::ux).data(); + const ParticleReal * const AMREX_RESTRICT uy = soa.GetRealData(PIdx::uy).data(); + const ParticleReal * const AMREX_RESTRICT uz = soa.GetRealData(PIdx::uz).data(); + const ParticleReal * const AMREX_RESTRICT ex = soa.GetRealData(PIdx::Ex).data(); + const ParticleReal * const AMREX_RESTRICT ey = soa.GetRealData(PIdx::Ey).data(); + const ParticleReal * const AMREX_RESTRICT ez = soa.GetRealData(PIdx::Ez).data(); + const ParticleReal * const AMREX_RESTRICT bx = soa.GetRealData(PIdx::Bx).data(); + const ParticleReal * const AMREX_RESTRICT by = soa.GetRealData(PIdx::By).data(); + const ParticleReal * const AMREX_RESTRICT bz = soa.GetRealData(PIdx::Bz).data(); int* ion_lev = soa.GetIntData(particle_icomps["ionization_level"]).data(); Real c = PhysConst::c; diff --git a/Source/Particles/Pusher/GetAndSetPosition.H b/Source/Particles/Pusher/GetAndSetPosition.H index 3c74baeb2..f0dfa4c83 100644 --- a/Source/Particles/Pusher/GetAndSetPosition.H +++ b/Source/Particles/Pusher/GetAndSetPosition.H @@ -11,7 +11,7 @@ * and stores them in the variables `x`, `y`, `z`. */ AMREX_GPU_HOST_DEVICE AMREX_INLINE void GetPosition( - amrex::Real& x, amrex::Real& y, amrex::Real& z, + amrex::ParticleReal& x, amrex::ParticleReal& y, amrex::ParticleReal& z, const WarpXParticleContainer::ParticleType& p) { #if (AMREX_SPACEDIM==3) @@ -20,7 +20,7 @@ void GetPosition( z = p.pos(2); #else x = p.pos(0); - y = std::numeric_limits<amrex::Real>::quiet_NaN(); + y = std::numeric_limits<amrex::ParticleReal>::quiet_NaN(); z = p.pos(1); #endif } @@ -30,7 +30,7 @@ void GetPosition( AMREX_GPU_HOST_DEVICE AMREX_INLINE void SetPosition( WarpXParticleContainer::ParticleType& p, - const amrex::Real x, const amrex::Real y, const amrex::Real z) + const amrex::ParticleReal x, const amrex::ParticleReal y, const amrex::ParticleReal z) { #if (AMREX_SPACEDIM==3) p.pos(0) = x; @@ -49,10 +49,10 @@ void SetPosition( * and store them in the variables `x`, `y`, `z` */ AMREX_GPU_HOST_DEVICE AMREX_INLINE void GetCartesianPositionFromCylindrical( - amrex::Real& x, amrex::Real& y, amrex::Real& z, - const WarpXParticleContainer::ParticleType& p, const amrex::Real theta) + amrex::ParticleReal& x, amrex::ParticleReal& y, amrex::ParticleReal& z, + const WarpXParticleContainer::ParticleType& p, const amrex::ParticleReal theta) { - const amrex::Real r = p.pos(0); + const amrex::ParticleReal r = p.pos(0); x = r*std::cos(theta); y = r*std::sin(theta); z = p.pos(1); @@ -63,8 +63,8 @@ void GetCartesianPositionFromCylindrical( * from the values of `x`, `y`, `z` */ AMREX_GPU_HOST_DEVICE AMREX_INLINE void SetCylindricalPositionFromCartesian( - WarpXParticleContainer::ParticleType& p, amrex::Real& theta, - const amrex::Real x, const amrex::Real y, const amrex::Real z ) + WarpXParticleContainer::ParticleType& p, amrex::ParticleReal& theta, + const amrex::ParticleReal x, const amrex::ParticleReal y, const amrex::ParticleReal z ) { theta = std::atan2(y, x); p.pos(0) = std::sqrt(x*x + y*y); diff --git a/Source/Particles/Pusher/UpdateMomentumBoris.H b/Source/Particles/Pusher/UpdateMomentumBoris.H index a33058347..205cc9a71 100644 --- a/Source/Particles/Pusher/UpdateMomentumBoris.H +++ b/Source/Particles/Pusher/UpdateMomentumBoris.H @@ -7,9 +7,9 @@ * given the value of its momenta `ux`, `uy`, `uz` */ AMREX_GPU_HOST_DEVICE AMREX_INLINE void UpdateMomentumBoris( - amrex::Real& ux, amrex::Real& uy, amrex::Real& uz, - const amrex::Real Ex, const amrex::Real Ey, const amrex::Real Ez, - const amrex::Real Bx, const amrex::Real By, const amrex::Real Bz, + 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::Real econst = 0.5*q*dt/m; diff --git a/Source/Particles/Pusher/UpdateMomentumVay.H b/Source/Particles/Pusher/UpdateMomentumVay.H index 1f0f19e63..433a891c5 100644 --- a/Source/Particles/Pusher/UpdateMomentumVay.H +++ b/Source/Particles/Pusher/UpdateMomentumVay.H @@ -9,9 +9,9 @@ * given the value of its momenta `ux`, `uy`, `uz` */ AMREX_GPU_HOST_DEVICE AMREX_INLINE void UpdateMomentumVay( - amrex::Real& ux, amrex::Real& uy, amrex::Real& uz, - const amrex::Real Ex, const amrex::Real Ey, const amrex::Real Ez, - const amrex::Real Bx, const amrex::Real By, const amrex::Real Bz, + 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 ) { // Constants diff --git a/Source/Particles/Pusher/UpdatePosition.H b/Source/Particles/Pusher/UpdatePosition.H index a9df63a30..da0e9cdf9 100644 --- a/Source/Particles/Pusher/UpdatePosition.H +++ b/Source/Particles/Pusher/UpdatePosition.H @@ -9,8 +9,8 @@ * given the value of its momenta `ux`, `uy`, `uz` */ AMREX_GPU_HOST_DEVICE AMREX_INLINE void UpdatePosition( - amrex::Real& x, amrex::Real& y, amrex::Real& z, - const amrex::Real ux, const amrex::Real uy, const amrex::Real uz, + amrex::ParticleReal& x, amrex::ParticleReal& y, amrex::ParticleReal& z, + const amrex::ParticleReal ux, const amrex::ParticleReal uy, const amrex::ParticleReal uz, const amrex::Real dt ) { diff --git a/Source/Particles/Pusher/UpdatePositionPhoton.H b/Source/Particles/Pusher/UpdatePositionPhoton.H index bd6e6cd21..f95c2b09d 100644 --- a/Source/Particles/Pusher/UpdatePositionPhoton.H +++ b/Source/Particles/Pusher/UpdatePositionPhoton.H @@ -10,8 +10,8 @@ * given the value of its momenta `ux`, `uy`, `uz` */ AMREX_GPU_HOST_DEVICE AMREX_INLINE void UpdatePositionPhoton( - amrex::Real& x, amrex::Real& y, amrex::Real& z, - const amrex::Real ux, const amrex::Real uy, const amrex::Real uz, + amrex::ParticleReal& x, amrex::ParticleReal& y, amrex::ParticleReal& z, + const amrex::ParticleReal ux, const amrex::ParticleReal uy, const amrex::ParticleReal uz, const amrex::Real dt ) { // Compute speed of light over inverse of momentum modulus diff --git a/Source/Particles/RigidInjectedParticleContainer.H b/Source/Particles/RigidInjectedParticleContainer.H index be3dd21f9..3abbb4afe 100644 --- a/Source/Particles/RigidInjectedParticleContainer.H +++ b/Source/Particles/RigidInjectedParticleContainer.H @@ -44,9 +44,9 @@ public: DtType a_dt_type=DtType::Full) override; virtual void PushPX(WarpXParIter& pti, - amrex::Cuda::ManagedDeviceVector<amrex::Real>& xp, - amrex::Cuda::ManagedDeviceVector<amrex::Real>& yp, - amrex::Cuda::ManagedDeviceVector<amrex::Real>& zp, + amrex::Cuda::ManagedDeviceVector<amrex::ParticleReal>& xp, + amrex::Cuda::ManagedDeviceVector<amrex::ParticleReal>& yp, + amrex::Cuda::ManagedDeviceVector<amrex::ParticleReal>& zp, amrex::Real dt, DtType a_dt_type=DtType::Full) override; virtual void PushP (int lev, amrex::Real dt, diff --git a/Source/Particles/RigidInjectedParticleContainer.cpp b/Source/Particles/RigidInjectedParticleContainer.cpp index 7d129fc01..891ade76d 100644 --- a/Source/Particles/RigidInjectedParticleContainer.cpp +++ b/Source/Particles/RigidInjectedParticleContainer.cpp @@ -76,7 +76,7 @@ RigidInjectedParticleContainer::RemapParticles() // Note that the particles are already in the boosted frame. // This value is saved to advance the particles not injected yet - Cuda::ManagedDeviceVector<Real> xp, yp, zp; + Cuda::ManagedDeviceVector<ParticleReal> xp, yp, zp; for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) { @@ -136,7 +136,7 @@ RigidInjectedParticleContainer::BoostandRemapParticles() #pragma omp parallel #endif { - Cuda::ManagedDeviceVector<Real> xp, yp, zp; + Cuda::ManagedDeviceVector<ParticleReal> xp, yp, zp; for (WarpXParIter pti(*this, 0); pti.isValid(); ++pti) { @@ -207,9 +207,9 @@ RigidInjectedParticleContainer::BoostandRemapParticles() void RigidInjectedParticleContainer::PushPX(WarpXParIter& pti, - Cuda::ManagedDeviceVector<Real>& xp, - Cuda::ManagedDeviceVector<Real>& yp, - Cuda::ManagedDeviceVector<Real>& zp, + Cuda::ManagedDeviceVector<ParticleReal>& xp, + Cuda::ManagedDeviceVector<ParticleReal>& yp, + Cuda::ManagedDeviceVector<ParticleReal>& zp, Real dt, DtType a_dt_type) { @@ -220,21 +220,21 @@ RigidInjectedParticleContainer::PushPX(WarpXParIter& pti, auto& uzp = attribs[PIdx::uz]; // Save the position and momenta, making copies - Cuda::ManagedDeviceVector<Real> xp_save, yp_save, zp_save; + Cuda::ManagedDeviceVector<ParticleReal> xp_save, yp_save, zp_save; RealVector uxp_save, uyp_save, uzp_save; - Real* const AMREX_RESTRICT x = xp.dataPtr(); - Real* const AMREX_RESTRICT y = yp.dataPtr(); - Real* const AMREX_RESTRICT z = zp.dataPtr(); - Real* const AMREX_RESTRICT ux = uxp.dataPtr(); - Real* const AMREX_RESTRICT uy = uyp.dataPtr(); - Real* const AMREX_RESTRICT uz = uzp.dataPtr(); - Real* const AMREX_RESTRICT Exp = attribs[PIdx::Ex].dataPtr(); - Real* const AMREX_RESTRICT Eyp = attribs[PIdx::Ey].dataPtr(); - Real* const AMREX_RESTRICT Ezp = attribs[PIdx::Ez].dataPtr(); - Real* const AMREX_RESTRICT Bxp = attribs[PIdx::Bx].dataPtr(); - Real* const AMREX_RESTRICT Byp = attribs[PIdx::By].dataPtr(); - Real* const AMREX_RESTRICT Bzp = attribs[PIdx::Bz].dataPtr(); + ParticleReal* const AMREX_RESTRICT x = xp.dataPtr(); + ParticleReal* const AMREX_RESTRICT y = yp.dataPtr(); + ParticleReal* const AMREX_RESTRICT z = zp.dataPtr(); + ParticleReal* const AMREX_RESTRICT ux = uxp.dataPtr(); + ParticleReal* const AMREX_RESTRICT uy = uyp.dataPtr(); + ParticleReal* const AMREX_RESTRICT uz = uzp.dataPtr(); + ParticleReal* const AMREX_RESTRICT Exp = attribs[PIdx::Ex].dataPtr(); + ParticleReal* const AMREX_RESTRICT Eyp = attribs[PIdx::Ey].dataPtr(); + ParticleReal* const AMREX_RESTRICT Ezp = attribs[PIdx::Ez].dataPtr(); + ParticleReal* const AMREX_RESTRICT Bxp = attribs[PIdx::Bx].dataPtr(); + ParticleReal* const AMREX_RESTRICT Byp = attribs[PIdx::By].dataPtr(); + ParticleReal* const AMREX_RESTRICT Bzp = attribs[PIdx::Bz].dataPtr(); if (!done_injecting_lev) { // If the old values are not already saved, create copies here. @@ -271,12 +271,12 @@ RigidInjectedParticleContainer::PushPX(WarpXParIter& pti, if (!done_injecting_lev) { - Real* AMREX_RESTRICT x_save = xp_save.dataPtr(); - Real* AMREX_RESTRICT y_save = yp_save.dataPtr(); - Real* AMREX_RESTRICT z_save = zp_save.dataPtr(); - Real* AMREX_RESTRICT ux_save = uxp_save.dataPtr(); - Real* AMREX_RESTRICT uy_save = uyp_save.dataPtr(); - Real* AMREX_RESTRICT uz_save = uzp_save.dataPtr(); + ParticleReal* AMREX_RESTRICT x_save = xp_save.dataPtr(); + ParticleReal* AMREX_RESTRICT y_save = yp_save.dataPtr(); + ParticleReal* AMREX_RESTRICT z_save = zp_save.dataPtr(); + ParticleReal* AMREX_RESTRICT ux_save = uxp_save.dataPtr(); + ParticleReal* AMREX_RESTRICT uy_save = uyp_save.dataPtr(); + ParticleReal* AMREX_RESTRICT uz_save = uzp_save.dataPtr(); // Undo the push for particles not injected yet. // The zp are advanced a fixed amount. @@ -415,16 +415,16 @@ RigidInjectedParticleContainer::PushP (int lev, Real dt, // This wraps the momentum advance so that inheritors can modify the call. // Extract pointers to the different particle quantities - const Real* const AMREX_RESTRICT zp = m_zp[thread_num].dataPtr(); - Real* const AMREX_RESTRICT uxpp = uxp.dataPtr(); - Real* const AMREX_RESTRICT uypp = uyp.dataPtr(); - Real* const AMREX_RESTRICT uzpp = uzp.dataPtr(); - const Real* const AMREX_RESTRICT Expp = Exp.dataPtr(); - const Real* const AMREX_RESTRICT Eypp = Eyp.dataPtr(); - const Real* const AMREX_RESTRICT Ezpp = Ezp.dataPtr(); - const Real* const AMREX_RESTRICT Bxpp = Bxp.dataPtr(); - const Real* const AMREX_RESTRICT Bypp = Byp.dataPtr(); - const Real* const AMREX_RESTRICT Bzpp = Bzp.dataPtr(); + const ParticleReal* const AMREX_RESTRICT zp = m_zp[thread_num].dataPtr(); + ParticleReal* const AMREX_RESTRICT uxpp = uxp.dataPtr(); + ParticleReal* const AMREX_RESTRICT uypp = uyp.dataPtr(); + ParticleReal* const AMREX_RESTRICT uzpp = uzp.dataPtr(); + const ParticleReal* const AMREX_RESTRICT Expp = Exp.dataPtr(); + const ParticleReal* const AMREX_RESTRICT Eypp = Eyp.dataPtr(); + const ParticleReal* const AMREX_RESTRICT Ezpp = Ezp.dataPtr(); + const ParticleReal* const AMREX_RESTRICT Bxpp = Bxp.dataPtr(); + const ParticleReal* const AMREX_RESTRICT Bypp = Byp.dataPtr(); + const ParticleReal* const AMREX_RESTRICT Bzpp = Bzp.dataPtr(); // Loop over the particles and update their momentum const Real q = this->charge; @@ -450,10 +450,10 @@ RigidInjectedParticleContainer::PushP (int lev, Real dt, // Undo the push for particles not injected yet. // It is assumed that PushP will only be called on the first and last steps // and that no particles will cross zinject_plane. - const Real* const AMREX_RESTRICT ux_save = uxp_save.dataPtr(); - const Real* const AMREX_RESTRICT uy_save = uyp_save.dataPtr(); - const Real* const AMREX_RESTRICT uz_save = uzp_save.dataPtr(); - const Real zz = zinject_plane_levels[lev]; + const ParticleReal* const AMREX_RESTRICT ux_save = uxp_save.dataPtr(); + const ParticleReal* const AMREX_RESTRICT uy_save = uyp_save.dataPtr(); + const ParticleReal* const AMREX_RESTRICT uz_save = uzp_save.dataPtr(); + const ParticleReal zz = zinject_plane_levels[lev]; amrex::ParallelFor( pti.numParticles(), [=] AMREX_GPU_DEVICE (long i) { if (zp[i] <= zz) { diff --git a/Source/Particles/Sorting/Make.package b/Source/Particles/Sorting/Make.package new file mode 100644 index 000000000..750d2607e --- /dev/null +++ b/Source/Particles/Sorting/Make.package @@ -0,0 +1,3 @@ +CEXE_sources += Partition.cpp +INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Particles/Sorting +VPATH_LOCATIONS += $(WARPX_HOME)/Source/Particles/Sorting diff --git a/Source/Particles/Sorting/Partition.cpp b/Source/Particles/Sorting/Partition.cpp new file mode 100644 index 000000000..683dbbd04 --- /dev/null +++ b/Source/Particles/Sorting/Partition.cpp @@ -0,0 +1,141 @@ +#include <WarpX.H> +#include <PhysicalParticleContainer.H> +#include <AMReX_Particles.H> + +using namespace amrex; + +/* \brief Determine which particles deposit/gather in the buffer, and + * and reorder the particle arrays accordingly + * + * More specifically: + * - Modify `nfine_current` and `nfine_gather` (in place) + * so that they correspond to the number of particles + * that deposit/gather in the fine patch respectively. + * - Reorder the particle arrays, + * so that the `nfine_current`/`nfine_gather` first particles + * deposit/gather in the fine patch + * and (thus) the `np-nfine_current`/`np-nfine_gather` last particles + * deposit/gather in the buffer + * + * \param nfine_current number of particles that deposit to the fine patch + * (modified by this function) + * \param nfine_gather number of particles that gather into the fine patch + * (modified by this function) + * \param np total number of particles in this tile + * \param pti object that holds the particle information for this tile + * \param lev current refinement level + * \param current_masks indicates, for each cell, whether that cell is + * in the deposition buffers or in the interior of the fine patch + * \param gather_masks indicates, for each cell, whether that cell is + * in the gather buffers or in the interior of the fine patch + * \param uxp, uyp, uzp, wp references to the particle momenta and weight + * (modified by this function) + */ +void +PhysicalParticleContainer::PartitionParticlesInBuffers( + long& nfine_current, long& nfine_gather, long const np, + WarpXParIter& pti, int const lev, + iMultiFab const* current_masks, + iMultiFab const* gather_masks, + RealVector& uxp, RealVector& uyp, RealVector& uzp, RealVector& wp) +{ + BL_PROFILE("PPC::Evolve::partition"); + + std::vector<bool> inexflag; + inexflag.resize(np); + + auto& aos = pti.GetArrayOfStructs(); + + // We need to partition the large buffer first + iMultiFab const* bmasks = + (WarpX::n_field_gather_buffer >= WarpX::n_current_deposition_buffer) ? + gather_masks : current_masks; + int i = 0; + const auto& msk = (*bmasks)[pti]; + for (const auto& p : aos) { + const IntVect& iv = Index(p, lev); + inexflag[i++] = msk(iv); + } + + Vector<long> pid; + pid.resize(np); + std::iota(pid.begin(), pid.end(), 0L); + auto sep = std::stable_partition(pid.begin(), pid.end(), + [&inexflag](long id) { return inexflag[id]; }); + + if (WarpX::n_current_deposition_buffer == WarpX::n_field_gather_buffer) { + nfine_current = nfine_gather = std::distance(pid.begin(), sep); + } else if (sep != pid.end()) { + int n_buf; + if (bmasks == gather_masks) { + nfine_gather = std::distance(pid.begin(), sep); + bmasks = current_masks; + n_buf = WarpX::n_current_deposition_buffer; + } else { + nfine_current = std::distance(pid.begin(), sep); + bmasks = gather_masks; + n_buf = WarpX::n_field_gather_buffer; + } + if (n_buf > 0) + { + const auto& msk2 = (*bmasks)[pti]; + for (auto it = sep; it != pid.end(); ++it) { + const long id = *it; + const IntVect& iv = Index(aos[id], lev); + inexflag[id] = msk2(iv); + } + + auto sep2 = std::stable_partition(sep, pid.end(), + [&inexflag](long id) {return inexflag[id];}); + if (bmasks == gather_masks) { + nfine_gather = std::distance(pid.begin(), sep2); + } else { + nfine_current = std::distance(pid.begin(), sep2); + } + } + } + + // only deposit / gather to coarsest grid + if (m_deposit_on_main_grid && lev > 0) { + nfine_current = 0; + } + if (m_gather_from_main_grid && lev > 0) { + nfine_gather = 0; + } + + if (nfine_current != np || nfine_gather != np) + { + + ParticleVector particle_tmp; + particle_tmp.resize(np); + + for (long ip = 0; ip < np; ++ip) { + particle_tmp[ip] = aos[pid[ip]]; + } + std::swap(aos(), particle_tmp); + + RealVector tmp; + tmp.resize(np); + + for (long ip = 0; ip < np; ++ip) { + tmp[ip] = wp[pid[ip]]; + } + std::swap(wp, tmp); + + for (long ip = 0; ip < np; ++ip) { + tmp[ip] = uxp[pid[ip]]; + } + std::swap(uxp, tmp); + + for (long ip = 0; ip < np; ++ip) { + tmp[ip] = uyp[pid[ip]]; + } + std::swap(uyp, tmp); + + for (long ip = 0; ip < np; ++ip) { + tmp[ip] = uzp[pid[ip]]; + } + std::swap(uzp, tmp); + } + +} diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index ee75bc511..7b0d2d1d0 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -68,12 +68,12 @@ public: WarpXParIter (ContainerType& pc, int level); #if (AMREX_SPACEDIM == 2) - void GetPosition (amrex::Cuda::ManagedDeviceVector<amrex::Real>& x, - amrex::Cuda::ManagedDeviceVector<amrex::Real>& y, - amrex::Cuda::ManagedDeviceVector<amrex::Real>& z) const; - void SetPosition (const amrex::Cuda::ManagedDeviceVector<amrex::Real>& x, - const amrex::Cuda::ManagedDeviceVector<amrex::Real>& y, - const amrex::Cuda::ManagedDeviceVector<amrex::Real>& z); + void GetPosition (amrex::Cuda::ManagedDeviceVector<amrex::ParticleReal>& x, + amrex::Cuda::ManagedDeviceVector<amrex::ParticleReal>& y, + amrex::Cuda::ManagedDeviceVector<amrex::ParticleReal>& z) const; + void SetPosition (const amrex::Cuda::ManagedDeviceVector<amrex::ParticleReal>& x, + const amrex::Cuda::ManagedDeviceVector<amrex::ParticleReal>& y, + const amrex::Cuda::ManagedDeviceVector<amrex::ParticleReal>& z); #endif const std::array<RealVector, PIdx::nattribs>& GetAttribs () const { return GetStructOfArrays().GetRealData(); @@ -104,7 +104,7 @@ class WarpXParticleContainer public: friend MultiParticleContainer; - // amrex::StructOfArrays with DiagIdx::nattribs amrex::Real components + // amrex::StructOfArrays with DiagIdx::nattribs amrex::ParticleReal components // and 0 int components for the particle data. using DiagnosticParticleData = amrex::StructOfArrays<DiagIdx::nattribs, 0>; // DiagnosticParticles is a vector, with one element per MR level. @@ -232,17 +232,17 @@ public: amrex::Real maxParticleVelocity(bool local = false); void AddNParticles (int lev, - int n, const amrex::Real* x, const amrex::Real* y, const amrex::Real* z, - const amrex::Real* vx, const amrex::Real* vy, const amrex::Real* vz, - int nattr, const amrex::Real* attr, int uniqueparticles, int id=-1); + int n, const amrex::ParticleReal* x, const amrex::ParticleReal* y, const amrex::ParticleReal* z, + const amrex::ParticleReal* vx, const amrex::ParticleReal* vy, const amrex::ParticleReal* vz, + int nattr, const amrex::ParticleReal* attr, int uniqueparticles, int id=-1); void AddOneParticle (int lev, int grid, int tile, - amrex::Real x, amrex::Real y, amrex::Real z, - std::array<amrex::Real,PIdx::nattribs>& attribs); + amrex::ParticleReal x, amrex::ParticleReal y, amrex::ParticleReal z, + std::array<amrex::ParticleReal,PIdx::nattribs>& attribs); void AddOneParticle (ParticleTileType& particle_tile, - amrex::Real x, amrex::Real y, amrex::Real z, - std::array<amrex::Real,PIdx::nattribs>& attribs); + amrex::ParticleReal x, amrex::ParticleReal y, amrex::ParticleReal z, + std::array<amrex::ParticleReal,PIdx::nattribs>& attribs); virtual void ReadHeader (std::istream& is); @@ -326,7 +326,7 @@ protected: amrex::Vector<amrex::FArrayBox> local_jy; amrex::Vector<amrex::FArrayBox> local_jz; - using DataContainer = amrex::Gpu::ManagedDeviceVector<amrex::Real>; + using DataContainer = amrex::Gpu::ManagedDeviceVector<amrex::ParticleReal>; using PairIndex = std::pair<int, int>; amrex::Vector<DataContainer> m_xp, m_yp, m_zp; diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 176c147da..65a82f233 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -25,7 +25,9 @@ WarpXParIter::WarpXParIter (ContainerType& pc, int level) #if (AMREX_SPACEDIM == 2) void -WarpXParIter::GetPosition (Cuda::ManagedDeviceVector<Real>& x, Cuda::ManagedDeviceVector<Real>& y, Cuda::ManagedDeviceVector<Real>& z) const +WarpXParIter::GetPosition (Gpu::ManagedDeviceVector<ParticleReal>& x, + Gpu::ManagedDeviceVector<ParticleReal>& y, + Gpu::ManagedDeviceVector<ParticleReal>& z) const { amrex::ParIter<0,0,PIdx::nattribs>::GetPosition(x, z); #ifdef WARPX_DIM_RZ @@ -38,17 +40,19 @@ WarpXParIter::GetPosition (Cuda::ManagedDeviceVector<Real>& x, Cuda::ManagedDevi x[i] = x[i]*std::cos(theta[i]); } #else - y.resize(x.size(), std::numeric_limits<Real>::quiet_NaN()); + y.resize(x.size(), std::numeric_limits<ParticleReal>::quiet_NaN()); #endif } void -WarpXParIter::SetPosition (const Cuda::ManagedDeviceVector<Real>& x, const Cuda::ManagedDeviceVector<Real>& y, const Cuda::ManagedDeviceVector<Real>& z) +WarpXParIter::SetPosition (const Gpu::ManagedDeviceVector<ParticleReal>& x, + const Gpu::ManagedDeviceVector<ParticleReal>& y, + const Gpu::ManagedDeviceVector<ParticleReal>& z) { #ifdef WARPX_DIM_RZ auto& attribs = GetAttribs(); auto& theta = attribs[PIdx::theta]; - Cuda::ManagedDeviceVector<Real> r(x.size()); + Gpu::ManagedDeviceVector<ParticleReal> r(x.size()); for (unsigned int i=0 ; i < x.size() ; i++) { theta[i] = std::atan2(y[i], x[i]); r[i] = std::sqrt(x[i]*x[i] + y[i]*y[i]); @@ -132,8 +136,8 @@ WarpXParticleContainer::AllocData () void WarpXParticleContainer::AddOneParticle (int lev, int grid, int tile, - Real x, Real y, Real z, - std::array<Real,PIdx::nattribs>& attribs) + ParticleReal x, ParticleReal y, ParticleReal z, + std::array<ParticleReal,PIdx::nattribs>& attribs) { auto& particle_tile = DefineAndReturnParticleTile(lev, grid, tile); AddOneParticle(particle_tile, x, y, z, attribs); @@ -141,8 +145,8 @@ WarpXParticleContainer::AddOneParticle (int lev, int grid, int tile, void WarpXParticleContainer::AddOneParticle (ParticleTileType& particle_tile, - Real x, Real y, Real z, - std::array<Real,PIdx::nattribs>& attribs) + ParticleReal x, ParticleReal y, ParticleReal z, + std::array<ParticleReal,PIdx::nattribs>& attribs) { ParticleType p; p.id() = ParticleType::NextID(); @@ -171,12 +175,12 @@ WarpXParticleContainer::AddOneParticle (ParticleTileType& particle_tile, void WarpXParticleContainer::AddNParticles (int lev, - int n, const Real* x, const Real* y, const Real* z, - const Real* vx, const Real* vy, const Real* vz, - int nattr, const Real* attr, int uniqueparticles, int id) + int n, const ParticleReal* x, const ParticleReal* y, const ParticleReal* z, + const ParticleReal* vx, const ParticleReal* vy, const ParticleReal* vz, + int nattr, const ParticleReal* attr, int uniqueparticles, int id) { BL_ASSERT(nattr == 1); - const Real* weight = attr; + const ParticleReal* weight = attr; int ibegin, iend; if (uniqueparticles) { @@ -204,7 +208,7 @@ WarpXParticleContainer::AddNParticles (int lev, std::size_t np = iend-ibegin; #ifdef WARPX_DIM_RZ - Vector<Real> theta(np); + Vector<ParticleReal> theta(np); #endif for (int i = ibegin; i < iend; ++i) @@ -253,7 +257,7 @@ WarpXParticleContainer::AddNParticles (int lev, { #ifdef WARPX_DIM_RZ if (comp == PIdx::theta) { - particle_tile.push_back_real(comp, theta.front(), theta.back()); + particle_tile.push_back_real(comp, theta.data(), theta.data() + np); } else { particle_tile.push_back_real(comp, np, 0.0); @@ -362,9 +366,9 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, // CPU, tiling: deposit into local_jx // (same for jx and jz) - Real* AMREX_RESTRICT xp = m_xp[thread_num].dataPtr() + offset; - Real* AMREX_RESTRICT zp = m_zp[thread_num].dataPtr() + offset; - Real* AMREX_RESTRICT yp = m_yp[thread_num].dataPtr() + offset; + ParticleReal* AMREX_RESTRICT xp = m_xp[thread_num].dataPtr() + offset; + ParticleReal* AMREX_RESTRICT zp = m_zp[thread_num].dataPtr() + offset; + ParticleReal* AMREX_RESTRICT yp = m_yp[thread_num].dataPtr() + offset; // Lower corner of tile box physical domain // Note that this includes guard cells since it is after tilebox.ngrow @@ -503,9 +507,9 @@ WarpXParticleContainer::DepositCharge (WarpXParIter& pti, RealVector& wp, // GPU, no tiling: deposit directly in rho // CPU, tiling: deposit into local_rho - Real* AMREX_RESTRICT xp = m_xp[thread_num].dataPtr() + offset; - Real* AMREX_RESTRICT zp = m_zp[thread_num].dataPtr() + offset; - Real* AMREX_RESTRICT yp = m_yp[thread_num].dataPtr() + offset; + ParticleReal* AMREX_RESTRICT xp = m_xp[thread_num].dataPtr() + offset; + ParticleReal* AMREX_RESTRICT zp = m_zp[thread_num].dataPtr() + offset; + ParticleReal* AMREX_RESTRICT yp = m_yp[thread_num].dataPtr() + offset; // Lower corner of tile box physical domain // Note that this includes guard cells since it is after tilebox.ngrow @@ -731,7 +735,7 @@ std::array<Real, 3> WarpXParticleContainer::meanParticleVelocity(bool local) { Real WarpXParticleContainer::maxParticleVelocity(bool local) { - amrex::Real max_v = 0.0; + amrex::ParticleReal max_v = 0.0; for (int lev = 0; lev <= finestLevel(); ++lev) { @@ -745,12 +749,12 @@ Real WarpXParticleContainer::maxParticleVelocity(bool local) { auto& uy = pti.GetAttribs(PIdx::uy); auto& uz = pti.GetAttribs(PIdx::uz); for (unsigned long i = 0; i < ux.size(); i++) { - max_v = std::max(max_v, sqrt(ux[i]*ux[i] + uy[i]*uy[i] + uz[i]*uz[i])); + max_v = std::max(max_v, std::sqrt(ux[i]*ux[i] + uy[i]*uy[i] + uz[i]*uz[i])); } } } - if (!local) ParallelDescriptor::ReduceRealMax(max_v); + if (!local) ParallelAllReduce::Max(max_v, ParallelDescriptor::Communicator()); return max_v; } @@ -819,17 +823,17 @@ WarpXParticleContainer::PushX (int lev, Real dt) ParticleType * AMREX_RESTRICT pstructs = &(pti.GetArrayOfStructs()[0]); // - momenta are stored as a struct of array, in `attribs` auto& attribs = pti.GetAttribs(); - Real* AMREX_RESTRICT ux = attribs[PIdx::ux].dataPtr(); - Real* AMREX_RESTRICT uy = attribs[PIdx::uy].dataPtr(); - Real* AMREX_RESTRICT uz = attribs[PIdx::uz].dataPtr(); + ParticleReal* AMREX_RESTRICT ux = attribs[PIdx::ux].dataPtr(); + ParticleReal* AMREX_RESTRICT uy = attribs[PIdx::uy].dataPtr(); + ParticleReal* AMREX_RESTRICT uz = attribs[PIdx::uz].dataPtr(); #ifdef WARPX_DIM_RZ - Real* AMREX_RESTRICT theta = attribs[PIdx::theta].dataPtr(); + ParticleReal* AMREX_RESTRICT theta = attribs[PIdx::theta].dataPtr(); #endif // Loop over the particles and update their position amrex::ParallelFor( pti.numParticles(), [=] AMREX_GPU_DEVICE (long i) { ParticleType& p = pstructs[i]; // Particle object that gets updated - Real x, y, z; // Temporary variables + ParticleReal x, y, z; // Temporary variables #ifndef WARPX_DIM_RZ GetPosition( x, y, z, p ); // Initialize x, y, z UpdatePosition( x, y, z, ux[i], uy[i], uz[i], dt); diff --git a/Source/Particles/deposit_cic.F90 b/Source/Particles/deposit_cic.F90 index 1fe80016f..2831ce96b 100644 --- a/Source/Particles/deposit_cic.F90 +++ b/Source/Particles/deposit_cic.F90 @@ -1,7 +1,7 @@ module warpx_ES_deposit_cic use iso_c_binding - use amrex_fort_module, only : amrex_real + use amrex_fort_module, only : amrex_real, amrex_particle_real implicit none @@ -28,8 +28,8 @@ contains ng) & bind(c,name='warpx_deposit_cic_3d') integer, value, intent(in) :: ns, np - real(amrex_real), intent(in) :: particles(ns,np) - real(amrex_real), intent(in) :: weights(np) + real(amrex_particle_real), intent(in) :: particles(ns,np) + real(amrex_particle_real), intent(in) :: weights(np) real(amrex_real), intent(in) :: charge integer, intent(in) :: lo(3) integer, intent(in) :: hi(3) @@ -86,8 +86,8 @@ contains ng) & bind(c,name='warpx_deposit_cic_2d') integer, value, intent(in) :: ns, np - real(amrex_real), intent(in) :: particles(ns,np) - real(amrex_real), intent(in) :: weights(np) + real(amrex_particle_real), intent(in) :: particles(ns,np) + real(amrex_particle_real), intent(in) :: weights(np) real(amrex_real), intent(in) :: charge integer, intent(in) :: lo(2) integer, intent(in) :: hi(2) diff --git a/Source/Particles/interpolate_cic.F90 b/Source/Particles/interpolate_cic.F90 index 005ab35f4..3eb361d2f 100644 --- a/Source/Particles/interpolate_cic.F90 +++ b/Source/Particles/interpolate_cic.F90 @@ -1,7 +1,7 @@ module warpx_ES_interpolate_cic use iso_c_binding - use amrex_fort_module, only : amrex_real + use amrex_fort_module, only : amrex_real, amrex_particle_real implicit none @@ -31,7 +31,7 @@ contains lo, hi, plo, dx, ng) & bind(c,name='warpx_interpolate_cic_3d') integer, value, intent(in) :: ns, np - real(amrex_real), intent(in) :: particles(ns,np) + real(amrex_particle_real), intent(in) :: particles(ns,np) real(amrex_real), intent(inout) :: Ex_p(np), Ey_p(np), Ez_p(np) integer, intent(in) :: ng integer, intent(in) :: lo(3) @@ -103,7 +103,7 @@ contains lo, hi, plo, dx, ng) & bind(c,name='warpx_interpolate_cic_2d') integer, value, intent(in) :: ns, np - real(amrex_real), intent(in) :: particles(ns,np) + real(amrex_particle_real), intent(in) :: particles(ns,np) real(amrex_real), intent(inout) :: Ex_p(np), Ey_p(np) integer, intent(in) :: ng integer, intent(in) :: lo(2) @@ -157,7 +157,7 @@ contains plo, ng, lev) & bind(c,name='warpx_interpolate_cic_two_levels_3d') integer, value, intent(in) :: ns, np - real(amrex_real), intent(in) :: particles(ns,np) + real(amrex_particle_real), intent(in) :: particles(ns,np) real(amrex_real), intent(inout) :: Ex_p(np), Ey_p(np), Ez_p(np) integer, intent(in) :: ng, lev integer, intent(in) :: lo(3), hi(3) @@ -290,7 +290,7 @@ contains plo, ng, lev) & bind(c,name='warpx_interpolate_cic_two_levels_2d') integer, value, intent(in) :: ns, np - real(amrex_real), intent(in) :: particles(ns,np) + real(amrex_particle_real), intent(in) :: particles(ns,np) real(amrex_real), intent(inout) :: Ex_p(np), Ey_p(np) integer, intent(in) :: ng, lev integer, intent(in) :: lo(2), hi(2) diff --git a/Source/Particles/push_particles_ES.F90 b/Source/Particles/push_particles_ES.F90 index 53dd3c181..b84f48d5f 100644 --- a/Source/Particles/push_particles_ES.F90 +++ b/Source/Particles/push_particles_ES.F90 @@ -1,7 +1,7 @@ module warpx_ES_push_particles use iso_c_binding - use amrex_fort_module, only : amrex_real + use amrex_fort_module, only : amrex_real, amrex_particle_real implicit none @@ -36,8 +36,8 @@ contains prob_lo, prob_hi) & bind(c,name='warpx_push_leapfrog_3d') integer, value, intent(in) :: ns, np - real(amrex_real), intent(inout) :: particles(ns,np) - real(amrex_real), intent(inout) :: vx_p(np), vy_p(np), vz_p(np) + real(amrex_particle_real), intent(inout) :: particles(ns,np) + real(amrex_particle_real), intent(inout) :: vx_p(np), vy_p(np), vz_p(np) real(amrex_real), intent(in) :: Ex_p(np), Ey_p(np), Ez_p(np) real(amrex_real), intent(in) :: charge real(amrex_real), intent(in) :: mass @@ -100,8 +100,8 @@ contains prob_lo, prob_hi) & bind(c,name='warpx_push_leapfrog_2d') integer, value, intent(in) :: ns, np - real(amrex_real), intent(inout) :: particles(ns,np) - real(amrex_real), intent(inout) :: vx_p(np), vy_p(np) + real(amrex_particle_real), intent(inout) :: particles(ns,np) + real(amrex_particle_real), intent(inout) :: vx_p(np), vy_p(np) real(amrex_real), intent(in) :: Ex_p(np), Ey_p(np) real(amrex_real), intent(in) :: charge real(amrex_real), intent(in) :: mass @@ -167,8 +167,8 @@ contains prob_lo, prob_hi) & bind(c,name='warpx_push_leapfrog_positions_3d') integer, value, intent(in) :: ns, np - real(amrex_real), intent(inout) :: particles(ns,np) - real(amrex_real), intent(inout) :: vx_p(np), vy_p(np), vz_p(np) + real(amrex_particle_real), intent(inout) :: particles(ns,np) + real(amrex_particle_real), intent(inout) :: vx_p(np), vy_p(np), vz_p(np) real(amrex_real), intent(in) :: dt real(amrex_real), intent(in) :: prob_lo(3), prob_hi(3) @@ -219,8 +219,8 @@ contains prob_lo, prob_hi) & bind(c,name='warpx_push_leapfrog_positions_2d') integer, value, intent(in) :: ns, np - real(amrex_real), intent(inout) :: particles(ns,np) - real(amrex_real), intent(inout) :: vx_p(np), vy_p(np) + real(amrex_particle_real), intent(inout) :: particles(ns,np) + real(amrex_particle_real), intent(inout) :: vx_p(np), vy_p(np) real(amrex_real), intent(in) :: dt real(amrex_real), intent(in) :: prob_lo(2), prob_hi(2) diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index a844e7aa0..116600c8b 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -103,7 +103,7 @@ IntVect WarpX::jz_nodal_flag(1,0); // z is the second dimension to 2D AMReX IntVect WarpX::filter_npass_each_dir(1); -int WarpX::n_field_gather_buffer = 0; +int WarpX::n_field_gather_buffer = -1; int WarpX::n_current_deposition_buffer = -1; int WarpX::do_nodal = false; @@ -149,8 +149,8 @@ WarpX::WarpX () #endif t_new.resize(nlevs_max, 0.0); - t_old.resize(nlevs_max, -1.e100); - dt.resize(nlevs_max, 1.e100); + t_old.resize(nlevs_max, std::numeric_limits<Real>::lowest()); + dt.resize(nlevs_max, std::numeric_limits<Real>::max()); // Particle Container mypc = std::unique_ptr<MultiParticleContainer> (new MultiParticleContainer(this)); @@ -730,11 +730,19 @@ WarpX::AllocLevelData (int lev, const BoxArray& ba, const DistributionMapping& d if (mypc->nSpeciesDepositOnMainGrid() && n_current_deposition_buffer == 0) { n_current_deposition_buffer = 1; + // This forces the allocation of buffers and allows the code associated + // with buffers to run. But the buffer size of `1` is in fact not used, + // `deposit_on_main_grid` forces all particles (whether or not they + // are in buffers) to deposition on the main grid. } if (n_current_deposition_buffer < 0) { n_current_deposition_buffer = ngJ.max(); } + if (n_field_gather_buffer < 0) { + // Field gather buffer should be larger than current deposition buffers + n_field_gather_buffer = n_current_deposition_buffer + 1; + } int ngF = (do_moving_window) ? 2 : 0; // CKC solver requires one additional guard cell @@ -988,7 +996,7 @@ WarpX::LowerCorner(const Box& bx, int lev) #if (AMREX_SPACEDIM == 3) return { xyzmin[0], xyzmin[1], xyzmin[2] }; #elif (AMREX_SPACEDIM == 2) - return { xyzmin[0], -1.e100, xyzmin[1] }; + return { xyzmin[0], std::numeric_limits<Real>::lowest(), xyzmin[1] }; #endif } @@ -1000,7 +1008,7 @@ WarpX::UpperCorner(const Box& bx, int lev) #if (AMREX_SPACEDIM == 3) return { xyzmax[0], xyzmax[1], xyzmax[2] }; #elif (AMREX_SPACEDIM == 2) - return { xyzmax[0], 1.e100, xyzmax[1] }; + return { xyzmax[0], std::numeric_limits<Real>::max(), xyzmax[1] }; #endif } diff --git a/Source/main.cpp b/Source/main.cpp index 551be9c30..cb183bc8d 100644 --- a/Source/main.cpp +++ b/Source/main.cpp @@ -12,6 +12,7 @@ using namespace amrex; int main(int argc, char* argv[]) { +#if defined AMREX_USE_MPI #if defined(_OPENMP) && defined(WARPX_USE_PSATD) int provided; MPI_Init_thread(&argc, &argv, MPI_THREAD_FUNNELED, &provided); @@ -19,6 +20,7 @@ int main(int argc, char* argv[]) #else MPI_Init(&argc, &argv); #endif +#endif amrex::Initialize(argc,argv); @@ -48,5 +50,7 @@ int main(int argc, char* argv[]) BL_PROFILE_VAR_STOP(pmain); amrex::Finalize(); +#if defined AMREX_USE_MPI MPI_Finalize(); +#endif } |