diff options
Diffstat (limited to 'Source/Particles/WarpXParticleContainer.cpp')
| -rw-r--r-- | Source/Particles/WarpXParticleContainer.cpp | 142 |
1 files changed, 75 insertions, 67 deletions
diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 7fb57500d..9f02da338 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -1,13 +1,13 @@ - #include <limits> #include <MultiParticleContainer.H> #include <WarpXParticleContainer.H> #include <AMReX_AmrParGDB.H> +#include <WarpXComm.H> #include <WarpX_f.H> #include <WarpX.H> #include <WarpXAlgorithmSelection.H> - +#include <WarpXComm.H> // Import low-level single-particle kernels #include <GetAndSetPosition.H> #include <UpdatePosition.H> @@ -17,6 +17,7 @@ using namespace amrex; int WarpXParticleContainer::do_not_push = 0; +int WarpXParticleContainer::do_not_deposit = 0; WarpXParIter::WarpXParIter (ContainerType& pc, int level) : ParIter(pc, level, MFItInfo().SetDynamic(WarpX::do_dynamic_scheduling)) @@ -270,11 +271,12 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, // If no particles, do not do anything if (np_to_depose == 0) return; + // If user decides not to deposit + if (do_not_deposit) return; + const long ngJ = jx->nGrow(); const std::array<Real,3>& dx = WarpX::CellSize(std::max(depos_lev,0)); - int j_is_nodal = jx->is_nodal() and jy->is_nodal() and jz->is_nodal(); Real q = this->charge; - 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); @@ -300,6 +302,9 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, #ifdef AMREX_USE_GPU // No tiling on GPU: jx_ptr points to the full // jx array (same for jy_ptr and jz_ptr). + auto & jx_fab = jx->get(pti); + auto & jy_fab = jy->get(pti); + auto & jz_fab = jz->get(pti); Array4<Real> const& jx_arr = jx->array(pti); Array4<Real> const& jy_arr = jy->array(pti); Array4<Real> const& jz_arr = jz->array(pti); @@ -319,6 +324,9 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, local_jy[thread_num].setVal(0.0); local_jz[thread_num].setVal(0.0); + auto & jx_fab = local_jx[thread_num]; + auto & jy_fab = local_jy[thread_num]; + auto & jz_fab = local_jz[thread_num]; Array4<Real> const& jx_arr = local_jx[thread_num].array(); Array4<Real> const& jy_arr = local_jy[thread_num].array(); Array4<Real> const& jz_arr = local_jz[thread_num].array(); @@ -333,13 +341,8 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, // Lower corner of tile box physical domain // Note that this includes guard cells since it is after tilebox.ngrow - const std::array<Real, 3>& xyzmin = WarpX::LowerCorner(tilebox, depos_lev); - // xyzmin is built on pti.tilebox(), so it does - // not include staggering, so the stagger_shift has to be done by hand. - // Alternatively, we could define xyzminx from tbx (and the same for 3 - // directions and for jx, jy, jz). This way, sx0 would not be needed. - // Better for memory? worth trying? const Dim3 lo = lbound(tilebox); + const std::array<Real, 3>& xyzmin = WarpX::LowerCorner(tilebox, depos_lev); BL_PROFILE_VAR_START(blp_deposit); if (WarpX::current_deposition_algo == CurrentDepositionAlgo::Esirkepov) { @@ -367,20 +370,20 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, 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); + jx_fab, jy_fab, jz_fab, np_to_depose, dt, dx, + xyzmin, lo, q); } else if (WarpX::nox == 2){ 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); + jx_fab, jy_fab, jz_fab, np_to_depose, dt, dx, + xyzmin, lo, q); } else if (WarpX::nox == 3){ 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); + jx_fab, jy_fab, jz_fab, np_to_depose, dt, dx, + xyzmin, lo, q); } } BL_PROFILE_VAR_STOP(blp_deposit); @@ -428,6 +431,9 @@ WarpXParticleContainer::DepositCharge (WarpXParIter& pti, RealVector& wp, // If no particles, do not do anything if (np_to_depose == 0) return; + // If user decides not to deposit + if (do_not_deposit) return; + const long ngRho = rho->nGrow(); const std::array<Real,3>& dx = WarpX::CellSize(std::max(depos_lev,0)); const Real q = this->charge; @@ -501,68 +507,68 @@ WarpXParticleContainer::DepositCharge (WarpXParIter& pti, RealVector& wp, } void -WarpXParticleContainer::DepositCharge (Vector<std::unique_ptr<MultiFab> >& rho, bool local) +WarpXParticleContainer::DepositCharge (Vector<std::unique_ptr<MultiFab> >& rho, + bool local, bool reset, + bool do_rz_volume_scaling) { + // Loop over the refinement levels + int const finest_level = rho.size() - 1; + for (int lev = 0; lev <= finest_level; ++lev) { - int num_levels = rho.size(); - int finest_level = num_levels - 1; - - // each level deposits it's own particles - const int ng = rho[0]->nGrow(); - for (int lev = 0; lev < num_levels; ++lev) { - - rho[lev]->setVal(0.0, ng); - - const auto& gm = m_gdb->Geom(lev); - const auto& ba = m_gdb->ParticleBoxArray(lev); + // Reset the `rho` array if `reset` is True + if (reset) rho[lev]->setVal(0.0, rho[lev]->nGrow()); - const Real* dx = gm.CellSize(); - const Real* plo = gm.ProbLo(); - BoxArray nba = ba; - nba.surroundingNodes(); + // Loop over particle tiles and deposit charge on each level +#ifdef _OPENMP + #pragma omp parallel + { + 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); - for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) { - const Box& box = nba[pti]; + pti.GetPosition(m_xp[thread_num], m_yp[thread_num], m_zp[thread_num]); - auto& wp = pti.GetAttribs(PIdx::w); - const auto& particles = pti.GetArrayOfStructs(); - int nstride = particles.dataShape().first; - const long np = pti.numParticles(); + int* AMREX_RESTRICT ion_lev; + if (do_field_ionization){ + ion_lev = pti.GetiAttribs(particle_icomps["ionization_level"]).dataPtr(); + } else { + ion_lev = nullptr; + } - FArrayBox& rhofab = (*rho[lev])[pti]; + DepositCharge(pti, wp, ion_lev, rho[lev].get(), 0, 0, np, thread_num, lev, lev); + } +#ifdef _OPENMP + } +#endif - WRPX_DEPOSIT_CIC(particles.dataPtr(), nstride, np, - wp.dataPtr(), &this->charge, - rhofab.dataPtr(), box.loVect(), box.hiVect(), - plo, dx, &ng); +#ifdef WARPX_DIM_RZ + if (do_rz_volume_scaling) { + WarpX::GetInstance().ApplyInverseVolumeScalingToChargeDensity(rho[lev].get(), lev); } +#endif - if (!local) rho[lev]->SumBoundary(gm.periodicity()); + // Exchange guard cells + if (!local) rho[lev]->SumBoundary( m_gdb->Geom(lev).periodicity() ); } - // now we average down fine to crse - std::unique_ptr<MultiFab> crse; + // Now that the charge has been deposited at each level, + // we average down from fine to crse for (int lev = finest_level - 1; lev >= 0; --lev) { - const BoxArray& fine_BA = rho[lev+1]->boxArray(); const DistributionMapping& fine_dm = rho[lev+1]->DistributionMap(); - BoxArray coarsened_fine_BA = fine_BA; + BoxArray coarsened_fine_BA = rho[lev+1]->boxArray(); coarsened_fine_BA.coarsen(m_gdb->refRatio(lev)); - 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 - - for (MFIter mfi(coarsened_fine_data); mfi.isValid(); ++mfi) { - const Box& bx = mfi.validbox(); - const Box& crse_box = coarsened_fine_data[mfi].box(); - const Box& fine_box = (*rho[lev+1])[mfi].box(); - WRPX_SUM_FINE_TO_CRSE_NODAL(bx.loVect(), bx.hiVect(), ratio.getVect(), - coarsened_fine_data[mfi].dataPtr(), crse_box.loVect(), crse_box.hiVect(), - (*rho[lev+1])[mfi].dataPtr(), fine_box.loVect(), fine_box.hiVect()); - } + int const refinement_ratio = 2; - rho[lev]->copy(coarsened_fine_data, m_gdb->Geom(lev).periodicity(), FabArrayBase::ADD); + interpolateDensityFineToCoarse( *rho[lev+1], coarsened_fine_data, refinement_ratio ); + rho[lev]->ParallelAdd( coarsened_fine_data, m_gdb->Geom(lev).periodicity() ); } } @@ -624,7 +630,8 @@ Real WarpXParticleContainer::sumParticleCharge(bool local) { amrex::Real total_charge = 0.0; - for (int lev = 0; lev < finestLevel(); ++lev) + const int nLevels = finestLevel(); + for (int lev = 0; lev < nLevels; ++lev) { #ifdef _OPENMP @@ -654,7 +661,8 @@ std::array<Real, 3> WarpXParticleContainer::meanParticleVelocity(bool local) { amrex::Real inv_clight_sq = 1.0/PhysConst::c/PhysConst::c; - for (int lev = 0; lev <= finestLevel(); ++lev) { + const int nLevels = finestLevel(); + for (int lev = 0; lev <= nLevels; ++lev) { #ifdef _OPENMP #pragma omp parallel reduction(+:vx_total, vy_total, vz_total, np_total) @@ -698,7 +706,8 @@ Real WarpXParticleContainer::maxParticleVelocity(bool local) { amrex::ParticleReal max_v = 0.0; - for (int lev = 0; lev <= finestLevel(); ++lev) + const int nLevels = finestLevel(); + for (int lev = 0; lev <= nLevels; ++lev) { #ifdef _OPENMP @@ -724,7 +733,7 @@ WarpXParticleContainer::PushXES (Real dt) { BL_PROFILE("WPC::PushXES()"); - int num_levels = finestLevel() + 1; + const int num_levels = finestLevel() + 1; for (int lev = 0; lev < num_levels; ++lev) { const auto& gm = m_gdb->Geom(lev); @@ -753,7 +762,8 @@ WarpXParticleContainer::PushXES (Real dt) void WarpXParticleContainer::PushX (Real dt) { - for (int lev = 0; lev <= finestLevel(); ++lev) { + const int nLevels = finestLevel(); + for (int lev = 0; lev <= nLevels; ++lev) { PushX(lev, dt); } } @@ -841,5 +851,3 @@ WarpXParticleContainer::particlePostLocate(ParticleType& p, if (pld.m_lev == lev-1){ } } - - |