From 94c7b75d8d8651dbf0260b7f5aab217991713578 Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Fri, 26 Jul 2019 13:20:18 -0700 Subject: use WarpX current deposition for laser particles as well --- Source/Laser/LaserParticleContainer.cpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) (limited to 'Source/Laser/LaserParticleContainer.cpp') diff --git a/Source/Laser/LaserParticleContainer.cpp b/Source/Laser/LaserParticleContainer.cpp index 3d3447a3c..515aa1f5d 100644 --- a/Source/Laser/LaserParticleContainer.cpp +++ b/Source/Laser/LaserParticleContainer.cpp @@ -504,15 +504,15 @@ LaserParticleContainer::Evolve (int lev, // Current Deposition // // Deposit inside domains - DepositCurrentFortran(pti, wp, uxp, uyp, uzp, &jx, &jy, &jz, - 0, np_current, thread_num, - lev, lev, dt); + DepositCurrent(pti, wp, uxp, uyp, uzp, &jx, &jy, &jz, + 0, np_current, thread_num, + lev, lev, dt); bool has_buffer = cjx; if (has_buffer){ // Deposit in buffers - DepositCurrentFortran(pti, wp, uxp, uyp, uzp, cjx, cjy, cjz, - np_current, np-np_current, thread_num, - lev, lev-1, dt); + DepositCurrent(pti, wp, uxp, uyp, uzp, cjx, cjy, cjz, + np_current, np-np_current, thread_num, + lev, lev-1, dt); } // -- cgit v1.2.3 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/Laser/LaserParticleContainer.cpp | 14 +- Source/Particles/Deposition/Make.package | 1 + Source/Particles/PhysicalParticleContainer.cpp | 14 +- Source/Particles/WarpXParticleContainer.H | 10 +- Source/Particles/WarpXParticleContainer.cpp | 243 +++++++------------------ 5 files changed, 100 insertions(+), 182 deletions(-) (limited to 'Source/Laser/LaserParticleContainer.cpp') diff --git a/Source/Laser/LaserParticleContainer.cpp b/Source/Laser/LaserParticleContainer.cpp index 515aa1f5d..786ebc622 100644 --- a/Source/Laser/LaserParticleContainer.cpp +++ b/Source/Laser/LaserParticleContainer.cpp @@ -453,7 +453,12 @@ LaserParticleContainer::Evolve (int lev, pti.GetPosition(m_xp[thread_num], m_yp[thread_num], m_zp[thread_num]); BL_PROFILE_VAR_STOP(blp_copy); - if (rho) DepositCharge(pti, wp, rho, crho, 0, np_current, np, thread_num, lev); + if (rho) { + DepositCharge(pti, wp, rho, 0, 0, np_current, thread_num, lev, lev); + if (crho) { + DepositCharge(pti, wp, crho, 0, np_current, np-np_current, thread_num, lev, lev-1); + } + } // // Particle Push @@ -522,7 +527,12 @@ LaserParticleContainer::Evolve (int lev, pti.SetPosition(m_xp[thread_num], m_yp[thread_num], m_zp[thread_num]); BL_PROFILE_VAR_STOP(blp_copy); - if (rho) DepositCharge(pti, wp, rho, crho, 1, np_current, np, thread_num, lev); + if (rho) { + DepositCharge(pti, wp, rho, 1, 0, np_current, thread_num, lev, lev); + if (crho) { + DepositCharge(pti, wp, crho, 1, np_current, np-np_current, thread_num, lev, lev-1); + } + } if (cost) { const Box& tbx = pti.tilebox(); diff --git a/Source/Particles/Deposition/Make.package b/Source/Particles/Deposition/Make.package index 0d5ebe2a7..e1aace998 100644 --- a/Source/Particles/Deposition/Make.package +++ b/Source/Particles/Deposition/Make.package @@ -1,3 +1,4 @@ CEXE_headers += CurrentDeposition.H +CEXE_headers += ChargeDeposition.H INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Particles/Deposition VPATH_LOCATIONS += $(WARPX_HOME)/Source/Particles/Deposition diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 59d3ce76f..d10390204 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1156,7 +1156,12 @@ PhysicalParticleContainer::Evolve (int lev, pti.GetPosition(m_xp[thread_num], m_yp[thread_num], m_zp[thread_num]); BL_PROFILE_VAR_STOP(blp_copy); - if (rho) DepositCharge(pti, wp, rho, crho, 0, np_current, np, thread_num, lev); + if (rho) { + DepositCharge(pti, wp, rho, 0, 0, np_current, thread_num, lev, lev); + if (has_buffer){ + DepositCharge(pti, wp, crho, 0, np_current, np-np_current, thread_num, lev, lev-1); + } + } if (! do_not_push) { @@ -1292,7 +1297,12 @@ PhysicalParticleContainer::Evolve (int lev, BL_PROFILE_VAR_STOP(blp_copy); } - if (rho) DepositCharge(pti, wp, rho, crho, 1, np_current, np, thread_num, lev); + if (rho) { + DepositCharge(pti, wp, rho, 1, 0, np_current, thread_num, lev, lev); + if (has_buffer){ + DepositCharge(pti, wp, crho, 1, np_current, np-np_current, thread_num, lev, lev-1); + } + } if (cost) { const Box& tbx = pti.tilebox(); diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index a609b4cb3..ac5b47ada 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -167,13 +167,13 @@ public: virtual void DepositCharge(WarpXParIter& pti, RealVector& wp, - amrex::MultiFab* rhomf, - amrex::MultiFab* crhomf, + amrex::MultiFab* rho, int icomp, - const long np_current, - const long np, + const long offset, + const long np_to_depose, int thread_num, - int lev ); + int lev, + int depos_lev); virtual void DepositCurrent(WarpXParIter& pti, RealVector& wp, 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