diff options
Diffstat (limited to 'Source/Particles/WarpXParticleContainer.cpp')
| -rw-r--r-- | Source/Particles/WarpXParticleContainer.cpp | 277 |
1 files changed, 142 insertions, 135 deletions
diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index a9e0a7418..c52e0a6d0 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -21,13 +21,35 @@ void WarpXParIter::GetPosition (Cuda::ManagedDeviceVector<Real>& x, Cuda::ManagedDeviceVector<Real>& y, Cuda::ManagedDeviceVector<Real>& z) const { amrex::ParIter<0,0,PIdx::nattribs>::GetPosition(x, z); +#ifdef WARPX_RZ + const auto& attribs = GetAttribs(); + const auto& theta = attribs[PIdx::theta]; + y.resize(x.size()); + for (unsigned int i=0 ; i < x.size() ; i++) { + // The x stored in the particles is actually the radius + y[i] = x[i]*std::sin(theta[i]); + x[i] = x[i]*std::cos(theta[i]); + } +#else y.resize(x.size(), std::numeric_limits<Real>::quiet_NaN()); +#endif } void WarpXParIter::SetPosition (const Cuda::ManagedDeviceVector<Real>& x, const Cuda::ManagedDeviceVector<Real>& y, const Cuda::ManagedDeviceVector<Real>& z) { +#ifdef WARPX_RZ + auto& attribs = GetAttribs(); + auto& theta = attribs[PIdx::theta]; + Cuda::DeviceVector<Real> 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]); + } + amrex::ParIter<0,0,PIdx::nattribs>::SetPosition(r, z); +#else amrex::ParIter<0,0,PIdx::nattribs>::SetPosition(x, z); +#endif } #endif @@ -90,7 +112,7 @@ WarpXParticleContainer::AllocData () void WarpXParticleContainer::AddOneParticle (int lev, int grid, int tile, Real x, Real y, Real z, - const std::array<Real,PIdx::nattribs>& attribs) + std::array<Real,PIdx::nattribs>& attribs) { auto& particle_tile = GetParticles(lev)[std::make_pair(grid,tile)]; AddOneParticle(particle_tile, x, y, z, attribs); @@ -99,7 +121,7 @@ WarpXParticleContainer::AddOneParticle (int lev, int grid, int tile, void WarpXParticleContainer::AddOneParticle (ParticleTileType& particle_tile, Real x, Real y, Real z, - const std::array<Real,PIdx::nattribs>& attribs) + std::array<Real,PIdx::nattribs>& attribs) { ParticleType p; p.id() = ParticleType::NextID(); @@ -109,6 +131,10 @@ WarpXParticleContainer::AddOneParticle (ParticleTileType& particle_tile, p.pos(1) = y; p.pos(2) = z; #elif (AMREX_SPACEDIM == 2) +#ifdef WARPX_RZ + attribs[PIdx::theta] = std::atan2(y, x); + x = std::sqrt(x*x + y*y); +#endif p.pos(0) = x; p.pos(1) = z; #endif @@ -149,6 +175,12 @@ WarpXParticleContainer::AddNParticles (int lev, std::pair<int,int> key {0,0}; auto& particle_tile = GetParticles(lev)[key]; + std::size_t np = iend-ibegin; + +#ifdef WARPX_RZ + Vector<Real> theta(np); +#endif + for (int i = ibegin; i < iend; ++i) { ParticleType p; @@ -164,14 +196,17 @@ WarpXParticleContainer::AddNParticles (int lev, p.pos(1) = y[i]; p.pos(2) = z[i]; #elif (AMREX_SPACEDIM == 2) +#ifdef WARPX_RZ + theta[i-ibegin] = std::atan2(y[i], x[i]); + p.pos(0) = std::sqrt(x[i]*x[i] + y[i]*y[i]); +#else p.pos(0) = x[i]; +#endif p.pos(1) = z[i]; #endif particle_tile.push_back(p); } - std::size_t np = iend-ibegin; - if (np > 0) { particle_tile.push_back_real(PIdx::w , weight + ibegin, weight + iend); @@ -181,7 +216,16 @@ WarpXParticleContainer::AddNParticles (int lev, for (int comp = PIdx::uz+1; comp < PIdx::nattribs; ++comp) { +#ifdef WARPX_RZ + if (comp == PIdx::theta) { + particle_tile.push_back_real(comp, theta.front(), theta.back()); + } + else { + particle_tile.push_back_real(comp, np, 0.0); + } +#else particle_tile.push_back_real(comp, np, 0.0); +#endif } } @@ -191,12 +235,12 @@ WarpXParticleContainer::AddNParticles (int lev, void WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, - RealVector& wp, RealVector& uxp, - RealVector& uyp, RealVector& uzp, - MultiFab& jx, MultiFab& jy, MultiFab& jz, - MultiFab* cjx, MultiFab* cjy, MultiFab* cjz, - const long np_current, const long np, - int thread_num, int lev, Real dt ) + RealVector& wp, RealVector& uxp, + RealVector& uyp, RealVector& uzp, + MultiFab& jx, MultiFab& jy, MultiFab& jz, + MultiFab* cjx, MultiFab* cjy, MultiFab* cjz, + const long np_current, const long np, + int thread_num, int lev, Real dt ) { Real *jx_ptr, *jy_ptr, *jz_ptr; const std::array<Real,3>& xyzmin_tile = WarpX::LowerCorner(pti.tilebox(), lev); @@ -219,7 +263,16 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, // Deposit charge for particles that are not in the current buffers if (np_current > 0) - { + { +#ifdef AMREX_USE_GPU + jx_ptr = jx[pti].dataPtr(); + jy_ptr = jy[pti].dataPtr(); + jz_ptr = jz[pti].dataPtr(); + + auto jxntot = jx[pti].length(); + auto jyntot = jy[pti].length(); + auto jzntot = jz[pti].length(); +#else tbx.grow(ngJ); tby.grow(ngJ); tbz.grow(ngJ); @@ -231,31 +284,15 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, jx_ptr = local_jx[thread_num].dataPtr(); jy_ptr = local_jy[thread_num].dataPtr(); jz_ptr = local_jz[thread_num].dataPtr(); - - auto jxarr = local_jx[thread_num].array(); - amrex::ParallelFor(tbx, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - jxarr(i,j,k) = 0.0; - }); - - auto jyarr = local_jy[thread_num].array(); - amrex::ParallelFor(tby, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - jyarr(i,j,k) = 0.0; - }); - - auto jzarr = local_jz[thread_num].array(); - amrex::ParallelFor(tbz, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - jzarr(i,j,k) = 0.0; - }); + + local_jx[thread_num].setVal(0.0); + local_jy[thread_num].setVal(0.0); + local_jz[thread_num].setVal(0.0); auto jxntot = local_jx[thread_num].length(); auto jyntot = local_jy[thread_num].length(); auto jzntot = local_jz[thread_num].length(); +#endif BL_PROFILE_VAR_START(blp_pxr_cd); if (j_is_nodal) { @@ -334,45 +371,46 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, &dt, &dx[0], &dx[1], &dx[2], &WarpX::nox,&WarpX::noy,&WarpX::noz, &lvect,&WarpX::current_deposition_algo); + +#ifdef WARPX_RZ + warpx_current_deposition_rz_volume_scaling( + jx_ptr, &ngJ, jxntot.getVect(), + jy_ptr, &ngJ, jyntot.getVect(), + jz_ptr, &ngJ, jzntot.getVect(), + &xyzmin[0], &dx[0]); +#endif } + BL_PROFILE_VAR_STOP(blp_pxr_cd); - + +#ifndef AMREX_USE_GPU BL_PROFILE_VAR_START(blp_accumulate); - - const auto local_jx_arr = local_jx[thread_num].array(); - auto global_jx_arr = jx.array(pti); - amrex::ParallelFor(tbx, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - Gpu::Atomic::Add(&global_jx_arr(i, j, k), local_jx_arr(i, j, k)); - }); - - const auto local_jy_arr = local_jy[thread_num].array(); - auto global_jy_arr = jy.array(pti); - amrex::ParallelFor(tby, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - Gpu::Atomic::Add(&global_jy_arr(i, j, k), local_jy_arr(i, j, k)); - }); - - const auto local_jz_arr = local_jz[thread_num].array(); - auto global_jz_arr = jz.array(pti); - amrex::ParallelFor(tbz, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - Gpu::Atomic::Add(&global_jz_arr(i, j, k), local_jz_arr(i, j, k)); - }); + + jx[pti].atomicAdd(local_jx[thread_num], tbx, tbx, 0, 0, 1); + jy[pti].atomicAdd(local_jy[thread_num], tby, tby, 0, 0, 1); + jz[pti].atomicAdd(local_jz[thread_num], tbz, tbz, 0, 0, 1); BL_PROFILE_VAR_STOP(blp_accumulate); - } +#endif + } // Deposit charge for particles that are in the current buffers if (np_current < np) - { + { const IntVect& ref_ratio = WarpX::RefRatio(lev-1); const Box& ctilebox = amrex::coarsen(pti.tilebox(),ref_ratio); const std::array<Real,3>& cxyzmin_tile = WarpX::LowerCorner(ctilebox, lev-1); +#ifdef AMREX_USE_GPU + jx_ptr = (*cjx)[pti].dataPtr(); + jy_ptr = (*cjy)[pti].dataPtr(); + jz_ptr = (*cjz)[pti].dataPtr(); + + auto jxntot = jx[pti].length(); + auto jyntot = jy[pti].length(); + auto jzntot = jz[pti].length(); +#else + tbx = amrex::convert(ctilebox, WarpX::jx_nodal_flag); tby = amrex::convert(ctilebox, WarpX::jy_nodal_flag); tbz = amrex::convert(ctilebox, WarpX::jz_nodal_flag); @@ -388,31 +426,15 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, jy_ptr = local_jy[thread_num].dataPtr(); jz_ptr = local_jz[thread_num].dataPtr(); - auto jxarr = local_jx[thread_num].array(); - amrex::ParallelFor(tbx, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - jxarr(i,j,k) = 0.0; - }); - - auto jyarr = local_jy[thread_num].array(); - amrex::ParallelFor(tby, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - jyarr(i,j,k) = 0.0; - }); - - auto jzarr = local_jz[thread_num].array(); - amrex::ParallelFor(tbz, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - jzarr(i,j,k) = 0.0; - }); + local_jx[thread_num].setVal(0.0); + local_jy[thread_num].setVal(0.0); + local_jz[thread_num].setVal(0.0); auto jxntot = local_jx[thread_num].length(); auto jyntot = local_jy[thread_num].length(); auto jzntot = local_jz[thread_num].length(); - +#endif + long ncrse = np - np_current; BL_PROFILE_VAR_START(blp_pxr_cd); if (j_is_nodal) { @@ -493,36 +515,26 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, &dt, &cdx[0], &cdx[1], &cdx[2], &WarpX::nox,&WarpX::noy,&WarpX::noz, &lvect,&WarpX::current_deposition_algo); +#ifdef WARPX_RZ + warpx_current_deposition_rz_volume_scaling( + jx_ptr, &ngJ, jxntot.getVect(), + jy_ptr, &ngJ, jyntot.getVect(), + jz_ptr, &ngJ, jzntot.getVect(), + &xyzmin[0], &dx[0]); +#endif } + BL_PROFILE_VAR_STOP(blp_pxr_cd); +#ifndef AMREX_USE_GPU BL_PROFILE_VAR_START(blp_accumulate); - const auto local_jx_arr = local_jx[thread_num].array(); - auto global_jx_arr = cjx->array(pti); - amrex::ParallelFor(tbx, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - Gpu::Atomic::Add(&global_jx_arr(i, j, k), local_jx_arr(i, j, k)); - }); - - const auto local_jy_arr = local_jy[thread_num].array(); - auto global_jy_arr = cjy->array(pti); - amrex::ParallelFor(tby, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - Gpu::Atomic::Add(&global_jy_arr(i, j, k), local_jy_arr(i, j, k)); - }); - - const auto local_jz_arr = local_jz[thread_num].array(); - auto global_jz_arr = cjz->array(pti); - amrex::ParallelFor(tbz, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - Gpu::Atomic::Add(&global_jz_arr(i, j, k), local_jz_arr(i, j, k)); - }); + (*cjx)[pti].atomicAdd(local_jx[thread_num], tbx, tbx, 0, 0, 1); + (*cjy)[pti].atomicAdd(local_jy[thread_num], tby, tby, 0, 0, 1); + (*cjz)[pti].atomicAdd(local_jz[thread_num], tbz, tbz, 0, 0, 1); BL_PROFILE_VAR_STOP(blp_accumulate); +#endif } }; @@ -549,20 +561,22 @@ WarpXParticleContainer::DepositCharge ( WarpXParIter& pti, RealVector& wp, // Deposit charge for particles that are not in the current buffers if (np_current > 0) - { + { const std::array<Real, 3>& xyzmin = xyzmin_tile; + +#ifdef AMREX_USE_GPU + data_ptr = (*rhomf)[pti].dataPtr(); + auto rholen = (*rhomf)[pti].length(); +#else tile_box.grow(ngRho); local_rho[thread_num].resize(tile_box); - auto rhoarr = local_rho[thread_num].array(); - amrex::ParallelFor(tile_box, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - rhoarr(i,j,k) = 0.0; - }); - data_ptr = local_rho[thread_num].dataPtr(); auto rholen = local_rho[thread_num].length(); + + local_rho[thread_num].setVal(0.0); +#endif + #if (AMREX_SPACEDIM == 3) const long nx = rholen[0]-1-2*ngRho; const long ny = rholen[1]-1-2*ngRho; @@ -586,39 +600,36 @@ WarpXParticleContainer::DepositCharge ( WarpXParIter& pti, RealVector& wp, &lvect, &WarpX::charge_deposition_algo); BL_PROFILE_VAR_STOP(blp_pxr_chd); +#ifndef AMREX_USE_GPU BL_PROFILE_VAR_START(blp_accumulate); - const auto local_rho_arr = local_rho[thread_num].array(); - auto global_rho_arr = rhomf->array(pti); - amrex::ParallelFor(tile_box, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - Gpu::Atomic::Add(&global_rho_arr(i, j, k, icomp), local_rho_arr(i, j, k)); - }); + (*rhomf)[pti].atomicAdd(local_rho[thread_num], tile_box, tile_box, 0, icomp, 1); BL_PROFILE_VAR_STOP(blp_accumulate); - } +#endif + } // Deposit charge for particles that are in the current buffers if (np_current < np) - { + { const IntVect& ref_ratio = WarpX::RefRatio(lev-1); const Box& ctilebox = amrex::coarsen(pti.tilebox(), ref_ratio); const std::array<Real,3>& cxyzmin_tile = WarpX::LowerCorner(ctilebox, lev-1); +#ifdef AMREX_USE_GPU + data_ptr = (*crhomf)[pti].dataPtr(); + auto rholen = (*crhomf)[pti].length(); +#else tile_box = amrex::convert(ctilebox, IntVect::TheUnitVector()); tile_box.grow(ngRho); local_rho[thread_num].resize(tile_box); - auto rhoarr = local_rho[thread_num].array(); - amrex::ParallelFor(tile_box, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - rhoarr(i,j,k) = 0.0; - }); - data_ptr = local_rho[thread_num].dataPtr(); auto rholen = local_rho[thread_num].length(); + + local_rho[thread_num].setVal(0.0); +#endif + #if (AMREX_SPACEDIM == 3) const long nx = rholen[0]-1-2*ngRho; const long ny = rholen[1]-1-2*ngRho; @@ -644,17 +655,13 @@ WarpXParticleContainer::DepositCharge ( WarpXParIter& pti, RealVector& wp, &lvect, &WarpX::charge_deposition_algo); BL_PROFILE_VAR_STOP(blp_pxr_chd); +#ifndef AMREX_USE_GPU BL_PROFILE_VAR_START(blp_accumulate); - const auto local_rho_arr = local_rho[thread_num].array(); - auto global_rho_arr = crhomf->array(pti); - amrex::ParallelFor(tile_box, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - Gpu::Atomic::Add(&global_rho_arr(i, j, k, icomp), local_rho_arr(i, j, k)); - }); - + (*crhomf)[pti].atomicAdd(local_rho[thread_num], tile_box, tile_box, 0, icomp, 1); + BL_PROFILE_VAR_STOP(blp_accumulate); +#endif } }; |