diff options
Diffstat (limited to 'Source/Particles/WarpXParticleContainer.cpp')
| -rw-r--r-- | Source/Particles/WarpXParticleContainer.cpp | 281 |
1 files changed, 90 insertions, 191 deletions
diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index befa5cfed..4e80374d8 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -85,17 +85,6 @@ WarpXParticleContainer::WarpXParticleContainer (AmrCore* amr_core, int ispecies) particle_comps["theta"] = PIdx::theta; #endif - if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) - { - particle_comps["xold"] = PIdx::nattribs; - particle_comps["yold"] = PIdx::nattribs+1; - particle_comps["zold"] = PIdx::nattribs+2; - particle_comps["uxold"] = PIdx::nattribs+3; - particle_comps["uyold"] = PIdx::nattribs+4; - particle_comps["uzold"] = PIdx::nattribs+5; - - } - // Initialize temporary local arrays for charge/current deposition int num_threads = 1; #ifdef _OPENMP @@ -238,12 +227,8 @@ WarpXParticleContainer::AddNParticles (int lev, p.pos(1) = z[i]; #endif - if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) - { - auto& ptile = DefineAndReturnParticleTile(0, 0, 0); - ptile.push_back_real(particle_comps["xold"], x[i]); - ptile.push_back_real(particle_comps["yold"], y[i]); - ptile.push_back_real(particle_comps["zold"], z[i]); + if ( (NumRuntimeRealComps()>0) || (NumRuntimeIntComps()>0) ){ + auto& ptile = DefineAndReturnParticleTile(0, 0, 0); } particle_tile.push_back(p); @@ -256,12 +241,8 @@ WarpXParticleContainer::AddNParticles (int lev, particle_tile.push_back_real(PIdx::uy, vy + ibegin, vy + iend); particle_tile.push_back_real(PIdx::uz, vz + ibegin, vz + iend); - if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) - { - auto& ptile = DefineAndReturnParticleTile(0, 0, 0); - ptile.push_back_real(particle_comps["uxold"], vx + ibegin, vx + iend); - ptile.push_back_real(particle_comps["uyold"], vy + ibegin, vy + iend); - ptile.push_back_real(particle_comps["uzold"], vz + ibegin, vz + iend); + if ( (NumRuntimeRealComps()>0) || (NumRuntimeIntComps()>0) ){ + auto& ptile = DefineAndReturnParticleTile(0, 0, 0); } for (int comp = PIdx::uz+1; comp < PIdx::nattribs; ++comp) @@ -282,136 +263,14 @@ WarpXParticleContainer::AddNParticles (int lev, Redistribute(); } -/* \brief Current Deposition for thread thread_num using PICSAR - * \param pti : Particle iterator - * \param wp : Array of particle weights - * \param uxp uyp uzp : Array of particle - * \param jx jy jz : Full array of current density - * \param offset : Index of first particle for which current is deposited - * \param np_to_depose: Number of particles for which current is deposited. - Particles [offset,offset+np_tp_depose] deposit their - current - * \param thread_num : Thread number (if tiling) - * \param lev : Level of box that contains particles - * \param depos_lev : Level on which particles deposit (if buffers are used) - * \param dt : Time step for particle level - */ -void -WarpXParticleContainer::DepositCurrentFortran(WarpXParIter& pti, - RealVector& wp, RealVector& uxp, - RealVector& uyp, RealVector& uzp, - MultiFab* jx, MultiFab* jy, MultiFab* jz, - const long offset, const long np_to_depose, - int thread_num, int lev, int depos_lev, - Real dt) -{ - AMREX_ALWAYS_ASSERT_WITH_MESSAGE((depos_lev==(lev-1)) || - (depos_lev==(lev )), - "Deposition buffers only work for lev-1"); - // If no particles, do not do anything - if (np_to_depose == 0) return; - - const long ngJ = jx->nGrow(); - const std::array<Real,3>& dx = WarpX::CellSize(std::max(depos_lev,0)); - const long lvect = 8; - int j_is_nodal = jx->is_nodal() and jy->is_nodal() and jz->is_nodal(); - - BL_PROFILE_VAR_NS("PICSAR::CurrentDeposition", blp_pxr_cd); - BL_PROFILE_VAR_NS("PPC::Evolve::Accumulate", blp_accumulate); - - // Get tile box where current is deposited. - // The tile box is different when depositing in the buffers (depos_lev<lev) - // or when depositing inside the level (depos_lev=lev) - Box tilebox; - if (lev == depos_lev) { - tilebox = pti.tilebox(); - } else { - const IntVect& ref_ratio = WarpX::RefRatio(depos_lev); - tilebox = amrex::coarsen(pti.tilebox(),ref_ratio); - } - - // Staggered tile boxes (different in each direction) - Box tbx = convert(tilebox, WarpX::jx_nodal_flag); - Box tby = convert(tilebox, WarpX::jy_nodal_flag); - Box tbz = convert(tilebox, WarpX::jz_nodal_flag); - - // Lower corner of tile box physical domain - const std::array<Real,3>& xyzmin_tile = WarpX::LowerCorner(tilebox, depos_lev); - const std::array<Real, 3>& xyzmin = xyzmin_tile; - -#ifdef AMREX_USE_GPU - // No tiling on GPU: jx_ptr points to the full - // jx array (same for jy_ptr and jz_ptr). - Real* jx_ptr = (*jx)[pti].dataPtr(); - Real* jy_ptr = (*jy)[pti].dataPtr(); - Real* jz_ptr = (*jz)[pti].dataPtr(); - - auto jxntot = (*jx)[pti].length(); - auto jyntot = (*jy)[pti].length(); - auto jzntot = (*jz)[pti].length(); -#else - // Tiling is on: jx_ptr points to local_jx[thread_num] - // (same for jy_ptr and jz_ptr) - tbx.grow(ngJ); - tby.grow(ngJ); - tbz.grow(ngJ); - - local_jx[thread_num].resize(tbx); - local_jy[thread_num].resize(tby); - local_jz[thread_num].resize(tbz); - - Real* jx_ptr = local_jx[thread_num].dataPtr(); - Real* jy_ptr = local_jy[thread_num].dataPtr(); - Real* jz_ptr = local_jz[thread_num].dataPtr(); - - // local_jx[thread_num] is set to zero - 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 - // GPU, no tiling: deposit directly in jx - // CPU, tiling: deposit into local_jx - // (same for jx and jz) - BL_PROFILE_VAR_START(blp_pxr_cd); - warpx_current_deposition( - jx_ptr, &ngJ, jxntot.getVect(), - jy_ptr, &ngJ, jyntot.getVect(), - jz_ptr, &ngJ, jzntot.getVect(), - &np_to_depose, - m_xp[thread_num].dataPtr() + offset, - m_yp[thread_num].dataPtr() + offset, - m_zp[thread_num].dataPtr() + offset, - uxp.dataPtr() + offset, - uyp.dataPtr() + offset, - uzp.dataPtr() + offset, - m_giv[thread_num].dataPtr() + offset, - wp.dataPtr() + offset, &this->charge, - &xyzmin[0], &xyzmin[1], &xyzmin[2], - &dt, &dx[0], &dx[1], &dx[2], - &WarpX::nox,&WarpX::noy,&WarpX::noz, &j_is_nodal, - &lvect,&WarpX::current_deposition_algo); - - BL_PROFILE_VAR_STOP(blp_pxr_cd); - -#ifndef AMREX_USE_GPU - BL_PROFILE_VAR_START(blp_accumulate); - // CPU, tiling: atomicAdd local_jx into jx - // (same for jx and jz) - (*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 -} - /* \brief Current Deposition for thread thread_num * \param pti : Particle iterator * \param wp : Array of particle weights * \param uxp uyp uzp : Array of particle + * \param ion_lev : Pointer to array of particle ionization level. This is + required to have the charge of each macroparticle + since q is a scalar. For non-ionizable species, + ion_lev is a null pointer. * \param jx jy jz : Full array of current density * \param offset : Index of first particle for which current is deposited * \param np_to_depose: Number of particles for which current is deposited. @@ -425,6 +284,7 @@ void WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, RealVector& wp, RealVector& uxp, RealVector& uyp, RealVector& uzp, + const int * const ion_lev, MultiFab* jx, MultiFab* jy, MultiFab* jz, const long offset, const long np_to_depose, int thread_num, int lev, int depos_lev, @@ -443,6 +303,8 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, const Real stagger_shift = j_is_nodal ? 0.0 : 0.5; BL_PROFILE_VAR_NS("PPC::Evolve::Accumulate", blp_accumulate); + BL_PROFILE_VAR_NS("PPC::CurrentDeposition", blp_deposit); + // Get tile box where current is deposited. // The tile box is different when depositing in the buffers (depos_lev<lev) @@ -474,9 +336,9 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, tby.grow(ngJ); tbz.grow(ngJ); - local_jx[thread_num].resize(tbx); - local_jy[thread_num].resize(tby); - local_jz[thread_num].resize(tbz); + local_jx[thread_num].resize(tbx, jx->nComp()); + local_jy[thread_num].resize(tby, jy->nComp()); + local_jz[thread_num].resize(tbz, jz->nComp()); // local_jx[thread_num] is set to zero local_jx[thread_num].setVal(0.0); @@ -505,55 +367,82 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, // Better for memory? worth trying? const Dim3 lo = lbound(tilebox); + BL_PROFILE_VAR_START(blp_deposit); if (WarpX::current_deposition_algo == CurrentDepositionAlgo::Esirkepov) { if (WarpX::nox == 1){ - doEsirkepovDepositionShapeN<1>(xp, yp, zp, wp.dataPtr() + offset, uxp.dataPtr() + offset, - uyp.dataPtr() + offset, uzp.dataPtr() + offset, jx_arr, jy_arr, - jz_arr, np_to_depose, dt, dx, - xyzmin, lo, q); + doEsirkepovDepositionShapeN<1>( + xp, yp, zp, wp.dataPtr() + offset, uxp.dataPtr() + offset, + uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev, + jx_arr, jy_arr, jz_arr, np_to_depose, dt, dx, xyzmin, lo, q, + WarpX::n_rz_azimuthal_modes); } else if (WarpX::nox == 2){ - doEsirkepovDepositionShapeN<2>(xp, yp, zp, wp.dataPtr() + offset, uxp.dataPtr() + offset, - uyp.dataPtr() + offset, uzp.dataPtr() + offset, jx_arr, jy_arr, - jz_arr, np_to_depose, dt, dx, - xyzmin, lo, q); + doEsirkepovDepositionShapeN<2>( + xp, yp, zp, wp.dataPtr() + offset, uxp.dataPtr() + offset, + uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev, + jx_arr, jy_arr, jz_arr, np_to_depose, dt, dx, xyzmin, lo, q, + WarpX::n_rz_azimuthal_modes); } else if (WarpX::nox == 3){ - doEsirkepovDepositionShapeN<3>(xp, yp, zp, wp.dataPtr() + offset, uxp.dataPtr() + offset, - uyp.dataPtr() + offset, uzp.dataPtr() + offset, jx_arr, jy_arr, - jz_arr, np_to_depose, dt, dx, - xyzmin, lo, q); + doEsirkepovDepositionShapeN<3>( + xp, yp, zp, wp.dataPtr() + offset, uxp.dataPtr() + offset, + uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev, + jx_arr, jy_arr, jz_arr, np_to_depose, dt, dx, xyzmin, lo, q, + WarpX::n_rz_azimuthal_modes); } } else { if (WarpX::nox == 1){ - doDepositionShapeN<1>(xp, yp, zp, wp.dataPtr() + offset, uxp.dataPtr() + offset, - uyp.dataPtr() + offset, uzp.dataPtr() + offset, jx_arr, jy_arr, - jz_arr, np_to_depose, dt, dx, - xyzmin, lo, stagger_shift, q); + doDepositionShapeN<1>( + xp, yp, zp, wp.dataPtr() + offset, uxp.dataPtr() + offset, + uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev, + jx_arr, jy_arr, jz_arr, np_to_depose, dt, dx, xyzmin, lo, + stagger_shift, q); } else if (WarpX::nox == 2){ - doDepositionShapeN<2>(xp, yp, zp, wp.dataPtr() + offset, uxp.dataPtr() + offset, - uyp.dataPtr() + offset, uzp.dataPtr() + offset, jx_arr, jy_arr, - jz_arr, np_to_depose, dt, dx, - xyzmin, lo, stagger_shift, q); + doDepositionShapeN<2>( + xp, yp, zp, wp.dataPtr() + offset, uxp.dataPtr() + offset, + uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev, + jx_arr, jy_arr, jz_arr, np_to_depose, dt, dx, xyzmin, lo, + stagger_shift, q); } else if (WarpX::nox == 3){ - doDepositionShapeN<3>(xp, yp, zp, wp.dataPtr() + offset, uxp.dataPtr() + offset, - uyp.dataPtr() + offset, uzp.dataPtr() + offset, jx_arr, jy_arr, - jz_arr, np_to_depose, dt, dx, - xyzmin, lo, stagger_shift, q); + doDepositionShapeN<3>( + xp, yp, zp, wp.dataPtr() + offset, uxp.dataPtr() + offset, + uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev, + jx_arr, jy_arr, jz_arr, np_to_depose, dt, dx, xyzmin, lo, + stagger_shift, q); } } + BL_PROFILE_VAR_STOP(blp_deposit); #ifndef AMREX_USE_GPU BL_PROFILE_VAR_START(blp_accumulate); // CPU, tiling: atomicAdd local_jx into jx // (same for jx and jz) - (*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); + (*jx)[pti].atomicAdd(local_jx[thread_num], tbx, tbx, 0, 0, jx->nComp()); + (*jy)[pti].atomicAdd(local_jy[thread_num], tby, tby, 0, 0, jy->nComp()); + (*jz)[pti].atomicAdd(local_jz[thread_num], tbz, tbz, 0, 0, jz->nComp()); BL_PROFILE_VAR_STOP(blp_accumulate); #endif } +/* \brief Charge Deposition for thread thread_num + * \param pti : Particle iterator + * \param wp : Array of particle weights + * \param ion_lev : Pointer to array of particle ionization level. This is + required to have the charge of each macroparticle + since q is a scalar. For non-ionizable species, + ion_lev is a null pointer. + * \param rho : Full array of charge density + * \param icomp : Component of rho into which charge is deposited. + 0: old value (before particle push). + 1: new value (after particle push). + * \param offset : Index of first particle for which charge is deposited + * \param np_to_depose: Number of particles for which charge is deposited. + Particles [offset,offset+np_tp_depose] deposit charge + * \param thread_num : Thread number (if tiling) + * \param lev : Level of box that contains particles + * \param depos_lev : Level on which particles deposit (if buffers are used) + */ void WarpXParticleContainer::DepositCharge (WarpXParIter& pti, RealVector& wp, + const int * const ion_lev, MultiFab* rho, int icomp, const long offset, const long np_to_depose, int thread_num, int lev, int depos_lev) @@ -585,15 +474,17 @@ WarpXParticleContainer::DepositCharge (WarpXParIter& pti, RealVector& wp, tilebox.grow(ngRho); + const int nc = (rho->nComp() == 1 ? 1 : rho->nComp()/2); + #ifdef AMREX_USE_GPU // No tiling on GPU: rho_arr points to the full rho array. - MultiFab rhoi(*rho, amrex::make_alias, icomp, 1); + MultiFab rhoi(*rho, amrex::make_alias, icomp*nc, nc); Array4<Real> 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()); - local_rho[thread_num].resize(tb); + local_rho[thread_num].resize(tb, nc); // local_rho[thread_num] is set to zero local_rho[thread_num].setVal(0.0); @@ -615,21 +506,21 @@ WarpXParticleContainer::DepositCharge (WarpXParIter& pti, RealVector& wp, 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); + doChargeDepositionShapeN<1>(xp, yp, zp, wp.dataPtr()+offset, ion_lev, + 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); + doChargeDepositionShapeN<2>(xp, yp, zp, wp.dataPtr()+offset, ion_lev, + 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); + doChargeDepositionShapeN<3>(xp, yp, zp, wp.dataPtr()+offset, ion_lev, + rho_arr, np_to_depose, dx, xyzmin, lo, q); } BL_PROFILE_VAR_STOP(blp_ppc_chd); #ifndef AMREX_USE_GPU BL_PROFILE_VAR_START(blp_accumulate); - (*rho)[pti].atomicAdd(local_rho[thread_num], tb, tb, 0, icomp, 1); + (*rho)[pti].atomicAdd(local_rho[thread_num], tb, tb, 0, icomp*nc, nc); BL_PROFILE_VAR_STOP(blp_accumulate); #endif @@ -683,7 +574,7 @@ WarpXParticleContainer::DepositCharge (Vector<std::unique_ptr<MultiFab> >& rho, BoxArray coarsened_fine_BA = fine_BA; coarsened_fine_BA.coarsen(m_gdb->refRatio(lev)); - MultiFab coarsened_fine_data(coarsened_fine_BA, fine_dm, 1, 0); + MultiFab coarsened_fine_data(coarsened_fine_BA, fine_dm, rho[lev+1]->nComp(), 0); coarsened_fine_data.setVal(0.0); IntVect ratio(AMREX_D_DECL(2, 2, 2)); // FIXME @@ -712,7 +603,7 @@ WarpXParticleContainer::GetChargeDensity (int lev, bool local) const int ng = WarpX::nox; - auto rho = std::unique_ptr<MultiFab>(new MultiFab(nba,dm,1,ng)); + auto rho = std::unique_ptr<MultiFab>(new MultiFab(nba,dm,WarpX::ncomps,ng)); rho->setVal(0.0); #ifdef _OPENMP @@ -732,7 +623,15 @@ WarpXParticleContainer::GetChargeDensity (int lev, bool local) pti.GetPosition(m_xp[thread_num], m_yp[thread_num], m_zp[thread_num]); - DepositCharge(pti, wp, rho.get(), 0, 0, np, thread_num, lev, lev); + int* AMREX_RESTRICT ion_lev; + if (do_field_ionization){ + ion_lev = pti.GetiAttribs(particle_icomps["ionization_level"]).dataPtr(); + } else { + ion_lev = nullptr; + } + + DepositCharge(pti, wp, ion_lev, rho.get(), 0, 0, np, + thread_num, lev, lev); } #ifdef _OPENMP } |