From 5cc13f04eaf1a41c22245129bf083f38e7c596ef Mon Sep 17 00:00:00 2001 From: Dave Grote Date: Thu, 8 Aug 2019 15:48:08 -0700 Subject: Converted DepositCharge to c++ (and removed duplication of code) --- Source/Particles/WarpXParticleContainer.cpp | 243 ++++++++-------------------- 1 file changed, 70 insertions(+), 173 deletions(-) (limited to 'Source/Particles/WarpXParticleContainer.cpp') diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index c5c6afc19..ec0f78564 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -12,6 +12,7 @@ #include #include #include +#include using namespace amrex; @@ -552,140 +553,91 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, } void -WarpXParticleContainer::DepositCharge ( WarpXParIter& pti, RealVector& wp, - MultiFab* rhomf, MultiFab* crhomf, int icomp, - const long np_current, - const long np, int thread_num, int lev ) +WarpXParticleContainer::DepositCharge (WarpXParIter& pti, RealVector& wp, + MultiFab* rho, int icomp, + const long offset, const long np_to_depose, + int thread_num, int lev, int depos_lev) { + AMREX_ALWAYS_ASSERT_WITH_MESSAGE((depos_lev==(lev-1)) || + (depos_lev==(lev )), + "Deposition buffers only work for lev-1"); - BL_PROFILE_VAR_NS("PICSAR::ChargeDeposition", blp_pxr_chd); - BL_PROFILE_VAR_NS("PPC::Evolve::Accumulate", blp_accumulate); - - const std::array& xyzmin_tile = WarpX::LowerCorner(pti.tilebox(), lev); - const long lvect = 8; + // If no particles, do not do anything + if (np_to_depose == 0) return; - long ngRho = rhomf->nGrow(); - Real* data_ptr; - Box tile_box = convert(pti.tilebox(), IntVect::TheUnitVector()); + const long ngRho = rho->nGrow(); + const std::array& dx = WarpX::CellSize(std::max(depos_lev,0)); + Real q = this->charge; - const std::array& dx = WarpX::CellSize(lev); - const std::array& cdx = WarpX::CellSize(std::max(lev-1,0)); + BL_PROFILE_VAR_NS("PPC::ChargeDeposition", blp_ppc_chd); + BL_PROFILE_VAR_NS("PPC::Evolve::Accumulate", blp_accumulate); - // Deposit charge for particles that are not in the current buffers - if (np_current > 0) - { - const std::array& xyzmin = xyzmin_tile; + // Get tile box where charge is deposited. + // The tile box is different when depositing in the buffers (depos_lev const& rho_arr = rho->array(pti); #else - tile_box.grow(ngRho); - local_rho[thread_num].resize(tile_box); + // Tiling is on: rho_ptr points to local_rho[thread_num] + Box tb = amrex::convert(tilebox, IntVect::TheUnitVector()); - data_ptr = local_rho[thread_num].dataPtr(); - auto rholen = local_rho[thread_num].length(); + local_rho[thread_num].resize(tb); - local_rho[thread_num].setVal(0.0); -#endif + // local_rho[thread_num] is set to zero + local_rho[thread_num].setVal(0.0); -#if (AMREX_SPACEDIM == 3) - const long nx = rholen[0]-1-2*ngRho; - const long ny = rholen[1]-1-2*ngRho; - const long nz = rholen[2]-1-2*ngRho; -#else - const long nx = rholen[0]-1-2*ngRho; - const long ny = 0; - const long nz = rholen[1]-1-2*ngRho; -#endif - BL_PROFILE_VAR_START(blp_pxr_chd); - warpx_charge_deposition(data_ptr, &np_current, - m_xp[thread_num].dataPtr(), - m_yp[thread_num].dataPtr(), - m_zp[thread_num].dataPtr(), - wp.dataPtr(), - &this->charge, - &xyzmin[0], &xyzmin[1], &xyzmin[2], - &dx[0], &dx[1], &dx[2], &nx, &ny, &nz, - &ngRho, &ngRho, &ngRho, - &WarpX::nox,&WarpX::noy,&WarpX::noz, - &lvect, &WarpX::charge_deposition_algo); -#ifdef WARPX_DIM_RZ - warpx_charge_deposition_rz_volume_scaling( - data_ptr, &ngRho, rholen.getVect(), - &xyzmin[0], &dx[0]); + Array4 const& rho_arr = local_rho[thread_num].array(); #endif - BL_PROFILE_VAR_STOP(blp_pxr_chd); - -#ifndef AMREX_USE_GPU - BL_PROFILE_VAR_START(blp_accumulate); - - (*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& cxyzmin_tile = WarpX::LowerCorner(ctilebox, lev-1); - -#ifdef AMREX_USE_GPU - data_ptr = (*crhomf)[pti].dataPtr(icomp); - auto rholen = (*crhomf)[pti].length(); -#else - tile_box = amrex::convert(ctilebox, IntVect::TheUnitVector()); - tile_box.grow(ngRho); - local_rho[thread_num].resize(tile_box); + // GPU, no tiling: deposit directly in rho + // CPU, tiling: deposit into local_rho - data_ptr = local_rho[thread_num].dataPtr(); - auto rholen = local_rho[thread_num].length(); - - local_rho[thread_num].setVal(0.0); -#endif + 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; -#if (AMREX_SPACEDIM == 3) - const long nx = rholen[0]-1-2*ngRho; - const long ny = rholen[1]-1-2*ngRho; - const long nz = rholen[2]-1-2*ngRho; -#else - const long nx = rholen[0]-1-2*ngRho; - const long ny = 0; - const long nz = rholen[1]-1-2*ngRho; -#endif + // Lower corner of tile box physical domain + // Note that this includes guard cells since it is after tilebox.ngrow + const std::array& xyzmin = WarpX::LowerCorner(tilebox, depos_lev); + // Indices of the lower bound + const Dim3 lo = lbound(tilebox); - long ncrse = np - np_current; - BL_PROFILE_VAR_START(blp_pxr_chd); - warpx_charge_deposition(data_ptr, &ncrse, - m_xp[thread_num].dataPtr() + np_current, - m_yp[thread_num].dataPtr() + np_current, - m_zp[thread_num].dataPtr() + np_current, - wp.dataPtr() + np_current, - &this->charge, - &cxyzmin_tile[0], &cxyzmin_tile[1], &cxyzmin_tile[2], - &cdx[0], &cdx[1], &cdx[2], &nx, &ny, &nz, - &ngRho, &ngRho, &ngRho, - &WarpX::nox,&WarpX::noy,&WarpX::noz, - &lvect, &WarpX::charge_deposition_algo); + BL_PROFILE_VAR_START(blp_ppc_chd); + if (WarpX::nox == 1){ + doChargeDepositionShapeN<1>(xp, yp, zp, wp.dataPtr()+offset, rho_arr, + np_to_depose, dx, xyzmin, lo, q); + } else if (WarpX::nox == 2){ + doChargeDepositionShapeN<2>(xp, yp, zp, wp.dataPtr()+offset, rho_arr, + np_to_depose, dx, xyzmin, lo, q); + } else if (WarpX::nox == 3){ + doChargeDepositionShapeN<3>(xp, yp, zp, wp.dataPtr()+offset, rho_arr, + np_to_depose, dx, xyzmin, lo, q); + } #ifdef WARPX_DIM_RZ - warpx_charge_deposition_rz_volume_scaling( - data_ptr, &ngRho, rholen.getVect(), - &cxyzmin_tile[0], &cdx[0]); + warpx_charge_deposition_rz_volume_scaling( + data_ptr, &ngRho, rholen.getVect(), + &xyzmin[0], &dx[0]); #endif - BL_PROFILE_VAR_STOP(blp_pxr_chd); + BL_PROFILE_VAR_STOP(blp_ppc_chd); #ifndef AMREX_USE_GPU - BL_PROFILE_VAR_START(blp_accumulate); + BL_PROFILE_VAR_START(blp_accumulate); - (*crhomf)[pti].atomicAdd(local_rho[thread_num], tile_box, tile_box, 0, icomp, 1); + (*rho)[pti].atomicAdd(local_rho[thread_num], tb, tb, 0, icomp, 1); - BL_PROFILE_VAR_STOP(blp_accumulate); + BL_PROFILE_VAR_STOP(blp_accumulate); #endif - } -}; +} void WarpXParticleContainer::DepositCharge (Vector >& rho, bool local) @@ -762,8 +714,6 @@ WarpXParticleContainer::GetChargeDensity (int lev, bool local) BoxArray nba = ba; nba.surroundingNodes(); - const std::array& dx = WarpX::CellSize(lev); - const int ng = WarpX::nox; auto rho = std::unique_ptr(new MultiFab(nba,dm,1,ng)); @@ -773,74 +723,21 @@ WarpXParticleContainer::GetChargeDensity (int lev, bool local) #pragma omp parallel { #endif - Cuda::ManagedDeviceVector xp, yp, zp; #ifdef _OPENMP - FArrayBox rho_loc; + int thread_num = omp_get_thread_num(); +#else + int thread_num = 0; #endif for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) { + const long np = pti.numParticles(); auto& wp = pti.GetAttribs(PIdx::w); - const long np = pti.numParticles(); + pti.GetPosition(m_xp[thread_num], m_yp[thread_num], m_zp[thread_num]); - pti.GetPosition(xp, yp, zp); - - // Data on the grid - Real* data_ptr; - FArrayBox& rhofab = (*rho)[pti]; -#ifdef _OPENMP - const std::array& xyzmin_tile = WarpX::LowerCorner(pti.tilebox(), lev); - Box tile_box = convert(pti.tilebox(), IntVect::TheUnitVector()); - const std::array& xyzmin = xyzmin_tile; - tile_box.grow(ng); - rho_loc.resize(tile_box); - rho_loc = 0.0; - data_ptr = rho_loc.dataPtr(); - auto rholen = rho_loc.length(); -#else - const Box& box = pti.validbox(); - const std::array& xyzmin_grid = WarpX::LowerCorner(box, lev); - const std::array& xyzmin = xyzmin_grid; - data_ptr = rhofab.dataPtr(); - auto rholen = rhofab.length(); -#endif - -#if (AMREX_SPACEDIM == 3) - const long nx = rholen[0]-1-2*ng; - const long ny = rholen[1]-1-2*ng; - const long nz = rholen[2]-1-2*ng; -#else - const long nx = rholen[0]-1-2*ng; - const long ny = 0; - const long nz = rholen[1]-1-2*ng; -#endif - - long nxg = ng; - long nyg = ng; - long nzg = ng; - long lvect = 8; - - warpx_charge_deposition(data_ptr, - &np, - xp.dataPtr(), - yp.dataPtr(), - zp.dataPtr(), wp.dataPtr(), - &this->charge, &xyzmin[0], &xyzmin[1], &xyzmin[2], - &dx[0], &dx[1], &dx[2], &nx, &ny, &nz, - &nxg, &nyg, &nzg, &WarpX::nox,&WarpX::noy,&WarpX::noz, - &lvect, &WarpX::charge_deposition_algo); -#ifdef WARPX_DIM_RZ - long ngRho = WarpX::nox; - warpx_charge_deposition_rz_volume_scaling( - data_ptr, &ngRho, rholen.getVect(), - &xyzmin[0], &dx[0]); -#endif - -#ifdef _OPENMP - rhofab.atomicAdd(rho_loc); + DepositCharge(pti, wp, rho.get(), 0, 0, np, thread_num, lev, lev); } -#endif } if (!local) rho->SumBoundary(gm.periodicity()); -- cgit v1.2.3 From 9fa3a16002c585720e1741b5fdbe5e6d5b974af9 Mon Sep 17 00:00:00 2001 From: Dave Grote Date: Thu, 8 Aug 2019 16:27:54 -0700 Subject: Implemented ApplyInverseVolumeScalingToChargeDensity --- Source/Evolve/WarpXEvolveEM.cpp | 10 +++++- Source/FieldSolver/WarpXPushFieldsEM.cpp | 49 +++++++++++++++++++++++++++++ Source/Particles/WarpXParticleContainer.cpp | 5 --- Source/WarpX.H | 4 +++ 4 files changed, 62 insertions(+), 6 deletions(-) (limited to 'Source/Particles/WarpXParticleContainer.cpp') diff --git a/Source/Evolve/WarpXEvolveEM.cpp b/Source/Evolve/WarpXEvolveEM.cpp index a097d8814..e05fde03e 100644 --- a/Source/Evolve/WarpXEvolveEM.cpp +++ b/Source/Evolve/WarpXEvolveEM.cpp @@ -483,11 +483,19 @@ WarpX::PushParticlesandDepose (int lev, Real cur_time) Bfield_cax[lev][0].get(), Bfield_cax[lev][1].get(), Bfield_cax[lev][2].get(), cur_time, dt[lev]); #ifdef WARPX_DIM_RZ - // This is called after all particles have deposited their current. + // This is called after all particles have deposited their current and charge. ApplyInverseVolumeScalingToCurrentDensity(current_fp[lev][0].get(), current_fp[lev][1].get(), current_fp[lev][2].get(), lev); if (current_buf[lev][0].get()) { ApplyInverseVolumeScalingToCurrentDensity(current_buf[lev][0].get(), current_buf[lev][1].get(), current_buf[lev][2].get(), lev-1); } + if (rho_fp[lev].get()) { + for (int icomp = 0 ; icomp < rho_fp[lev].get()->nComp() ; icomp++) { + ApplyInverseVolumeScalingToChargeDensity(rho_fp[lev].get(), icomp, lev); + if (charge_buf[lev].get()) { + ApplyInverseVolumeScalingToChargeDensity(charge_buf[lev].get(), icomp, lev-1); + } + } + } #endif } diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index 21000424d..f199e0660 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -671,4 +671,53 @@ WarpX::ApplyInverseVolumeScalingToCurrentDensity (MultiFab* Jx, MultiFab* Jy, Mu }); } } + +void +WarpX::ApplyInverseVolumeScalingToChargeDensity (MultiFab* Rho, int icomp, int lev) +{ + const long ngRho = Rho->nGrow(); + const std::array& dx = WarpX::CellSize(lev); + const Real dr = dx[0]; + + Box tilebox; + + for ( MFIter mfi(*Rho, TilingIfNotGPU()); mfi.isValid(); ++mfi ) + { + + Array4 const& Rho_arr = Rho->array(mfi); + + tilebox = mfi.tilebox(); + Box tb = convert(tilebox, IntVect::TheUnitVector()); + + // Lower corner of tile box physical domain + // Note that this is done before the tilebox.grow so that + // these do not include the guard cells. + const std::array& xyzmin = WarpX::LowerCorner(tilebox, lev); + const Dim3 lo = lbound(tilebox); + const Real rmin = xyzmin[0]; + const int irmin = lo.x; + + // Rescale charge in r-z mode since the inverse volume factor was not + // included in the charge deposition. + amrex::ParallelFor(tb, + [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + // Wrap the charge density deposited in the guard cells around + // to the cells above the axis. + // Rho is located on the boundary + if (rmin == 0. && 0 < i && i <= ngRho) { + Rho_arr(i,j,0,icomp) += Rho_arr(-i,j,0,icomp); + } + + // Apply the inverse volume scaling + const amrex::Real r = std::abs(rmin + (i - irmin)*dr); + if (r == 0.) { + // Verboncoeur JCP 164, 421-427 (2001) : corrected volume on axis + Rho_arr(i,j,0,icomp) /= (MathConst::pi*dr/3.); + } else { + Rho_arr(i,j,0,icomp) /= (2.*MathConst::pi*r); + } + }); + } +} #endif diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index ec0f78564..32abae3ad 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -623,11 +623,6 @@ WarpXParticleContainer::DepositCharge (WarpXParIter& pti, RealVector& wp, doChargeDepositionShapeN<3>(xp, yp, zp, wp.dataPtr()+offset, rho_arr, np_to_depose, dx, xyzmin, lo, q); } -#ifdef WARPX_DIM_RZ - warpx_charge_deposition_rz_volume_scaling( - data_ptr, &ngRho, rholen.getVect(), - &xyzmin[0], &dx[0]); -#endif BL_PROFILE_VAR_STOP(blp_ppc_chd); #ifndef AMREX_USE_GPU diff --git a/Source/WarpX.H b/Source/WarpX.H index d8b1b9ce0..67d2366d4 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -183,6 +183,10 @@ public: amrex::MultiFab* Jy, amrex::MultiFab* Jz, int lev); + + void ApplyInverseVolumeScalingToChargeDensity(amrex::MultiFab* Rho, + int icomp, + int lev); #endif void DampPML (); -- cgit v1.2.3 From c45bde5edf520cd808a620e243e426ace445441d Mon Sep 17 00:00:00 2001 From: Dave Grote Date: Thu, 8 Aug 2019 17:07:04 -0700 Subject: Minor clean up in charge deposition conversion --- Source/Particles/WarpXParticleContainer.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) (limited to 'Source/Particles/WarpXParticleContainer.cpp') diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 32abae3ad..4ccddaedf 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -567,7 +567,7 @@ WarpXParticleContainer::DepositCharge (WarpXParIter& pti, RealVector& wp, const long ngRho = rho->nGrow(); const std::array& dx = WarpX::CellSize(std::max(depos_lev,0)); - Real q = this->charge; + const Real q = this->charge; BL_PROFILE_VAR_NS("PPC::ChargeDeposition", blp_ppc_chd); BL_PROFILE_VAR_NS("PPC::Evolve::Accumulate", blp_accumulate); @@ -586,11 +586,11 @@ WarpXParticleContainer::DepositCharge (WarpXParIter& pti, RealVector& wp, tilebox.grow(ngRho); #ifdef AMREX_USE_GPU - // No tiling on GPU: rho_ptr points to the full rho array. + // No tiling on GPU: rho_arr points to the full rho array. Array4 const& rho_arr = rho->array(pti); #else - // Tiling is on: rho_ptr points to local_rho[thread_num] - Box tb = amrex::convert(tilebox, IntVect::TheUnitVector()); + // Tiling is on: rho_arr points to local_rho[thread_num] + const Box tb = amrex::convert(tilebox, IntVect::TheUnitVector()); local_rho[thread_num].resize(tb); -- cgit v1.2.3 From 00cd6d4cee3c5f01b5e68b080019fcd85dc73714 Mon Sep 17 00:00:00 2001 From: Dave Grote Date: Thu, 8 Aug 2019 17:57:58 -0700 Subject: In GetChargeDensity, added call to ApplyInverseVolumeScalingToChargeDensity --- Source/Particles/WarpXParticleContainer.cpp | 4 ++++ 1 file changed, 4 insertions(+) (limited to 'Source/Particles/WarpXParticleContainer.cpp') diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 4ccddaedf..c4e670bf8 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -735,6 +735,10 @@ WarpXParticleContainer::GetChargeDensity (int lev, bool local) } } +#ifdef WARPX_DIM_RZ + WarpX::GetInstance().ApplyInverseVolumeScalingToChargeDensity(rho.get(), lev); +#endif + if (!local) rho->SumBoundary(gm.periodicity()); return rho; -- cgit v1.2.3 From 8b455ff9c46aa8d8a378d8c262eea8edcdf10e3f Mon Sep 17 00:00:00 2001 From: Dave Grote Date: Fri, 9 Aug 2019 16:37:29 -0700 Subject: Fix charge deposition for GPU --- Source/Particles/WarpXParticleContainer.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) (limited to 'Source/Particles/WarpXParticleContainer.cpp') diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index c4e670bf8..aa74f9c4a 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -587,7 +587,8 @@ WarpXParticleContainer::DepositCharge (WarpXParIter& pti, RealVector& wp, #ifdef AMREX_USE_GPU // No tiling on GPU: rho_arr points to the full rho array. - Array4 const& rho_arr = rho->array(pti); + MultiFab rhoi(*rho, amrex::make_alias, icomp, 1); + Array4 const& rho_arr = rhoi.array(pti); #else // Tiling is on: rho_arr points to local_rho[thread_num] const Box tb = amrex::convert(tilebox, IntVect::TheUnitVector()); -- cgit v1.2.3 From a38e581301c7be825a6b7ca0d742fd714aef882d Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Fri, 9 Aug 2019 20:21:33 -0700 Subject: Correct syntax error in the absence of OPENMP --- Source/Particles/WarpXParticleContainer.cpp | 2 ++ 1 file changed, 2 insertions(+) (limited to 'Source/Particles/WarpXParticleContainer.cpp') diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index aa74f9c4a..3bbb65d1b 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -734,7 +734,9 @@ WarpXParticleContainer::GetChargeDensity (int lev, bool local) DepositCharge(pti, wp, rho.get(), 0, 0, np, thread_num, lev, lev); } +#ifdef _OPENMP } +#endif #ifdef WARPX_DIM_RZ WarpX::GetInstance().ApplyInverseVolumeScalingToChargeDensity(rho.get(), lev); -- cgit v1.2.3