Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Clean up charge deposition code duplication #4113

Merged
merged 2 commits into from
Aug 2, 2023
Merged

Clean up charge deposition code duplication #4113

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
439e09c
clean up charge deposition code duplication
roelof-groenewald Jul 24, 2023
f32a33f
Fix for clang-tidy test
roelof-groenewald Aug 2, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 4 additions & 0 deletions Source/Particles/WarpXParticleContainer.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -207,6 +207,10 @@ public:
const bool apply_boundary_and_scale_volume = false,
const bool interpolate_across_levels = true,
const int icomp = 0);
void DepositCharge (std::unique_ptr<amrex::MultiFab>& rho, const int lev,
const bool local = false, const bool reset = false,
const bool apply_boundary_and_scale_volume = false,
const int icomp = 0);

std::unique_ptr<amrex::MultiFab> GetChargeDensity(int lev, bool local = false);

Expand Down
153 changes: 61 additions & 92 deletions Source/Particles/WarpXParticleContainer.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -976,58 +976,9 @@ WarpXParticleContainer::DepositCharge (amrex::Vector<std::unique_ptr<amrex::Mult
int const finest_level = rho.size() - 1;
for (int lev = 0; lev <= finest_level; ++lev)
{
// Reset the rho array if reset is True
int const nc = WarpX::ncomps;
if (reset) rho[lev]->setVal(0., icomp*nc, nc, rho[lev]->nGrowVect());

// Loop over particle tiles and deposit charge on each level
#ifdef AMREX_USE_OMP
#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
{
const int thread_num = omp_get_thread_num();
#else
const int thread_num = 0;
#endif
for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti)
{
const long np = pti.numParticles();
auto const & wp = pti.GetAttribs(PIdx::w);

int* AMREX_RESTRICT ion_lev = nullptr;
if (do_field_ionization)
{
ion_lev = pti.GetiAttribs(particle_icomps["ionizationLevel"]).dataPtr();
}

DepositCharge(pti, wp, ion_lev, rho[lev].get(), icomp, 0, np, thread_num, lev, lev);
}
#ifdef AMREX_USE_OMP
}
#endif

#ifdef WARPX_DIM_RZ
if (apply_boundary_and_scale_volume)
{
WarpX::GetInstance().ApplyInverseVolumeScalingToChargeDensity(rho[lev].get(), lev);
}
#endif

// Exchange guard cells
if (!local) {
// Possible performance optimization:
// pass less than `rho[lev]->nGrowVect()` in the fifth input variable `dst_ng`
ablastr::utils::communication::SumBoundary(
*rho[lev], 0, rho[lev]->nComp(), rho[lev]->nGrowVect(), rho[lev]->nGrowVect(),
WarpX::do_single_precision_comms, m_gdb->Geom(lev).periodicity());
}

#ifndef WARPX_DIM_RZ
if (apply_boundary_and_scale_volume)
{
// Reflect density over PEC boundaries, if needed.
WarpX::GetInstance().ApplyRhofieldBoundary(lev, rho[lev].get(), PatchType::fine);
}
#endif
DepositCharge (
rho[lev], lev, local, reset, apply_boundary_and_scale_volume, icomp
);
}

// Now that the charge has been deposited at each level,
Expand All @@ -1053,72 +1004,90 @@ WarpXParticleContainer::DepositCharge (amrex::Vector<std::unique_ptr<amrex::Mult
}
}

std::unique_ptr<MultiFab>
WarpXParticleContainer::GetChargeDensity (int lev, bool local)
void
WarpXParticleContainer::DepositCharge (std::unique_ptr<amrex::MultiFab>& rho,
const int lev, const bool local, const bool reset,
const bool apply_boundary_and_scale_volume,
const int icomp)
{
const auto& gm = m_gdb->Geom(lev);
const auto& ba = m_gdb->ParticleBoxArray(lev);
const auto& dm = m_gdb->DistributionMap(lev);
BoxArray nba = ba;

#ifdef WARPX_DIM_RZ
const bool is_PSATD_RZ =
(WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::PSATD);
#else
const bool is_PSATD_RZ = false;
#endif
if( !is_PSATD_RZ )
nba.surroundingNodes();

// Number of guard cells for local deposition of rho
const WarpX& warpx = WarpX::GetInstance();
const int ng_rho = warpx.get_ng_depos_rho().max();

auto rho = std::make_unique<MultiFab>(nba,dm,WarpX::ncomps,ng_rho);
rho->setVal(0.0);
// Reset the rho array if reset is True
int const nc = WarpX::ncomps;
if (reset) rho->setVal(0., icomp*nc, nc, rho->nGrowVect());

// Loop over particle tiles and deposit charge on each level
#ifdef AMREX_USE_OMP
#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
{
#endif
#ifdef AMREX_USE_OMP
const int thread_num = omp_get_thread_num();
const int thread_num = omp_get_thread_num();
#else
const int thread_num = 0;
const int thread_num = 0;
#endif
for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti)
{
const long np = pti.numParticles();
auto const & wp = pti.GetAttribs(PIdx::w);

for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti)
int* AMREX_RESTRICT ion_lev = nullptr;
if (do_field_ionization)
{
const long np = pti.numParticles();
auto& wp = pti.GetAttribs(PIdx::w);

const int* const AMREX_RESTRICT ion_lev = (do_field_ionization)?
pti.GetiAttribs(particle_icomps["ionizationLevel"]).dataPtr():nullptr;

DepositCharge(pti, wp, ion_lev, rho.get(), 0, 0, np,
thread_num, lev, lev);
ion_lev = pti.GetiAttribs(particle_icomps["ionizationLevel"]).dataPtr();
}

DepositCharge(pti, wp, ion_lev, rho.get(), icomp, 0, np, thread_num, lev, lev);
}
#ifdef AMREX_USE_OMP
}
#endif

#ifdef WARPX_DIM_RZ
WarpX::GetInstance().ApplyInverseVolumeScalingToChargeDensity(rho.get(), lev);
if (apply_boundary_and_scale_volume)
{
WarpX::GetInstance().ApplyInverseVolumeScalingToChargeDensity(rho.get(), lev);
}
#endif

if (!local) {
// Exchange guard cells
if (local == false) {
// Possible performance optimization:
// pass less than `rho->nGrowVect()` in the fifth input variable `dst_ng`
ablastr::utils::communication::SumBoundary(
*rho, 0, rho->nComp(), rho->nGrowVect(), rho->nGrowVect(),
WarpX::do_single_precision_comms, gm.periodicity());
WarpX::do_single_precision_comms,
m_gdb->Geom(lev).periodicity()
);
}

#ifndef WARPX_DIM_RZ
// Reflect density over PEC boundaries, if needed.
WarpX::GetInstance().ApplyRhofieldBoundary(lev, rho.get(), PatchType::fine);
if (apply_boundary_and_scale_volume)
{
// Reflect density over PEC boundaries, if needed.
WarpX::GetInstance().ApplyRhofieldBoundary(lev, rho.get(), PatchType::fine);
}
#endif
}

std::unique_ptr<MultiFab>
WarpXParticleContainer::GetChargeDensity (int lev, bool local)
{
const auto& ba = m_gdb->ParticleBoxArray(lev);
const auto& dm = m_gdb->DistributionMap(lev);
BoxArray nba = ba;

#ifdef WARPX_DIM_RZ
const bool is_PSATD_RZ =
(WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::PSATD);
#else
const bool is_PSATD_RZ = false;
#endif
if( !is_PSATD_RZ )
nba.surroundingNodes();

// Number of guard cells for local deposition of rho
const WarpX& warpx = WarpX::GetInstance();
const int ng_rho = warpx.get_ng_depos_rho().max();

auto rho = std::make_unique<MultiFab>(nba, dm, WarpX::ncomps,ng_rho);
DepositCharge(rho, lev, local, true, true, 0);
return rho;
}

Expand Down