From eb07f18b380cb6b7205c05e7828f00a0fa2b8161 Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Tue, 12 Feb 2019 10:42:13 -0800 Subject: re-order tree structure of source code --- Source/Particles/WarpXParticleContainer.cpp | 905 ++++++++++++++++++++++++++++ 1 file changed, 905 insertions(+) create mode 100644 Source/Particles/WarpXParticleContainer.cpp (limited to 'Source/Particles/WarpXParticleContainer.cpp') diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp new file mode 100644 index 000000000..f195d420b --- /dev/null +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -0,0 +1,905 @@ + +#include + +#include +#include +#include +#include +#include + +using namespace amrex; + +int WarpXParticleContainer::do_not_push = 0; + +WarpXParIter::WarpXParIter (ContainerType& pc, int level) + : ParIter(pc, level, MFItInfo().SetDynamic(WarpX::do_dynamic_scheduling)) +{ +} + +#if (AMREX_SPACEDIM == 2) +void +WarpXParIter::GetPosition (Cuda::DeviceVector& x, Cuda::DeviceVector& y, Cuda::DeviceVector& z) const +{ + amrex::ParIter<0,0,PIdx::nattribs>::GetPosition(x, z); + y.resize(x.size(), std::numeric_limits::quiet_NaN()); +} + +void +WarpXParIter::SetPosition (const Cuda::DeviceVector& x, const Cuda::DeviceVector& y, const Cuda::DeviceVector& z) +{ + amrex::ParIter<0,0,PIdx::nattribs>::SetPosition(x, z); +} +#endif + +WarpXParticleContainer::WarpXParticleContainer (AmrCore* amr_core, int ispecies) + : ParticleContainer<0,0,PIdx::nattribs>(amr_core->GetParGDB()) + , species_id(ispecies) +{ + for (unsigned int i = PIdx::Ex; i <= PIdx::Bz; ++i) { + communicate_real_comp[i] = false; // Don't need to communicate E and B. + } + SetParticleSize(); + ReadParameters(); + + // Initialize temporary local arrays for charge/current deposition + int num_threads = 1; + #ifdef _OPENMP + #pragma omp parallel + #pragma omp single + num_threads = omp_get_num_threads(); + #endif + local_rho.resize(num_threads); + local_jx.resize(num_threads); + local_jy.resize(num_threads); + local_jz.resize(num_threads); + m_xp.resize(num_threads); + m_yp.resize(num_threads); + m_zp.resize(num_threads); + m_giv.resize(num_threads); + for (int i = 0; i < num_threads; ++i) + { + local_rho[i].reset(nullptr); + local_jx[i].reset(nullptr); + local_jy[i].reset(nullptr); + local_jz[i].reset(nullptr); + } + +} + +void +WarpXParticleContainer::ReadParameters () +{ + static bool initialized = false; + if (!initialized) + { + ParmParse pp("particles"); + +#ifdef AMREX_USE_GPU + do_tiling = false; // By default, tiling is off on GPU +#else + do_tiling = true; +#endif + pp.query("do_tiling", do_tiling); + pp.query("do_not_push", do_not_push); + + initialized = true; + } +} + +void +WarpXParticleContainer::AllocData () +{ + // have to resize here, not in the constructor because grids have not + // been built when constructor was called. + reserveData(); + resizeData(); +} + +void +WarpXParticleContainer::AddOneParticle (int lev, int grid, int tile, + Real x, Real y, Real z, + const std::array& attribs) +{ + auto& particle_tile = GetParticles(lev)[std::make_pair(grid,tile)]; + AddOneParticle(particle_tile, x, y, z, attribs); +} + +void +WarpXParticleContainer::AddOneParticle (ParticleTileType& particle_tile, + Real x, Real y, Real z, + const std::array& attribs) +{ + ParticleType p; + p.id() = ParticleType::NextID(); + p.cpu() = ParallelDescriptor::MyProc(); +#if (AMREX_SPACEDIM == 3) + p.pos(0) = x; + p.pos(1) = y; + p.pos(2) = z; +#elif (AMREX_SPACEDIM == 2) + p.pos(0) = x; + p.pos(1) = z; +#endif + + particle_tile.push_back(p); + particle_tile.push_back_real(attribs); +} + +void +WarpXParticleContainer::AddNParticles (int lev, + int n, const Real* x, const Real* y, const Real* z, + const Real* vx, const Real* vy, const Real* vz, + int nattr, const Real* attr, int uniqueparticles, int id) +{ + BL_ASSERT(nattr == 1); + const Real* weight = attr; + + int ibegin, iend; + if (uniqueparticles) { + ibegin = 0; + iend = n; + } else { + int myproc = ParallelDescriptor::MyProc(); + int nprocs = ParallelDescriptor::NProcs(); + int navg = n/nprocs; + int nleft = n - navg * nprocs; + if (myproc < nleft) { + ibegin = myproc*(navg+1); + iend = ibegin + navg+1; + } else { + ibegin = myproc*navg + nleft; + iend = ibegin + navg; + } + } + + // Add to grid 0 and tile 0 + // Redistribute() will move them to proper places. + std::pair key {0,0}; + auto& particle_tile = GetParticles(lev)[key]; + + for (int i = ibegin; i < iend; ++i) + { + ParticleType p; + if (id==-1) + { + p.id() = ParticleType::NextID(); + } else { + p.id() = id; + } + p.cpu() = ParallelDescriptor::MyProc(); +#if (AMREX_SPACEDIM == 3) + p.pos(0) = x[i]; + p.pos(1) = y[i]; + p.pos(2) = z[i]; +#elif (AMREX_SPACEDIM == 2) + p.pos(0) = x[i]; + p.pos(1) = z[i]; +#endif + particle_tile.push_back(p); + } + + std::size_t np = iend-ibegin; + + if (np > 0) + { + particle_tile.push_back_real(PIdx::w , weight + ibegin, weight + iend); + particle_tile.push_back_real(PIdx::ux, vx + ibegin, vx + iend); + particle_tile.push_back_real(PIdx::uy, vy + ibegin, vy + iend); + particle_tile.push_back_real(PIdx::uz, vz + ibegin, vz + iend); + + for (int comp = PIdx::uz+1; comp < PIdx::nattribs; ++comp) + { + particle_tile.push_back_real(comp, np, 0.0); + } + } + + Redistribute(); +} + + +void +WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, + RealVector& wp, RealVector& uxp, + RealVector& uyp, RealVector& uzp, + MultiFab& jx, MultiFab& jy, MultiFab& jz, + MultiFab* cjx, MultiFab* cjy, MultiFab* cjz, + const long np_current, const long np, + int thread_num, int lev, Real dt ) +{ + Real *jx_ptr, *jy_ptr, *jz_ptr; + const std::array& xyzmin_tile = WarpX::LowerCorner(pti.tilebox(), lev); + const std::array& dx = WarpX::CellSize(lev); + const std::array& cdx = WarpX::CellSize(std::max(lev-1,0)); + const std::array& xyzmin = xyzmin_tile; + const long lvect = 8; + + BL_PROFILE_VAR_NS("PICSAR::CurrentDeposition", blp_pxr_cd); + BL_PROFILE_VAR_NS("PPC::Evolve::Accumulate", blp_accumulate); + + Box tbx = convert(pti.tilebox(), WarpX::jx_nodal_flag); + Box tby = convert(pti.tilebox(), WarpX::jy_nodal_flag); + Box tbz = convert(pti.tilebox(), WarpX::jz_nodal_flag); + + // WarpX assumes the same number of guard cells for Jx, Jy, Jz + long ngJ = jx.nGrow(); + + // Deposit charge for particles that are not in the current buffers + if (np_current > 0) + { + 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); + + jx_ptr = local_jx[thread_num]->dataPtr(); + jy_ptr = local_jy[thread_num]->dataPtr(); + jz_ptr = local_jz[thread_num]->dataPtr(); + + FArrayBox* local_jx_ptr = local_jx[thread_num].get(); + AMREX_LAUNCH_HOST_DEVICE_LAMBDA(tbx, b, + { + local_jx_ptr->setVal(0.0, b, 0, 1); + }); + + FArrayBox* local_jy_ptr = local_jy[thread_num].get(); + AMREX_LAUNCH_HOST_DEVICE_LAMBDA(tby, b, + { + local_jy_ptr->setVal(0.0, b, 0, 1); + }); + + FArrayBox* local_jz_ptr = local_jz[thread_num].get(); + AMREX_LAUNCH_HOST_DEVICE_LAMBDA(tbz, b, + { + local_jz_ptr->setVal(0.0, b, 0, 1); + }); + + auto jxntot = local_jx[thread_num]->length(); + auto jyntot = local_jy[thread_num]->length(); + auto jzntot = local_jz[thread_num]->length(); + + 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_current, + m_xp[thread_num].dataPtr(), + m_yp[thread_num].dataPtr(), + m_zp[thread_num].dataPtr(), + uxp.dataPtr(), uyp.dataPtr(), uzp.dataPtr(), + m_giv[thread_num].dataPtr(), + wp.dataPtr(), &this->charge, + &xyzmin[0], &xyzmin[1], &xyzmin[2], + &dt, &dx[0], &dx[1], &dx[2], + &WarpX::nox,&WarpX::noy,&WarpX::noz, + &lvect,&WarpX::current_deposition_algo); + BL_PROFILE_VAR_STOP(blp_pxr_cd); + + BL_PROFILE_VAR_START(blp_accumulate); + + FArrayBox const* local_jx_const_ptr = local_jx[thread_num].get(); + FArrayBox* global_jx_ptr = jx.fabPtr(pti); + AMREX_LAUNCH_HOST_DEVICE_LAMBDA(tbx, thread_bx, + { + global_jx_ptr->atomicAdd(*local_jx_const_ptr, thread_bx, thread_bx, 0, 0, 1); + }); + + FArrayBox const* local_jy_const_ptr = local_jy[thread_num].get(); + FArrayBox* global_jy_ptr = jy.fabPtr(pti); + AMREX_LAUNCH_HOST_DEVICE_LAMBDA(tby, thread_bx, + { + global_jy_ptr->atomicAdd(*local_jy_const_ptr, thread_bx, thread_bx, 0, 0, 1); + }); + + FArrayBox const* local_jz_const_ptr = local_jz[thread_num].get(); + FArrayBox* global_jz_ptr = jz.fabPtr(pti); + AMREX_LAUNCH_HOST_DEVICE_LAMBDA(tbz, thread_bx, + { + global_jz_ptr->atomicAdd(*local_jz_const_ptr, thread_bx, thread_bx, 0, 0, 1); + }); + + BL_PROFILE_VAR_STOP(blp_accumulate); + } + + // 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); + + tbx = amrex::convert(ctilebox, WarpX::jx_nodal_flag); + tby = amrex::convert(ctilebox, WarpX::jy_nodal_flag); + tbz = amrex::convert(ctilebox, WarpX::jz_nodal_flag); + 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); + + jx_ptr = local_jx[thread_num]->dataPtr(); + jy_ptr = local_jy[thread_num]->dataPtr(); + jz_ptr = local_jz[thread_num]->dataPtr(); + + FArrayBox* local_jx_ptr = local_jx[thread_num].get(); + AMREX_LAUNCH_HOST_DEVICE_LAMBDA(tbx, b, + { + local_jx_ptr->setVal(0.0, b, 0, 1); + }); + + FArrayBox* local_jy_ptr = local_jy[thread_num].get(); + AMREX_LAUNCH_HOST_DEVICE_LAMBDA(tby, b, + { + local_jy_ptr->setVal(0.0, b, 0, 1); + }); + + FArrayBox* local_jz_ptr = local_jz[thread_num].get(); + AMREX_LAUNCH_HOST_DEVICE_LAMBDA(tbz, b, + { + local_jz_ptr->setVal(0.0, b, 0, 1); + }); + auto jxntot = local_jx[thread_num]->length(); + auto jyntot = local_jy[thread_num]->length(); + auto jzntot = local_jz[thread_num]->length(); + + long ncrse = np - np_current; + 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(), + &ncrse, + m_xp[thread_num].dataPtr() +np_current, + m_yp[thread_num].dataPtr() +np_current, + m_zp[thread_num].dataPtr() +np_current, + uxp.dataPtr()+np_current, + uyp.dataPtr()+np_current, + uzp.dataPtr()+np_current, + m_giv[thread_num].dataPtr()+np_current, + wp.dataPtr()+np_current, &this->charge, + &cxyzmin_tile[0], &cxyzmin_tile[1], &cxyzmin_tile[2], + &dt, &cdx[0], &cdx[1], &cdx[2], + &WarpX::nox,&WarpX::noy,&WarpX::noz, + &lvect,&WarpX::current_deposition_algo); + BL_PROFILE_VAR_STOP(blp_pxr_cd); + + BL_PROFILE_VAR_START(blp_accumulate); + + FArrayBox const* local_jx_const_ptr = local_jx[thread_num].get(); + FArrayBox* global_jx_ptr = cjx->fabPtr(pti); + AMREX_LAUNCH_HOST_DEVICE_LAMBDA(tbx, thread_bx, + { + global_jx_ptr->atomicAdd(*local_jx_const_ptr, thread_bx, thread_bx, 0, 0, 1); + }); + + FArrayBox const* local_jy_const_ptr = local_jy[thread_num].get(); + FArrayBox* global_jy_ptr = cjy->fabPtr(pti); + AMREX_LAUNCH_HOST_DEVICE_LAMBDA(tby, thread_bx, + { + global_jy_ptr->atomicAdd(*local_jy_const_ptr, thread_bx, thread_bx, 0, 0, 1); + }); + + FArrayBox const* local_jz_const_ptr = local_jz[thread_num].get(); + FArrayBox* global_jz_ptr = cjz->fabPtr(pti); + AMREX_LAUNCH_HOST_DEVICE_LAMBDA(tbz, thread_bx, + { + global_jz_ptr->atomicAdd(*local_jz_const_ptr, thread_bx, thread_bx, 0, 0, 1); + }); + + BL_PROFILE_VAR_STOP(blp_accumulate); + } +}; + + +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 ) +{ + + 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; + + long ngRho = rhomf->nGrow(); + Real* data_ptr; + Box tile_box = convert(pti.tilebox(), IntVect::TheUnitVector()); + + const std::array& dx = WarpX::CellSize(lev); + const std::array& cdx = WarpX::CellSize(std::max(lev-1,0)); + + // Deposit charge for particles that are not in the current buffers + if (np_current > 0) + { + const std::array& xyzmin = xyzmin_tile; + tile_box.grow(ngRho); + local_rho[thread_num]->resize(tile_box); + FArrayBox* local_rho_ptr = local_rho[thread_num].get(); + AMREX_LAUNCH_HOST_DEVICE_LAMBDA(tile_box, b, + { + local_rho_ptr->setVal(0.0, b, 0, 1); + }); + + data_ptr = local_rho[thread_num]->dataPtr(); + auto rholen = local_rho[thread_num]->length(); +#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); + BL_PROFILE_VAR_STOP(blp_pxr_chd); + + const int ncomp = 1; + FArrayBox const* local_fab = local_rho[thread_num].get(); + FArrayBox* global_fab = rhomf->fabPtr(pti); + BL_PROFILE_VAR_START(blp_accumulate); + AMREX_LAUNCH_HOST_DEVICE_LAMBDA(tile_box, tbx, + { + global_fab->atomicAdd(*local_fab, tbx, tbx, 0, icomp, ncomp); + }); + BL_PROFILE_VAR_STOP(blp_accumulate); + } + + // 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); + + tile_box = amrex::convert(ctilebox, IntVect::TheUnitVector()); + tile_box.grow(ngRho); + local_rho[thread_num]->resize(tile_box); + FArrayBox* local_rho_ptr = local_rho[thread_num].get(); + AMREX_LAUNCH_HOST_DEVICE_LAMBDA(tile_box, b, + { + local_rho_ptr->setVal(0.0, b, 0, 1); + }); + + data_ptr = local_rho[thread_num]->dataPtr(); + auto rholen = local_rho[thread_num]->length(); +#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 + + 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_STOP(blp_pxr_chd); + + const int ncomp = 1; + FArrayBox const* local_fab = local_rho[thread_num].get(); + FArrayBox* global_fab = crhomf->fabPtr(pti); + BL_PROFILE_VAR_START(blp_accumulate); + AMREX_LAUNCH_HOST_DEVICE_LAMBDA(tile_box, tbx, + { + global_fab->atomicAdd(*local_fab, tbx, tbx, 0, icomp, ncomp); + }); + BL_PROFILE_VAR_STOP(blp_accumulate); + } +}; + +void +WarpXParticleContainer::DepositCharge (Vector >& rho, bool local) +{ + + 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); + const auto& dm = m_gdb->DistributionMap(lev); + + const Real* dx = gm.CellSize(); + const Real* plo = gm.ProbLo(); + BoxArray nba = ba; + nba.surroundingNodes(); + + for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) { + const Box& box = nba[pti]; + + auto& wp = pti.GetAttribs(PIdx::w); + const auto& particles = pti.GetArrayOfStructs(); + int nstride = particles.dataShape().first; + const long np = pti.numParticles(); + + FArrayBox& rhofab = (*rho[lev])[pti]; + + WRPX_DEPOSIT_CIC(particles.dataPtr(), nstride, np, + wp.dataPtr(), &this->charge, + rhofab.dataPtr(), box.loVect(), box.hiVect(), + plo, dx, &ng); + } + + if (!local) rho[lev]->SumBoundary(gm.periodicity()); + } + + // now we average down fine to crse + std::unique_ptr 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; + coarsened_fine_BA.coarsen(m_gdb->refRatio(lev)); + + MultiFab coarsened_fine_data(coarsened_fine_BA, fine_dm, 1, 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()); + } + + rho[lev]->copy(coarsened_fine_data, m_gdb->Geom(lev).periodicity(), FabArrayBase::ADD); + } +} + +std::unique_ptr +WarpXParticleContainer::GetChargeDensity (int lev, bool local) +{ + const auto& gm = m_gdb->Geom(lev); + const auto& ba = m_gdb->ParticleBoxArray(lev); + const auto& dm = m_gdb->DistributionMap(lev); + 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)); + rho->setVal(0.0); + +#ifdef _OPENMP +#pragma omp parallel +#endif + { + Cuda::DeviceVector xp, yp, zp; + FArrayBox local_rho; + + for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) + { + const Box& box = pti.validbox(); + + auto& wp = pti.GetAttribs(PIdx::w); + + const long np = pti.numParticles(); + + pti.GetPosition(xp, yp, zp); + + const std::array& xyzmin_tile = WarpX::LowerCorner(pti.tilebox(), lev); + const std::array& xyzmin_grid = WarpX::LowerCorner(box, lev); + + // Data on the grid + Real* data_ptr; + FArrayBox& rhofab = (*rho)[pti]; +#ifdef _OPENMP + Box tile_box = convert(pti.tilebox(), IntVect::TheUnitVector()); + const std::array& xyzmin = xyzmin_tile; + tile_box.grow(ng); + local_rho.resize(tile_box); + local_rho = 0.0; + data_ptr = local_rho.dataPtr(); + auto rholen = local_rho.length(); +#else + 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 _OPENMP + rhofab.atomicAdd(local_rho); +#endif + } + + } + + if (!local) rho->SumBoundary(gm.periodicity()); + + return rho; +} + +Real WarpXParticleContainer::sumParticleCharge(bool local) { + + amrex::Real total_charge = 0.0; + + for (int lev = 0; lev < finestLevel(); ++lev) + { + +#ifdef _OPENMP +#pragma omp parallel reduction(+:total_charge) +#endif + for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) + { + auto& wp = pti.GetAttribs(PIdx::w); + for (unsigned long i = 0; i < wp.size(); i++) { + total_charge += wp[i]; + } + } + } + + if (!local) ParallelDescriptor::ReduceRealSum(total_charge); + total_charge *= this->charge; + return total_charge; +} + +std::array WarpXParticleContainer::meanParticleVelocity(bool local) { + + amrex::Real vx_total = 0.0; + amrex::Real vy_total = 0.0; + amrex::Real vz_total = 0.0; + + long np_total = 0; + + amrex::Real inv_clight_sq = 1.0/PhysConst::c/PhysConst::c; + + for (int lev = 0; lev <= finestLevel(); ++lev) { + +#ifdef _OPENMP +#pragma omp parallel reduction(+:vx_total, vy_total, vz_total, np_total) +#endif + for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) + { + auto& ux = pti.GetAttribs(PIdx::ux); + auto& uy = pti.GetAttribs(PIdx::uy); + auto& uz = pti.GetAttribs(PIdx::uz); + + np_total += pti.numParticles(); + + for (unsigned long i = 0; i < ux.size(); i++) { + Real usq = (ux[i]*ux[i] + uy[i]*uy[i] + uz[i]*uz[i])*inv_clight_sq; + Real gaminv = 1.0/std::sqrt(1.0 + usq); + vx_total += ux[i]*gaminv; + vy_total += uy[i]*gaminv; + vz_total += uz[i]*gaminv; + } + } + } + + if (!local) { + ParallelDescriptor::ReduceRealSum(vx_total); + ParallelDescriptor::ReduceRealSum(vy_total); + ParallelDescriptor::ReduceRealSum(vz_total); + ParallelDescriptor::ReduceLongSum(np_total); + } + + std::array mean_v; + if (np_total > 0) { + mean_v[0] = vx_total / np_total; + mean_v[1] = vy_total / np_total; + mean_v[2] = vz_total / np_total; + } + + return mean_v; +} + +Real WarpXParticleContainer::maxParticleVelocity(bool local) { + + amrex::Real max_v = 0.0; + + for (int lev = 0; lev <= finestLevel(); ++lev) + { + +#ifdef _OPENMP +#pragma omp parallel reduction(max:max_v) +#endif + for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) + { + auto& ux = pti.GetAttribs(PIdx::ux); + auto& uy = pti.GetAttribs(PIdx::uy); + auto& uz = pti.GetAttribs(PIdx::uz); + for (unsigned long i = 0; i < ux.size(); i++) { + max_v = std::max(max_v, sqrt(ux[i]*ux[i] + uy[i]*uy[i] + uz[i]*uz[i])); + } + } + } + + if (!local) ParallelDescriptor::ReduceRealMax(max_v); + return max_v; +} + +void +WarpXParticleContainer::PushXES (Real dt) +{ + BL_PROFILE("WPC::PushXES()"); + + int num_levels = finestLevel() + 1; + + for (int lev = 0; lev < num_levels; ++lev) { + const auto& gm = m_gdb->Geom(lev); + const RealBox& prob_domain = gm.ProbDomain(); + for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) { + auto& particles = pti.GetArrayOfStructs(); + int nstride = particles.dataShape().first; + const long np = pti.numParticles(); + + auto& attribs = pti.GetAttribs(); + auto& uxp = attribs[PIdx::ux]; + auto& uyp = attribs[PIdx::uy]; + auto& uzp = attribs[PIdx::uz]; + + WRPX_PUSH_LEAPFROG_POSITIONS(particles.dataPtr(), nstride, np, + uxp.dataPtr(), uyp.dataPtr(), +#if AMREX_SPACEDIM == 3 + uzp.dataPtr(), +#endif + &dt, + prob_domain.lo(), prob_domain.hi()); + } + } +} + +void +WarpXParticleContainer::PushX (Real dt) +{ + for (int lev = 0; lev <= finestLevel(); ++lev) { + PushX(lev, dt); + } +} + +void +WarpXParticleContainer::PushX (int lev, Real dt) +{ + BL_PROFILE("WPC::PushX()"); + BL_PROFILE_VAR_NS("WPC::PushX::Copy", blp_copy); + BL_PROFILE_VAR_NS("WPC:PushX::Push", blp_pxr_pp); + + if (do_not_push) return; + + MultiFab* cost = WarpX::getCosts(lev); + +#ifdef _OPENMP +#pragma omp parallel +#endif + { + Cuda::DeviceVector xp, yp, zp, giv; + + for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) + { + Real wt = amrex::second(); + + auto& attribs = pti.GetAttribs(); + + auto& uxp = attribs[PIdx::ux]; + auto& uyp = attribs[PIdx::uy]; + auto& uzp = attribs[PIdx::uz]; + + const long np = pti.numParticles(); + + giv.resize(np); + + // + // copy data from particle container to temp arrays + // + BL_PROFILE_VAR_START(blp_copy); + pti.GetPosition(xp, yp, zp); + BL_PROFILE_VAR_STOP(blp_copy); + + // + // Particle Push + // + BL_PROFILE_VAR_START(blp_pxr_pp); + warpx_particle_pusher_positions(&np, + xp.dataPtr(), + yp.dataPtr(), + zp.dataPtr(), + uxp.dataPtr(), uyp.dataPtr(), uzp.dataPtr(), + giv.dataPtr(), &dt); + BL_PROFILE_VAR_STOP(blp_pxr_pp); + + // + // copy particle data back + // + BL_PROFILE_VAR_START(blp_copy); + pti.SetPosition(xp, yp, zp); + BL_PROFILE_VAR_STOP(blp_copy); + + if (cost) { + const Box& tbx = pti.tilebox(); + wt = (amrex::second() - wt) / tbx.d_numPts(); + FArrayBox* costfab = cost->fabPtr(pti); + AMREX_LAUNCH_HOST_DEVICE_LAMBDA ( tbx, work_box, + { + costfab->plus(wt, work_box); + }); + } + } + } +} + +// This function is called in Redistribute, just after locate +void +WarpXParticleContainer::particlePostLocate(ParticleType& p, + const ParticleLocData& pld, + const int lev) +{ + // Tag particle if goes to higher level. + // It will be split later in the loop + if (pld.m_lev == lev+1 + and p.m_idata.id != NoSplitParticleID + and p.m_idata.id >= 0) + { + p.m_idata.id = DoSplitParticleID; + } + // For the moment, do not do anything if particles goes + // to lower level. + if (pld.m_lev == lev-1){ + } +} -- cgit v1.2.3 From 8282a703e7ab08d1f3f5cf197aa242fabf717c57 Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Fri, 1 Mar 2019 10:25:55 -0800 Subject: move warpxwrappers. rename ParticleContainer files to MultiParticleContainer files --- Source/Diagnostics/BoostedFrameDiagnostic.H | 2 +- Source/Diagnostics/ParticleIO.cpp | 2 +- Source/Laser/LaserParticleContainer.cpp | 2 +- Source/Particles/Make.package | 4 +- Source/Particles/MultiParticleContainer.H | 209 ++++++++++++ Source/Particles/MultiParticleContainer.cpp | 389 +++++++++++++++++++++ Source/Particles/ParticleContainer.H | 209 ------------ Source/Particles/ParticleContainer.cpp | 389 --------------------- Source/Particles/PhysicalParticleContainer.cpp | 2 +- Source/Particles/WarpXParticleContainer.cpp | 2 +- Source/Python/Make.package | 3 +- Source/Python/WarpXWrappers.cpp | 446 +++++++++++++++++++++++++ Source/Python/WarpXWrappers.h | 120 +++++++ Source/Utils/Make.package | 2 - Source/Utils/WarpXConst.cpp | 2 +- Source/Utils/WarpXWrappers.cpp | 446 ------------------------- Source/Utils/WarpXWrappers.h | 120 ------- 17 files changed, 1174 insertions(+), 1175 deletions(-) create mode 100644 Source/Particles/MultiParticleContainer.H create mode 100644 Source/Particles/MultiParticleContainer.cpp delete mode 100644 Source/Particles/ParticleContainer.H delete mode 100644 Source/Particles/ParticleContainer.cpp create mode 100644 Source/Python/WarpXWrappers.cpp create mode 100644 Source/Python/WarpXWrappers.h delete mode 100644 Source/Utils/WarpXWrappers.cpp delete mode 100644 Source/Utils/WarpXWrappers.h (limited to 'Source/Particles/WarpXParticleContainer.cpp') diff --git a/Source/Diagnostics/BoostedFrameDiagnostic.H b/Source/Diagnostics/BoostedFrameDiagnostic.H index 452062137..17f2b161c 100644 --- a/Source/Diagnostics/BoostedFrameDiagnostic.H +++ b/Source/Diagnostics/BoostedFrameDiagnostic.H @@ -7,7 +7,7 @@ #include #include -#include "ParticleContainer.H" +#include "MultiParticleContainer.H" #include "WarpXConst.H" /// diff --git a/Source/Diagnostics/ParticleIO.cpp b/Source/Diagnostics/ParticleIO.cpp index ced2d6b9a..be9809bac 100644 --- a/Source/Diagnostics/ParticleIO.cpp +++ b/Source/Diagnostics/ParticleIO.cpp @@ -1,5 +1,5 @@ -#include +#include #include using namespace amrex; diff --git a/Source/Laser/LaserParticleContainer.cpp b/Source/Laser/LaserParticleContainer.cpp index d5665e053..08d0ae861 100644 --- a/Source/Laser/LaserParticleContainer.cpp +++ b/Source/Laser/LaserParticleContainer.cpp @@ -7,7 +7,7 @@ #include #include #include -#include +#include using namespace amrex; diff --git a/Source/Particles/Make.package b/Source/Particles/Make.package index 93f7ca515..acbfe3390 100644 --- a/Source/Particles/Make.package +++ b/Source/Particles/Make.package @@ -2,10 +2,10 @@ F90EXE_sources += deposit_cic.F90 F90EXE_sources += interpolate_cic.F90 F90EXE_sources += push_particles_ES.F90 CEXE_sources += PhysicalParticleContainer.cpp -CEXE_sources += ParticleContainer.cpp +CEXE_sources += MultiParticleContainer.cpp CEXE_sources += WarpXParticleContainer.cpp CEXE_sources += RigidInjectedParticleContainer.cpp -CEXE_headers += ParticleContainer.H +CEXE_headers += MultiParticleContainer.H CEXE_headers += WarpXParticleContainer.H CEXE_headers += RigidInjectedParticleContainer.H CEXE_headers += PhysicalParticleContainer.H diff --git a/Source/Particles/MultiParticleContainer.H b/Source/Particles/MultiParticleContainer.H new file mode 100644 index 000000000..b21247ff6 --- /dev/null +++ b/Source/Particles/MultiParticleContainer.H @@ -0,0 +1,209 @@ + +#ifndef WARPX_ParticleContainer_H_ +#define WARPX_ParticleContainer_H_ + +#include +#include +#include + +#include + +#include +#include +#include +#include + +// +// MultiParticleContainer holds multiple (nspecies or npsecies+1 when +// using laser) WarpXParticleContainer. WarpXParticleContainer is +// derived from amrex::ParticleContainer, and it is +// polymorphic. PhysicalParticleContainer and LaserParticleContainer are +// concrete class dervied from WarpXParticlContainer. +// + +class MultiParticleContainer +{ + +public: + + MultiParticleContainer (amrex::AmrCore* amr_core); + + ~MultiParticleContainer() {} + + WarpXParticleContainer& GetParticleContainer (int ispecies) { + return *allcontainers[ispecies]; + } + + std::array meanParticleVelocity(int ispecies) { + return allcontainers[ispecies]->meanParticleVelocity(); + } + + void AllocData (); + + void InitData (); + +#ifdef WARPX_DO_ELECTROSTATIC + /// + /// Performs the field gather operation using the input field E, for all the species + /// in the MultiParticleContainer. This is the electrostatic version of the field gather. + /// + void FieldGatherES (const amrex::Vector, 3> >& E, + const amrex::Vector > > >& masks); + + /// + /// This evolves all the particles by one PIC time step, including charge deposition, the + /// field solve, and pushing the particles, for all the species in the MultiParticleContainer. + /// This is the electrostatic version. + /// + void EvolveES (const amrex::Vector, 3> >& E, + amrex::Vector >& rho, + amrex::Real t, amrex::Real dt); + + /// + /// This pushes the particle positions by one half time step for all the species in the + /// MultiParticleContainer. It is used to desynchronize the particles after initializaton + /// or when restarting from a checkpoint. This is the electrostatic version. + /// + void PushXES (amrex::Real dt); + + /// + /// This deposits the particle charge onto rho, accumulating the value for all the species + /// in the MultiParticleContainer. rho is assumed to contain node-centered multifabs. + /// This version is hard-coded for CIC deposition. + /// + void DepositCharge(amrex::Vector >& rho, bool local = false); + + /// + /// This returns the total particle charge for all the species in this MultiParticleContainer. + /// This is needed to subtract the offset for periodic boundary conditions. + /// + amrex::Real sumParticleCharge(bool local = false); +#endif // WARPX_DO_ELECTROSTATIC + + /// + /// Performs the field gather operation using the input fields E and B, for all the species + /// in the MultiParticleContainer. This is the electromagnetic version of the field gather. + /// + void FieldGather (int lev, + const amrex::MultiFab& Ex, const amrex::MultiFab& Ey, const amrex::MultiFab& Ez, + const amrex::MultiFab& Bx, const amrex::MultiFab& By, const amrex::MultiFab& Bz); + + /// + /// This evolves all the particles by one PIC time step, including current deposition, the + /// field solve, and pushing the particles, for all the species in the MultiParticleContainer. + /// This is the electromagnetic version. + /// + void Evolve (int lev, + const amrex::MultiFab& Ex, const amrex::MultiFab& Ey, const amrex::MultiFab& Ez, + const amrex::MultiFab& Bx, const amrex::MultiFab& By, const amrex::MultiFab& Bz, + amrex::MultiFab& jx, amrex::MultiFab& jy, amrex::MultiFab& jz, + amrex::MultiFab* cjx, amrex::MultiFab* cjy, amrex::MultiFab* cjz, + amrex::MultiFab* rho, amrex::MultiFab* crho, + const amrex::MultiFab* cEx, const amrex::MultiFab* cEy, const amrex::MultiFab* cEz, + const amrex::MultiFab* cBx, const amrex::MultiFab* cBy, const amrex::MultiFab* cBz, + amrex::Real t, amrex::Real dt); + + /// + /// This pushes the particle positions by one half time step for all the species in the + /// MultiParticleContainer. It is used to desynchronize the particles after initializaton + /// or when restarting from a checkpoint. This is the electromagnetic version. + /// + void PushX (amrex::Real dt); + + /// + /// This pushes the particle momenta by dt for all the species in the + /// MultiParticleContainer. It is used to desynchronize the particles after initializaton + /// or when restarting from a checkpoint. It is also used to synchronize particles at the + /// the end of the run. This is the electromagnetic version. + /// + void PushP (int lev, amrex::Real dt, + const amrex::MultiFab& Ex, const amrex::MultiFab& Ey, const amrex::MultiFab& Ez, + const amrex::MultiFab& Bx, const amrex::MultiFab& By, const amrex::MultiFab& Bz); + + /// + /// This deposits the particle charge onto a node-centered MultiFab and returns a unique ptr + /// to it. The charge density is accumulated over all the particles in the MultiParticleContainer + /// This version uses PICSAR routines to perform the deposition. + /// + std::unique_ptr GetChargeDensity(int lev, bool local = false); + + void Checkpoint (const std::string& dir, + bool is_checkpoint, + const amrex::Vector& varnames = amrex::Vector()) const; + + void Restart (const std::string& dir); + + void PostRestart (); + + void ReadHeader (std::istream& is); + + void WriteHeader (std::ostream& os) const; + + void SortParticlesByCell (); + + void Redistribute (); + + void RedistributeLocal (const int num_ghost); + + amrex::Vector NumberOfParticlesInGrid(int lev) const; + + void Increment (amrex::MultiFab& mf, int lev); + + void SetParticleBoxArray (int lev, amrex::BoxArray& new_ba); + void SetParticleDistributionMap (int lev, amrex::DistributionMapping& new_dm); + + int nSpecies() const {return nspecies;} + + int nSpeciesDepositOnMainGrid () const { + int r = 0; + for (int i : deposit_on_main_grid) { + if (i) ++r; + } + return r; + } + + void GetLabFrameData(const std::string& snapshot_name, + const int i_lab, const int direction, + const amrex::Real z_old, const amrex::Real z_new, + const amrex::Real t_boost, const amrex::Real t_lab, const amrex::Real dt, + amrex::Vector& parts) const; + + // + // Parameters for the Cherenkov corrector in the FDTD solver. + // Both stencils are calculated ar runtime. + // + // Number of coefficients for the stencil of the NCI corrector. + // The stencil is applied in the z direction only. + static constexpr int nstencilz_fdtd_nci_corr=5; + + amrex::Vector > fdtd_nci_stencilz_ex; + amrex::Vector > fdtd_nci_stencilz_by; + + std::vector GetSpeciesNames() const { return species_names; } + + PhysicalParticleContainer& GetPCtmp () { return *pc_tmp; } + +protected: + + // Particle container types + enum struct PCTypes {Physical, RigidInjected}; + + std::vector species_names; + + std::vector deposit_on_main_grid; + + std::vector species_types; + +private: + + // physical particles (+ laser) + amrex::Vector > allcontainers; + // Temporary particle container, used e.g. for particle splitting. + std::unique_ptr pc_tmp; + + void ReadParameters (); + + // runtime parameters + int nspecies = 1; // physical particles only. If WarpX::use_laser, nspecies+1 == allcontainers.size(). +}; +#endif /*WARPX_ParticleContainer_H_*/ diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp new file mode 100644 index 000000000..1b644b543 --- /dev/null +++ b/Source/Particles/MultiParticleContainer.cpp @@ -0,0 +1,389 @@ +#include +#include +#include + +#include +#include +#include + +using namespace amrex; + +constexpr int MultiParticleContainer::nstencilz_fdtd_nci_corr; + +MultiParticleContainer::MultiParticleContainer (AmrCore* amr_core) +{ + ReadParameters(); + + int n = WarpX::use_laser ? nspecies+1 : nspecies; + allcontainers.resize(n); + for (int i = 0; i < nspecies; ++i) { + if (species_types[i] == PCTypes::Physical) { + allcontainers[i].reset(new PhysicalParticleContainer(amr_core, i, species_names[i])); + } + else if (species_types[i] == PCTypes::RigidInjected) { + allcontainers[i].reset(new RigidInjectedParticleContainer(amr_core, i, species_names[i])); + } + allcontainers[i]->deposit_on_main_grid = deposit_on_main_grid[i]; + } + if (WarpX::use_laser) { + allcontainers[n-1].reset(new LaserParticleContainer(amr_core,n-1)); + } + pc_tmp.reset(new PhysicalParticleContainer(amr_core)); +} + +void +MultiParticleContainer::ReadParameters () +{ + static bool initialized = false; + if (!initialized) + { + ParmParse pp("particles"); + + pp.query("nspecies", nspecies); + BL_ASSERT(nspecies >= 0); + + if (nspecies > 0) { + pp.getarr("species_names", species_names); + BL_ASSERT(species_names.size() == nspecies); + + deposit_on_main_grid.resize(nspecies, 0); + std::vector tmp; + pp.queryarr("deposit_on_main_grid", tmp); + for (auto const& name : tmp) { + auto it = std::find(species_names.begin(), species_names.end(), name); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(it != species_names.end(), "ERROR: species in particles.deposit_on_main_grid must be part of particles.species_names"); + int i = std::distance(species_names.begin(), it); + deposit_on_main_grid[i] = 1; + } + + species_types.resize(nspecies, PCTypes::Physical); + + std::vector rigid_injected_species; + pp.queryarr("rigid_injected_species", rigid_injected_species); + + if (!rigid_injected_species.empty()) { + for (auto const& name : rigid_injected_species) { + auto it = std::find(species_names.begin(), species_names.end(), name); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(it != species_names.end(), "ERROR: species in particles.rigid_injected_species must be part of particles.species_names"); + int i = std::distance(species_names.begin(), it); + species_types[i] = PCTypes::RigidInjected; + } + } + } + pp.query("use_fdtd_nci_corr", WarpX::use_fdtd_nci_corr); + pp.query("l_lower_order_in_v", WarpX::l_lower_order_in_v); + initialized = true; + } +} + +void +MultiParticleContainer::AllocData () +{ + for (auto& pc : allcontainers) { + pc->AllocData(); + } + pc_tmp->AllocData(); +} + +void +MultiParticleContainer::InitData () +{ + for (auto& pc : allcontainers) { + pc->InitData(); + } + pc_tmp->InitData(); +} + + +#ifdef WARPX_DO_ELECTROSTATIC +void +MultiParticleContainer::FieldGatherES (const Vector, 3> >& E, + const amrex::Vector > > >& masks) +{ + for (auto& pc : allcontainers) { + pc->FieldGatherES(E, masks); + } +} + +void +MultiParticleContainer::EvolveES (const Vector, 3> >& E, + Vector >& rho, + Real t, Real dt) +{ + + int nlevs = rho.size(); + int ng = rho[0]->nGrow(); + + for (unsigned i = 0; i < nlevs; i++) { + rho[i]->setVal(0.0, ng); + } + + for (auto& pc : allcontainers) { + pc->EvolveES(E, rho, t, dt); + } + + for (unsigned i = 0; i < nlevs; i++) { + const Geometry& gm = allcontainers[0]->Geom(i); + rho[i]->SumBoundary(gm.periodicity()); + } +} + +void +MultiParticleContainer::Evolve (int lev, + const MultiFab& Ex, const MultiFab& Ey, const MultiFab& Ez, + const MultiFab& Bx, const MultiFab& By, const MultiFab& Bz, + MultiFab& jx, MultiFab& jy, MultiFab& jz, + MultiFab* cjx, MultiFab* cjy, MultiFab* cjz, + MultiFab* rho, + const MultiFab* cEx, const MultiFab* cEy, const MultiFab* cEz, + const MultiFab* cBx, const MultiFab* cBy, const MultiFab* cBz, + Real t, Real dt) +{ + jx.setVal(0.0); + jy.setVal(0.0); + jz.setVal(0.0); + if (cjx) cjx->setVal(0.0); + if (cjy) cjy->setVal(0.0); + if (cjz) cjz->setVal(0.0); + if (rho) rho->setVal(0.0); + for (auto& pc : allcontainers) { + pc->Evolve(lev, Ex, Ey, Ez, Bx, By, Bz, jx, jy, jz, cjx, cjy, cjz, + rho, cEx, cEy, cEz, cBx, cBy, cBz, t, dt); + } +} + +void +MultiParticleContainer::PushXES (Real dt) +{ + for (auto& pc : allcontainers) { + pc->PushXES(dt); + } +} + +void +MultiParticleContainer:: +DepositCharge (Vector >& rho, bool local) +{ + int nlevs = rho.size(); + int ng = rho[0]->nGrow(); + + for (unsigned i = 0; i < nlevs; i++) { + rho[i]->setVal(0.0, ng); + } + + for (unsigned i = 0, n = allcontainers.size(); i < n; ++i) { + allcontainers[i]->DepositCharge(rho, true); + } + + if (!local) { + for (unsigned i = 0; i < nlevs; i++) { + const Geometry& gm = allcontainers[0]->Geom(i); + rho[i]->SumBoundary(gm.periodicity()); + } + } +} + +amrex::Real +MultiParticleContainer::sumParticleCharge (bool local) +{ + amrex::Real total_charge = allcontainers[0]->sumParticleCharge(local); + for (unsigned i = 1, n = allcontainers.size(); i < n; ++i) { + total_charge += allcontainers[i]->sumParticleCharge(local); + } + return total_charge; +} + +#endif // WARPX_DO_ELECTROSTATIC + +void +MultiParticleContainer::FieldGather (int lev, + const MultiFab& Ex, const MultiFab& Ey, const MultiFab& Ez, + const MultiFab& Bx, const MultiFab& By, const MultiFab& Bz) +{ + for (auto& pc : allcontainers) { + pc->FieldGather(lev, Ex, Ey, Ez, Bx, By, Bz); + } +} + +void +MultiParticleContainer::Evolve (int lev, + const MultiFab& Ex, const MultiFab& Ey, const MultiFab& Ez, + const MultiFab& Bx, const MultiFab& By, const MultiFab& Bz, + MultiFab& jx, MultiFab& jy, MultiFab& jz, + MultiFab* cjx, MultiFab* cjy, MultiFab* cjz, + MultiFab* rho, MultiFab* crho, + const MultiFab* cEx, const MultiFab* cEy, const MultiFab* cEz, + const MultiFab* cBx, const MultiFab* cBy, const MultiFab* cBz, + Real t, Real dt) +{ + jx.setVal(0.0); + jy.setVal(0.0); + jz.setVal(0.0); + if (cjx) cjx->setVal(0.0); + if (cjy) cjy->setVal(0.0); + if (cjz) cjz->setVal(0.0); + if (rho) rho->setVal(0.0); + if (crho) crho->setVal(0.0); + for (auto& pc : allcontainers) { + pc->Evolve(lev, Ex, Ey, Ez, Bx, By, Bz, jx, jy, jz, cjx, cjy, cjz, + rho, crho, cEx, cEy, cEz, cBx, cBy, cBz, t, dt); + } +} + +void +MultiParticleContainer::PushX (Real dt) +{ + for (auto& pc : allcontainers) { + pc->PushX(dt); + } +} + +void +MultiParticleContainer::PushP (int lev, Real dt, + const MultiFab& Ex, const MultiFab& Ey, const MultiFab& Ez, + const MultiFab& Bx, const MultiFab& By, const MultiFab& Bz) +{ + for (auto& pc : allcontainers) { + pc->PushP(lev, dt, Ex, Ey, Ez, Bx, By, Bz); + } +} + +std::unique_ptr +MultiParticleContainer::GetChargeDensity (int lev, bool local) +{ + std::unique_ptr rho = allcontainers[0]->GetChargeDensity(lev, true); + for (unsigned i = 1, n = allcontainers.size(); i < n; ++i) { + std::unique_ptr rhoi = allcontainers[i]->GetChargeDensity(lev, true); + MultiFab::Add(*rho, *rhoi, 0, 0, 1, rho->nGrow()); + } + if (!local) { + const Geometry& gm = allcontainers[0]->Geom(lev); + rho->SumBoundary(gm.periodicity()); + } + return rho; +} + +void +MultiParticleContainer::SortParticlesByCell () +{ + for (auto& pc : allcontainers) { + pc->SortParticlesByCell(); + } +} + +void +MultiParticleContainer::Redistribute () +{ + for (auto& pc : allcontainers) { + pc->Redistribute(); + } +} + +void +MultiParticleContainer::RedistributeLocal (const int num_ghost) +{ + for (auto& pc : allcontainers) { + pc->Redistribute(0, 0, 0, num_ghost); + } +} + +Vector +MultiParticleContainer::NumberOfParticlesInGrid(int lev) const +{ + const bool only_valid=true, only_local=true; + Vector r = allcontainers[0]->NumberOfParticlesInGrid(lev,only_valid,only_local); + for (unsigned i = 1, n = allcontainers.size(); i < n; ++i) { + const auto& ri = allcontainers[i]->NumberOfParticlesInGrid(lev,only_valid,only_local); + for (unsigned j=0, m=ri.size(); jIncrement(mf,lev); + } +} + +void +MultiParticleContainer::SetParticleBoxArray (int lev, BoxArray& new_ba) +{ + for (auto& pc : allcontainers) { + pc->SetParticleBoxArray(lev,new_ba); + } +} + +void +MultiParticleContainer::SetParticleDistributionMap (int lev, DistributionMapping& new_dm) +{ + for (auto& pc : allcontainers) { + pc->SetParticleDistributionMap(lev,new_dm); + } +} + +void +MultiParticleContainer::PostRestart () +{ + for (auto& pc : allcontainers) { + pc->PostRestart(); + } + pc_tmp->PostRestart(); +} + +void +MultiParticleContainer +::GetLabFrameData(const std::string& snapshot_name, + const int i_lab, const int direction, + const Real z_old, const Real z_new, + const Real t_boost, const Real t_lab, const Real dt, + Vector& parts) const +{ + + BL_PROFILE("MultiParticleContainer::GetLabFrameData"); + + for (int i = 0; i < nspecies; ++i) + { + WarpXParticleContainer* pc = allcontainers[i].get(); + WarpXParticleContainer::DiagnosticParticles diagnostic_particles; + pc->GetParticleSlice(direction, z_old, z_new, t_boost, t_lab, dt, diagnostic_particles); + + for (int lev = 0; lev <= pc->finestLevel(); ++lev) + { + for (auto it = diagnostic_particles[lev].begin(); it != diagnostic_particles[lev].end(); ++it) + { + parts[i].GetRealData(DiagIdx::w).insert( parts[i].GetRealData(DiagIdx::w ).end(), + it->second.GetRealData(DiagIdx::w ).begin(), + it->second.GetRealData(DiagIdx::w ).end()); + + parts[i].GetRealData(DiagIdx::x).insert( parts[i].GetRealData(DiagIdx::x ).end(), + it->second.GetRealData(DiagIdx::x ).begin(), + it->second.GetRealData(DiagIdx::x ).end()); + + parts[i].GetRealData(DiagIdx::y).insert( parts[i].GetRealData(DiagIdx::y ).end(), + it->second.GetRealData(DiagIdx::y ).begin(), + it->second.GetRealData(DiagIdx::y ).end()); + + parts[i].GetRealData(DiagIdx::z).insert( parts[i].GetRealData(DiagIdx::z ).end(), + it->second.GetRealData(DiagIdx::z ).begin(), + it->second.GetRealData(DiagIdx::z ).end()); + + parts[i].GetRealData(DiagIdx::ux).insert( parts[i].GetRealData(DiagIdx::ux).end(), + it->second.GetRealData(DiagIdx::ux).begin(), + it->second.GetRealData(DiagIdx::ux).end()); + + parts[i].GetRealData(DiagIdx::uy).insert( parts[i].GetRealData(DiagIdx::uy).end(), + it->second.GetRealData(DiagIdx::uy).begin(), + it->second.GetRealData(DiagIdx::uy).end()); + + parts[i].GetRealData(DiagIdx::uz).insert( parts[i].GetRealData(DiagIdx::uz).end(), + it->second.GetRealData(DiagIdx::uz).begin(), + it->second.GetRealData(DiagIdx::uz).end()); + } + } + } +} diff --git a/Source/Particles/ParticleContainer.H b/Source/Particles/ParticleContainer.H deleted file mode 100644 index b21247ff6..000000000 --- a/Source/Particles/ParticleContainer.H +++ /dev/null @@ -1,209 +0,0 @@ - -#ifndef WARPX_ParticleContainer_H_ -#define WARPX_ParticleContainer_H_ - -#include -#include -#include - -#include - -#include -#include -#include -#include - -// -// MultiParticleContainer holds multiple (nspecies or npsecies+1 when -// using laser) WarpXParticleContainer. WarpXParticleContainer is -// derived from amrex::ParticleContainer, and it is -// polymorphic. PhysicalParticleContainer and LaserParticleContainer are -// concrete class dervied from WarpXParticlContainer. -// - -class MultiParticleContainer -{ - -public: - - MultiParticleContainer (amrex::AmrCore* amr_core); - - ~MultiParticleContainer() {} - - WarpXParticleContainer& GetParticleContainer (int ispecies) { - return *allcontainers[ispecies]; - } - - std::array meanParticleVelocity(int ispecies) { - return allcontainers[ispecies]->meanParticleVelocity(); - } - - void AllocData (); - - void InitData (); - -#ifdef WARPX_DO_ELECTROSTATIC - /// - /// Performs the field gather operation using the input field E, for all the species - /// in the MultiParticleContainer. This is the electrostatic version of the field gather. - /// - void FieldGatherES (const amrex::Vector, 3> >& E, - const amrex::Vector > > >& masks); - - /// - /// This evolves all the particles by one PIC time step, including charge deposition, the - /// field solve, and pushing the particles, for all the species in the MultiParticleContainer. - /// This is the electrostatic version. - /// - void EvolveES (const amrex::Vector, 3> >& E, - amrex::Vector >& rho, - amrex::Real t, amrex::Real dt); - - /// - /// This pushes the particle positions by one half time step for all the species in the - /// MultiParticleContainer. It is used to desynchronize the particles after initializaton - /// or when restarting from a checkpoint. This is the electrostatic version. - /// - void PushXES (amrex::Real dt); - - /// - /// This deposits the particle charge onto rho, accumulating the value for all the species - /// in the MultiParticleContainer. rho is assumed to contain node-centered multifabs. - /// This version is hard-coded for CIC deposition. - /// - void DepositCharge(amrex::Vector >& rho, bool local = false); - - /// - /// This returns the total particle charge for all the species in this MultiParticleContainer. - /// This is needed to subtract the offset for periodic boundary conditions. - /// - amrex::Real sumParticleCharge(bool local = false); -#endif // WARPX_DO_ELECTROSTATIC - - /// - /// Performs the field gather operation using the input fields E and B, for all the species - /// in the MultiParticleContainer. This is the electromagnetic version of the field gather. - /// - void FieldGather (int lev, - const amrex::MultiFab& Ex, const amrex::MultiFab& Ey, const amrex::MultiFab& Ez, - const amrex::MultiFab& Bx, const amrex::MultiFab& By, const amrex::MultiFab& Bz); - - /// - /// This evolves all the particles by one PIC time step, including current deposition, the - /// field solve, and pushing the particles, for all the species in the MultiParticleContainer. - /// This is the electromagnetic version. - /// - void Evolve (int lev, - const amrex::MultiFab& Ex, const amrex::MultiFab& Ey, const amrex::MultiFab& Ez, - const amrex::MultiFab& Bx, const amrex::MultiFab& By, const amrex::MultiFab& Bz, - amrex::MultiFab& jx, amrex::MultiFab& jy, amrex::MultiFab& jz, - amrex::MultiFab* cjx, amrex::MultiFab* cjy, amrex::MultiFab* cjz, - amrex::MultiFab* rho, amrex::MultiFab* crho, - const amrex::MultiFab* cEx, const amrex::MultiFab* cEy, const amrex::MultiFab* cEz, - const amrex::MultiFab* cBx, const amrex::MultiFab* cBy, const amrex::MultiFab* cBz, - amrex::Real t, amrex::Real dt); - - /// - /// This pushes the particle positions by one half time step for all the species in the - /// MultiParticleContainer. It is used to desynchronize the particles after initializaton - /// or when restarting from a checkpoint. This is the electromagnetic version. - /// - void PushX (amrex::Real dt); - - /// - /// This pushes the particle momenta by dt for all the species in the - /// MultiParticleContainer. It is used to desynchronize the particles after initializaton - /// or when restarting from a checkpoint. It is also used to synchronize particles at the - /// the end of the run. This is the electromagnetic version. - /// - void PushP (int lev, amrex::Real dt, - const amrex::MultiFab& Ex, const amrex::MultiFab& Ey, const amrex::MultiFab& Ez, - const amrex::MultiFab& Bx, const amrex::MultiFab& By, const amrex::MultiFab& Bz); - - /// - /// This deposits the particle charge onto a node-centered MultiFab and returns a unique ptr - /// to it. The charge density is accumulated over all the particles in the MultiParticleContainer - /// This version uses PICSAR routines to perform the deposition. - /// - std::unique_ptr GetChargeDensity(int lev, bool local = false); - - void Checkpoint (const std::string& dir, - bool is_checkpoint, - const amrex::Vector& varnames = amrex::Vector()) const; - - void Restart (const std::string& dir); - - void PostRestart (); - - void ReadHeader (std::istream& is); - - void WriteHeader (std::ostream& os) const; - - void SortParticlesByCell (); - - void Redistribute (); - - void RedistributeLocal (const int num_ghost); - - amrex::Vector NumberOfParticlesInGrid(int lev) const; - - void Increment (amrex::MultiFab& mf, int lev); - - void SetParticleBoxArray (int lev, amrex::BoxArray& new_ba); - void SetParticleDistributionMap (int lev, amrex::DistributionMapping& new_dm); - - int nSpecies() const {return nspecies;} - - int nSpeciesDepositOnMainGrid () const { - int r = 0; - for (int i : deposit_on_main_grid) { - if (i) ++r; - } - return r; - } - - void GetLabFrameData(const std::string& snapshot_name, - const int i_lab, const int direction, - const amrex::Real z_old, const amrex::Real z_new, - const amrex::Real t_boost, const amrex::Real t_lab, const amrex::Real dt, - amrex::Vector& parts) const; - - // - // Parameters for the Cherenkov corrector in the FDTD solver. - // Both stencils are calculated ar runtime. - // - // Number of coefficients for the stencil of the NCI corrector. - // The stencil is applied in the z direction only. - static constexpr int nstencilz_fdtd_nci_corr=5; - - amrex::Vector > fdtd_nci_stencilz_ex; - amrex::Vector > fdtd_nci_stencilz_by; - - std::vector GetSpeciesNames() const { return species_names; } - - PhysicalParticleContainer& GetPCtmp () { return *pc_tmp; } - -protected: - - // Particle container types - enum struct PCTypes {Physical, RigidInjected}; - - std::vector species_names; - - std::vector deposit_on_main_grid; - - std::vector species_types; - -private: - - // physical particles (+ laser) - amrex::Vector > allcontainers; - // Temporary particle container, used e.g. for particle splitting. - std::unique_ptr pc_tmp; - - void ReadParameters (); - - // runtime parameters - int nspecies = 1; // physical particles only. If WarpX::use_laser, nspecies+1 == allcontainers.size(). -}; -#endif /*WARPX_ParticleContainer_H_*/ diff --git a/Source/Particles/ParticleContainer.cpp b/Source/Particles/ParticleContainer.cpp deleted file mode 100644 index 8a46dd835..000000000 --- a/Source/Particles/ParticleContainer.cpp +++ /dev/null @@ -1,389 +0,0 @@ -#include -#include -#include - -#include -#include -#include - -using namespace amrex; - -constexpr int MultiParticleContainer::nstencilz_fdtd_nci_corr; - -MultiParticleContainer::MultiParticleContainer (AmrCore* amr_core) -{ - ReadParameters(); - - int n = WarpX::use_laser ? nspecies+1 : nspecies; - allcontainers.resize(n); - for (int i = 0; i < nspecies; ++i) { - if (species_types[i] == PCTypes::Physical) { - allcontainers[i].reset(new PhysicalParticleContainer(amr_core, i, species_names[i])); - } - else if (species_types[i] == PCTypes::RigidInjected) { - allcontainers[i].reset(new RigidInjectedParticleContainer(amr_core, i, species_names[i])); - } - allcontainers[i]->deposit_on_main_grid = deposit_on_main_grid[i]; - } - if (WarpX::use_laser) { - allcontainers[n-1].reset(new LaserParticleContainer(amr_core,n-1)); - } - pc_tmp.reset(new PhysicalParticleContainer(amr_core)); -} - -void -MultiParticleContainer::ReadParameters () -{ - static bool initialized = false; - if (!initialized) - { - ParmParse pp("particles"); - - pp.query("nspecies", nspecies); - BL_ASSERT(nspecies >= 0); - - if (nspecies > 0) { - pp.getarr("species_names", species_names); - BL_ASSERT(species_names.size() == nspecies); - - deposit_on_main_grid.resize(nspecies, 0); - std::vector tmp; - pp.queryarr("deposit_on_main_grid", tmp); - for (auto const& name : tmp) { - auto it = std::find(species_names.begin(), species_names.end(), name); - AMREX_ALWAYS_ASSERT_WITH_MESSAGE(it != species_names.end(), "ERROR: species in particles.deposit_on_main_grid must be part of particles.species_names"); - int i = std::distance(species_names.begin(), it); - deposit_on_main_grid[i] = 1; - } - - species_types.resize(nspecies, PCTypes::Physical); - - std::vector rigid_injected_species; - pp.queryarr("rigid_injected_species", rigid_injected_species); - - if (!rigid_injected_species.empty()) { - for (auto const& name : rigid_injected_species) { - auto it = std::find(species_names.begin(), species_names.end(), name); - AMREX_ALWAYS_ASSERT_WITH_MESSAGE(it != species_names.end(), "ERROR: species in particles.rigid_injected_species must be part of particles.species_names"); - int i = std::distance(species_names.begin(), it); - species_types[i] = PCTypes::RigidInjected; - } - } - } - pp.query("use_fdtd_nci_corr", WarpX::use_fdtd_nci_corr); - pp.query("l_lower_order_in_v", WarpX::l_lower_order_in_v); - initialized = true; - } -} - -void -MultiParticleContainer::AllocData () -{ - for (auto& pc : allcontainers) { - pc->AllocData(); - } - pc_tmp->AllocData(); -} - -void -MultiParticleContainer::InitData () -{ - for (auto& pc : allcontainers) { - pc->InitData(); - } - pc_tmp->InitData(); -} - - -#ifdef WARPX_DO_ELECTROSTATIC -void -MultiParticleContainer::FieldGatherES (const Vector, 3> >& E, - const amrex::Vector > > >& masks) -{ - for (auto& pc : allcontainers) { - pc->FieldGatherES(E, masks); - } -} - -void -MultiParticleContainer::EvolveES (const Vector, 3> >& E, - Vector >& rho, - Real t, Real dt) -{ - - int nlevs = rho.size(); - int ng = rho[0]->nGrow(); - - for (unsigned i = 0; i < nlevs; i++) { - rho[i]->setVal(0.0, ng); - } - - for (auto& pc : allcontainers) { - pc->EvolveES(E, rho, t, dt); - } - - for (unsigned i = 0; i < nlevs; i++) { - const Geometry& gm = allcontainers[0]->Geom(i); - rho[i]->SumBoundary(gm.periodicity()); - } -} - -void -MultiParticleContainer::Evolve (int lev, - const MultiFab& Ex, const MultiFab& Ey, const MultiFab& Ez, - const MultiFab& Bx, const MultiFab& By, const MultiFab& Bz, - MultiFab& jx, MultiFab& jy, MultiFab& jz, - MultiFab* cjx, MultiFab* cjy, MultiFab* cjz, - MultiFab* rho, - const MultiFab* cEx, const MultiFab* cEy, const MultiFab* cEz, - const MultiFab* cBx, const MultiFab* cBy, const MultiFab* cBz, - Real t, Real dt) -{ - jx.setVal(0.0); - jy.setVal(0.0); - jz.setVal(0.0); - if (cjx) cjx->setVal(0.0); - if (cjy) cjy->setVal(0.0); - if (cjz) cjz->setVal(0.0); - if (rho) rho->setVal(0.0); - for (auto& pc : allcontainers) { - pc->Evolve(lev, Ex, Ey, Ez, Bx, By, Bz, jx, jy, jz, cjx, cjy, cjz, - rho, cEx, cEy, cEz, cBx, cBy, cBz, t, dt); - } -} - -void -MultiParticleContainer::PushXES (Real dt) -{ - for (auto& pc : allcontainers) { - pc->PushXES(dt); - } -} - -void -MultiParticleContainer:: -DepositCharge (Vector >& rho, bool local) -{ - int nlevs = rho.size(); - int ng = rho[0]->nGrow(); - - for (unsigned i = 0; i < nlevs; i++) { - rho[i]->setVal(0.0, ng); - } - - for (unsigned i = 0, n = allcontainers.size(); i < n; ++i) { - allcontainers[i]->DepositCharge(rho, true); - } - - if (!local) { - for (unsigned i = 0; i < nlevs; i++) { - const Geometry& gm = allcontainers[0]->Geom(i); - rho[i]->SumBoundary(gm.periodicity()); - } - } -} - -amrex::Real -MultiParticleContainer::sumParticleCharge (bool local) -{ - amrex::Real total_charge = allcontainers[0]->sumParticleCharge(local); - for (unsigned i = 1, n = allcontainers.size(); i < n; ++i) { - total_charge += allcontainers[i]->sumParticleCharge(local); - } - return total_charge; -} - -#endif // WARPX_DO_ELECTROSTATIC - -void -MultiParticleContainer::FieldGather (int lev, - const MultiFab& Ex, const MultiFab& Ey, const MultiFab& Ez, - const MultiFab& Bx, const MultiFab& By, const MultiFab& Bz) -{ - for (auto& pc : allcontainers) { - pc->FieldGather(lev, Ex, Ey, Ez, Bx, By, Bz); - } -} - -void -MultiParticleContainer::Evolve (int lev, - const MultiFab& Ex, const MultiFab& Ey, const MultiFab& Ez, - const MultiFab& Bx, const MultiFab& By, const MultiFab& Bz, - MultiFab& jx, MultiFab& jy, MultiFab& jz, - MultiFab* cjx, MultiFab* cjy, MultiFab* cjz, - MultiFab* rho, MultiFab* crho, - const MultiFab* cEx, const MultiFab* cEy, const MultiFab* cEz, - const MultiFab* cBx, const MultiFab* cBy, const MultiFab* cBz, - Real t, Real dt) -{ - jx.setVal(0.0); - jy.setVal(0.0); - jz.setVal(0.0); - if (cjx) cjx->setVal(0.0); - if (cjy) cjy->setVal(0.0); - if (cjz) cjz->setVal(0.0); - if (rho) rho->setVal(0.0); - if (crho) crho->setVal(0.0); - for (auto& pc : allcontainers) { - pc->Evolve(lev, Ex, Ey, Ez, Bx, By, Bz, jx, jy, jz, cjx, cjy, cjz, - rho, crho, cEx, cEy, cEz, cBx, cBy, cBz, t, dt); - } -} - -void -MultiParticleContainer::PushX (Real dt) -{ - for (auto& pc : allcontainers) { - pc->PushX(dt); - } -} - -void -MultiParticleContainer::PushP (int lev, Real dt, - const MultiFab& Ex, const MultiFab& Ey, const MultiFab& Ez, - const MultiFab& Bx, const MultiFab& By, const MultiFab& Bz) -{ - for (auto& pc : allcontainers) { - pc->PushP(lev, dt, Ex, Ey, Ez, Bx, By, Bz); - } -} - -std::unique_ptr -MultiParticleContainer::GetChargeDensity (int lev, bool local) -{ - std::unique_ptr rho = allcontainers[0]->GetChargeDensity(lev, true); - for (unsigned i = 1, n = allcontainers.size(); i < n; ++i) { - std::unique_ptr rhoi = allcontainers[i]->GetChargeDensity(lev, true); - MultiFab::Add(*rho, *rhoi, 0, 0, 1, rho->nGrow()); - } - if (!local) { - const Geometry& gm = allcontainers[0]->Geom(lev); - rho->SumBoundary(gm.periodicity()); - } - return rho; -} - -void -MultiParticleContainer::SortParticlesByCell () -{ - for (auto& pc : allcontainers) { - pc->SortParticlesByCell(); - } -} - -void -MultiParticleContainer::Redistribute () -{ - for (auto& pc : allcontainers) { - pc->Redistribute(); - } -} - -void -MultiParticleContainer::RedistributeLocal (const int num_ghost) -{ - for (auto& pc : allcontainers) { - pc->Redistribute(0, 0, 0, num_ghost); - } -} - -Vector -MultiParticleContainer::NumberOfParticlesInGrid(int lev) const -{ - const bool only_valid=true, only_local=true; - Vector r = allcontainers[0]->NumberOfParticlesInGrid(lev,only_valid,only_local); - for (unsigned i = 1, n = allcontainers.size(); i < n; ++i) { - const auto& ri = allcontainers[i]->NumberOfParticlesInGrid(lev,only_valid,only_local); - for (unsigned j=0, m=ri.size(); jIncrement(mf,lev); - } -} - -void -MultiParticleContainer::SetParticleBoxArray (int lev, BoxArray& new_ba) -{ - for (auto& pc : allcontainers) { - pc->SetParticleBoxArray(lev,new_ba); - } -} - -void -MultiParticleContainer::SetParticleDistributionMap (int lev, DistributionMapping& new_dm) -{ - for (auto& pc : allcontainers) { - pc->SetParticleDistributionMap(lev,new_dm); - } -} - -void -MultiParticleContainer::PostRestart () -{ - for (auto& pc : allcontainers) { - pc->PostRestart(); - } - pc_tmp->PostRestart(); -} - -void -MultiParticleContainer -::GetLabFrameData(const std::string& snapshot_name, - const int i_lab, const int direction, - const Real z_old, const Real z_new, - const Real t_boost, const Real t_lab, const Real dt, - Vector& parts) const -{ - - BL_PROFILE("MultiParticleContainer::GetLabFrameData"); - - for (int i = 0; i < nspecies; ++i) - { - WarpXParticleContainer* pc = allcontainers[i].get(); - WarpXParticleContainer::DiagnosticParticles diagnostic_particles; - pc->GetParticleSlice(direction, z_old, z_new, t_boost, t_lab, dt, diagnostic_particles); - - for (int lev = 0; lev <= pc->finestLevel(); ++lev) - { - for (auto it = diagnostic_particles[lev].begin(); it != diagnostic_particles[lev].end(); ++it) - { - parts[i].GetRealData(DiagIdx::w).insert( parts[i].GetRealData(DiagIdx::w ).end(), - it->second.GetRealData(DiagIdx::w ).begin(), - it->second.GetRealData(DiagIdx::w ).end()); - - parts[i].GetRealData(DiagIdx::x).insert( parts[i].GetRealData(DiagIdx::x ).end(), - it->second.GetRealData(DiagIdx::x ).begin(), - it->second.GetRealData(DiagIdx::x ).end()); - - parts[i].GetRealData(DiagIdx::y).insert( parts[i].GetRealData(DiagIdx::y ).end(), - it->second.GetRealData(DiagIdx::y ).begin(), - it->second.GetRealData(DiagIdx::y ).end()); - - parts[i].GetRealData(DiagIdx::z).insert( parts[i].GetRealData(DiagIdx::z ).end(), - it->second.GetRealData(DiagIdx::z ).begin(), - it->second.GetRealData(DiagIdx::z ).end()); - - parts[i].GetRealData(DiagIdx::ux).insert( parts[i].GetRealData(DiagIdx::ux).end(), - it->second.GetRealData(DiagIdx::ux).begin(), - it->second.GetRealData(DiagIdx::ux).end()); - - parts[i].GetRealData(DiagIdx::uy).insert( parts[i].GetRealData(DiagIdx::uy).end(), - it->second.GetRealData(DiagIdx::uy).begin(), - it->second.GetRealData(DiagIdx::uy).end()); - - parts[i].GetRealData(DiagIdx::uz).insert( parts[i].GetRealData(DiagIdx::uz).end(), - it->second.GetRealData(DiagIdx::uz).begin(), - it->second.GetRealData(DiagIdx::uz).end()); - } - } - } -} diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 9cfd30862..afc1412d5 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1,7 +1,7 @@ #include #include -#include +#include #include #include #include diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index ca7d18d19..ad80f7c4f 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -1,7 +1,7 @@ #include -#include +#include #include #include #include diff --git a/Source/Python/Make.package b/Source/Python/Make.package index c79049f02..71bd4ebe8 100644 --- a/Source/Python/Make.package +++ b/Source/Python/Make.package @@ -2,8 +2,9 @@ ifeq ($(USE_PYTHON_MAIN),TRUE) CEXE_sources += WarpXWrappers.cpp CEXE_headers += WarpXWrappers.h endif - +CEXE_sources += WarpXWrappers.cpp CEXE_sources += WarpX_py.cpp +CEXE_headers += WarpXWrappers.h CEXE_headers += WarpX_py.H INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Python diff --git a/Source/Python/WarpXWrappers.cpp b/Source/Python/WarpXWrappers.cpp new file mode 100644 index 000000000..6ee06a7d5 --- /dev/null +++ b/Source/Python/WarpXWrappers.cpp @@ -0,0 +1,446 @@ + +#include +#include + +#include +#include +#include +#include + +namespace +{ + double** getMultiFabPointers(const amrex::MultiFab& mf, int *num_boxes, int *ngrow, int **shapes) + { + *ngrow = mf.nGrow(); + *num_boxes = mf.local_size(); + *shapes = (int*) malloc(AMREX_SPACEDIM * (*num_boxes) * sizeof(int)); + double** data = (double**) malloc((*num_boxes) * sizeof(double*)); + + int i = 0; +#ifdef _OPENMP +#pragma omp parallel +#endif + for ( amrex::MFIter mfi(mf, false); mfi.isValid(); ++mfi, ++i ) { + data[i] = (double*) mf[mfi].dataPtr(); + for (int j = 0; j < AMREX_SPACEDIM; ++j) { + (*shapes)[AMREX_SPACEDIM*i+j] = mf[mfi].box().length(j); + } + } + return data; + } + int* getMultiFabLoVects(const amrex::MultiFab& mf, int *num_boxes, int *ngrow) + { + *ngrow = mf.nGrow(); + *num_boxes = mf.local_size(); + int *loVects = (int*) malloc((*num_boxes)*AMREX_SPACEDIM * sizeof(int)); + + int i = 0; + for ( amrex::MFIter mfi(mf, false); mfi.isValid(); ++mfi, ++i ) { + const int* loVect = mf[mfi].loVect(); + for (int j = 0; j < AMREX_SPACEDIM; ++j) { + loVects[AMREX_SPACEDIM*i+j] = loVect[j]; + } + } + return loVects; + } +} + +extern "C" +{ + + int warpx_nSpecies() + { + auto & mypc = WarpX::GetInstance().GetPartContainer(); + return mypc.nSpecies(); + } + + bool warpx_use_fdtd_nci_corr() + { + return WarpX::use_fdtd_nci_corr; + } + + int warpx_l_lower_order_in_v() + { + return WarpX::l_lower_order_in_v; + } + + int warpx_nComps() + { + return PIdx::nattribs; + } + + int warpx_SpaceDim() + { + return AMREX_SPACEDIM; + } + + void amrex_init (int argc, char* argv[]) + { + amrex::Initialize(argc,argv); + } + +#ifdef BL_USE_MPI + void amrex_init_with_inited_mpi (int argc, char* argv[], MPI_Comm mpicomm) + { + amrex::Initialize(argc,argv,true,mpicomm); + } +#endif + + void amrex_finalize (int finalize_mpi) + { + amrex::Finalize(finalize_mpi); + } + + void warpx_init () + { + WarpX& warpx = WarpX::GetInstance(); + warpx.InitData(); + if (warpx_py_afterinit) warpx_py_afterinit(); + if (warpx_py_particleloader) warpx_py_particleloader(); + } + + void warpx_finalize () + { + WarpX::ResetInstance(); + } + + void warpx_set_callback_py_afterinit (WARPX_CALLBACK_PY_FUNC_0 callback) + { + warpx_py_afterinit = callback; + } + void warpx_set_callback_py_beforeEsolve (WARPX_CALLBACK_PY_FUNC_0 callback) + { + warpx_py_beforeEsolve = callback; + } + void warpx_set_callback_py_afterEsolve (WARPX_CALLBACK_PY_FUNC_0 callback) + { + warpx_py_afterEsolve = callback; + } + void warpx_set_callback_py_beforedeposition (WARPX_CALLBACK_PY_FUNC_0 callback) + { + warpx_py_beforedeposition = callback; + } + void warpx_set_callback_py_afterdeposition (WARPX_CALLBACK_PY_FUNC_0 callback) + { + warpx_py_afterdeposition = callback; + } + void warpx_set_callback_py_particlescraper (WARPX_CALLBACK_PY_FUNC_0 callback) + { + warpx_py_particlescraper = callback; + } + void warpx_set_callback_py_particleloader (WARPX_CALLBACK_PY_FUNC_0 callback) + { + warpx_py_particleloader = callback; + } + void warpx_set_callback_py_beforestep (WARPX_CALLBACK_PY_FUNC_0 callback) + { + warpx_py_beforestep = callback; + } + void warpx_set_callback_py_afterstep (WARPX_CALLBACK_PY_FUNC_0 callback) + { + warpx_py_afterstep = callback; + } + void warpx_set_callback_py_afterrestart (WARPX_CALLBACK_PY_FUNC_0 callback) + { + warpx_py_afterrestart = callback; + } + void warpx_set_callback_py_particleinjection (WARPX_CALLBACK_PY_FUNC_0 callback) + { + warpx_py_particleinjection = callback; + } + void warpx_set_callback_py_appliedfields (WARPX_CALLBACK_PY_FUNC_0 callback) + { + warpx_py_appliedfields = callback; + } + + void warpx_evolve (int numsteps) + { + WarpX& warpx = WarpX::GetInstance(); + warpx.Evolve(numsteps); + } + + void warpx_addNParticles(int speciesnumber, int lenx, + double* x, double* y, double* z, + double* vx, double* vy, double* vz, + int nattr, double* attr, int uniqueparticles) + { + auto & mypc = WarpX::GetInstance().GetPartContainer(); + auto & myspc = mypc.GetParticleContainer(speciesnumber); + const int lev = 0; + myspc.AddNParticles(lev, lenx, x, y, z, vx, vy, vz, nattr, attr, uniqueparticles); + } + + double warpx_getProbLo(int dir) + { + WarpX& warpx = WarpX::GetInstance(); + const amrex::Geometry& geom = warpx.Geom(0); + return geom.ProbLo(dir); + } + + double warpx_getProbHi(int dir) + { + WarpX& warpx = WarpX::GetInstance(); + const amrex::Geometry& geom = warpx.Geom(0); + return geom.ProbHi(dir); + } + + long warpx_getNumParticles(int speciesnumber) { + auto & mypc = WarpX::GetInstance().GetPartContainer(); + auto & myspc = mypc.GetParticleContainer(speciesnumber); + return myspc.TotalNumberOfParticles(); + } + + double** warpx_getEfield(int lev, int direction, + int *return_size, int *ngrow, int **shapes) { + auto & mf = WarpX::GetInstance().getEfield(lev, direction); + return getMultiFabPointers(mf, return_size, ngrow, shapes); + } + + int* warpx_getEfieldLoVects(int lev, int direction, + int *return_size, int *ngrow) { + auto & mf = WarpX::GetInstance().getEfield(lev, direction); + return getMultiFabLoVects(mf, return_size, ngrow); + } + + double** warpx_getEfieldCP(int lev, int direction, + int *return_size, int *ngrow, int **shapes) { + auto & mf = WarpX::GetInstance().getEfield_cp(lev, direction); + return getMultiFabPointers(mf, return_size, ngrow, shapes); + } + + int* warpx_getEfieldCPLoVects(int lev, int direction, + int *return_size, int *ngrow) { + auto & mf = WarpX::GetInstance().getEfield_cp(lev, direction); + return getMultiFabLoVects(mf, return_size, ngrow); + } + + double** warpx_getEfieldFP(int lev, int direction, + int *return_size, int *ngrow, int **shapes) { + auto & mf = WarpX::GetInstance().getEfield_fp(lev, direction); + return getMultiFabPointers(mf, return_size, ngrow, shapes); + } + + int* warpx_getEfieldFPLoVects(int lev, int direction, + int *return_size, int *ngrow) { + auto & mf = WarpX::GetInstance().getEfield_fp(lev, direction); + return getMultiFabLoVects(mf, return_size, ngrow); + } + + double** warpx_getBfield(int lev, int direction, + int *return_size, int *ngrow, int **shapes) { + auto & mf = WarpX::GetInstance().getBfield(lev, direction); + return getMultiFabPointers(mf, return_size, ngrow, shapes); + } + + int* warpx_getBfieldLoVects(int lev, int direction, + int *return_size, int *ngrow) { + auto & mf = WarpX::GetInstance().getBfield(lev, direction); + return getMultiFabLoVects(mf, return_size, ngrow); + } + + double** warpx_getBfieldCP(int lev, int direction, + int *return_size, int *ngrow, int **shapes) { + auto & mf = WarpX::GetInstance().getBfield_cp(lev, direction); + return getMultiFabPointers(mf, return_size, ngrow, shapes); + } + + int* warpx_getBfieldCPLoVects(int lev, int direction, + int *return_size, int *ngrow) { + auto & mf = WarpX::GetInstance().getBfield_cp(lev, direction); + return getMultiFabLoVects(mf, return_size, ngrow); + } + + double** warpx_getBfieldFP(int lev, int direction, + int *return_size, int *ngrow, int **shapes) { + auto & mf = WarpX::GetInstance().getBfield_fp(lev, direction); + return getMultiFabPointers(mf, return_size, ngrow, shapes); + } + + int* warpx_getBfieldFPLoVects(int lev, int direction, + int *return_size, int *ngrow) { + auto & mf = WarpX::GetInstance().getBfield_fp(lev, direction); + return getMultiFabLoVects(mf, return_size, ngrow); + } + + double** warpx_getCurrentDensity(int lev, int direction, + int *return_size, int *ngrow, int **shapes) { + auto & mf = WarpX::GetInstance().getcurrent(lev, direction); + return getMultiFabPointers(mf, return_size, ngrow, shapes); + } + + int* warpx_getCurrentDensityLoVects(int lev, int direction, + int *return_size, int *ngrow) { + auto & mf = WarpX::GetInstance().getcurrent(lev, direction); + return getMultiFabLoVects(mf, return_size, ngrow); + } + + double** warpx_getCurrentDensityCP(int lev, int direction, + int *return_size, int *ngrow, int **shapes) { + auto & mf = WarpX::GetInstance().getcurrent_cp(lev, direction); + return getMultiFabPointers(mf, return_size, ngrow, shapes); + } + + int* warpx_getCurrentDensityCPLoVects(int lev, int direction, + int *return_size, int *ngrow) { + auto & mf = WarpX::GetInstance().getcurrent_cp(lev, direction); + return getMultiFabLoVects(mf, return_size, ngrow); + } + + double** warpx_getCurrentDensityFP(int lev, int direction, + int *return_size, int *ngrow, int **shapes) { + auto & mf = WarpX::GetInstance().getcurrent_fp(lev, direction); + return getMultiFabPointers(mf, return_size, ngrow, shapes); + } + + int* warpx_getCurrentDensityFPLoVects(int lev, int direction, + int *return_size, int *ngrow) { + auto & mf = WarpX::GetInstance().getcurrent_fp(lev, direction); + return getMultiFabLoVects(mf, return_size, ngrow); + } + + double** warpx_getParticleStructs(int speciesnumber, + int* num_tiles, int** particles_per_tile) { + auto & mypc = WarpX::GetInstance().GetPartContainer(); + auto & myspc = mypc.GetParticleContainer(speciesnumber); + + const int level = 0; + + int i = 0; + for (WarpXParIter pti(myspc, level); pti.isValid(); ++pti, ++i) {} + + // *num_tiles = myspc.numLocalTilesAtLevel(level); + *num_tiles = i; + *particles_per_tile = (int*) malloc(*num_tiles*sizeof(int)); + + double** data = (double**) malloc(*num_tiles*sizeof(typename WarpXParticleContainer::ParticleType*)); + i = 0; + for (WarpXParIter pti(myspc, level); pti.isValid(); ++pti, ++i) { + auto& aos = pti.GetArrayOfStructs(); + data[i] = (double*) aos.data(); + (*particles_per_tile)[i] = pti.numParticles(); + } + return data; + } + + double** warpx_getParticleArrays(int speciesnumber, int comp, + int* num_tiles, int** particles_per_tile) { + auto & mypc = WarpX::GetInstance().GetPartContainer(); + auto & myspc = mypc.GetParticleContainer(speciesnumber); + + const int level = 0; + + int i = 0; + for (WarpXParIter pti(myspc, level); pti.isValid(); ++pti, ++i) {} + + // *num_tiles = myspc.numLocalTilesAtLevel(level); + *num_tiles = i; + *particles_per_tile = (int*) malloc(*num_tiles*sizeof(int)); + + double** data = (double**) malloc(*num_tiles*sizeof(double*)); + i = 0; + for (WarpXParIter pti(myspc, level); pti.isValid(); ++pti, ++i) { + auto& soa = pti.GetStructOfArrays(); + data[i] = (double*) soa.GetRealData(comp).dataPtr(); + (*particles_per_tile)[i] = pti.numParticles(); + } + return data; + } + + void warpx_ComputeDt () { + WarpX& warpx = WarpX::GetInstance(); + warpx.ComputeDt (); + } + void warpx_MoveWindow () { + WarpX& warpx = WarpX::GetInstance(); + warpx.MoveWindow (true); + } + + void warpx_EvolveE (double dt) { + WarpX& warpx = WarpX::GetInstance(); + warpx.EvolveE (dt); + } + void warpx_EvolveB (double dt) { + WarpX& warpx = WarpX::GetInstance(); + warpx.EvolveB (dt); + } + void warpx_FillBoundaryE () { + WarpX& warpx = WarpX::GetInstance(); + warpx.FillBoundaryE (); + } + void warpx_FillBoundaryB () { + WarpX& warpx = WarpX::GetInstance(); + warpx.FillBoundaryB (); + } + void warpx_SyncCurrent () { + WarpX& warpx = WarpX::GetInstance(); + warpx.SyncCurrent (); + } + void warpx_UpdateAuxilaryData () { + WarpX& warpx = WarpX::GetInstance(); + warpx.UpdateAuxilaryData (); + } + void warpx_PushParticlesandDepose (double cur_time) { + WarpX& warpx = WarpX::GetInstance(); + warpx.PushParticlesandDepose (cur_time); + } + + int warpx_getistep (int lev) { + WarpX& warpx = WarpX::GetInstance(); + return warpx.getistep (lev); + } + void warpx_setistep (int lev, int ii) { + WarpX& warpx = WarpX::GetInstance(); + warpx.setistep (lev, ii); + } + double warpx_gett_new (int lev) { + WarpX& warpx = WarpX::GetInstance(); + return warpx.gett_new (lev); + } + void warpx_sett_new (int lev, double time) { + WarpX& warpx = WarpX::GetInstance(); + warpx.sett_new (lev, time); + } + double warpx_getdt (int lev) { + WarpX& warpx = WarpX::GetInstance(); + return warpx.getdt (lev); + } + + int warpx_maxStep () { + WarpX& warpx = WarpX::GetInstance(); + return warpx.maxStep (); + } + double warpx_stopTime () { + WarpX& warpx = WarpX::GetInstance(); + return warpx.stopTime (); + } + + int warpx_checkInt () { + WarpX& warpx = WarpX::GetInstance(); + return warpx.checkInt (); + } + int warpx_plotInt () { + WarpX& warpx = WarpX::GetInstance(); + return warpx.plotInt (); + } + + void warpx_WriteCheckPointFile () { + WarpX& warpx = WarpX::GetInstance(); + warpx.WriteCheckPointFile (); + } + void warpx_WritePlotFile () { + WarpX& warpx = WarpX::GetInstance(); + warpx.WritePlotFile (); + } + + int warpx_finestLevel () { + WarpX& warpx = WarpX::GetInstance(); + return warpx.finestLevel (); + } + + void mypc_Redistribute () { + auto & mypc = WarpX::GetInstance().GetPartContainer(); + mypc.Redistribute(); + } + +} + diff --git a/Source/Python/WarpXWrappers.h b/Source/Python/WarpXWrappers.h new file mode 100644 index 000000000..07d6f80f7 --- /dev/null +++ b/Source/Python/WarpXWrappers.h @@ -0,0 +1,120 @@ +#ifndef WARPX_WRAPPERS_H_ +#define WARPX_WRAPPERS_H_ + +#ifdef BL_USE_MPI +#include +#endif + +#ifdef __cplusplus +extern "C" { +#endif + + int warpx_nSpecies(); + + bool warpx_use_fdtd_nci_corr(); + + int warpx_l_lower_order_in_v(); + + int warpx_nComps(); + + int warpx_SpaceDim(); + + void amrex_init (int argc, char* argv[]); + +#ifdef BL_USE_MPI + void amrex_init_with_inited_mpi (int argc, char* argv[], MPI_Comm mpicomm); +#endif + + void amrex_finalize (int finalize_mpi); + + void warpx_init (); + + void warpx_finalize (); + + typedef void(*WARPX_CALLBACK_PY_FUNC_0)(); + + void warpx_set_callback_py_afterinit (WARPX_CALLBACK_PY_FUNC_0); + void warpx_set_callback_py_beforeEsolve (WARPX_CALLBACK_PY_FUNC_0); + void warpx_set_callback_py_afterEsolve (WARPX_CALLBACK_PY_FUNC_0); + void warpx_set_callback_py_beforedeposition (WARPX_CALLBACK_PY_FUNC_0); + void warpx_set_callback_py_afterdeposition (WARPX_CALLBACK_PY_FUNC_0); + void warpx_set_callback_py_particlescraper (WARPX_CALLBACK_PY_FUNC_0); + void warpx_set_callback_py_particleloader (WARPX_CALLBACK_PY_FUNC_0); + void warpx_set_callback_py_beforestep (WARPX_CALLBACK_PY_FUNC_0); + void warpx_set_callback_py_afterstep (WARPX_CALLBACK_PY_FUNC_0); + void warpx_set_callback_py_afterrestart (WARPX_CALLBACK_PY_FUNC_0); + void warpx_set_callback_py_particleinjection (WARPX_CALLBACK_PY_FUNC_0); + void warpx_set_callback_py_appliedfields (WARPX_CALLBACK_PY_FUNC_0); + + void warpx_evolve (int numsteps); // -1 means the inputs parameter will be used. + + void warpx_addNParticles(int speciesnumber, int lenx, + double* x, double* y, double* z, + double* vx, double* vy, double* vz, + int nattr, double* attr, int uniqueparticles); + + double warpx_getProbLo(int dir); + + double warpx_getProbHi(int dir); + + long warpx_getNumParticles(int speciesnumber); + + double** warpx_getEfield(int lev, int direction, + int *return_size, int* ngrow, int **shapes); + + int* warpx_getEfieldLoVects(int lev, int direction, + int *return_size, int* ngrow); + + double** warpx_getBfield(int lev, int direction, + int *return_size, int* ngrow, int **shapes); + + int* warpx_getBfieldLoVects(int lev, int direction, + int *return_size, int* ngrow); + + double** warpx_getCurrentDensity(int lev, int direction, + int *return_size, int* ngrow, int **shapes); + + int* warpx_getCurrentDensityLoVects(int lev, int direction, + int *return_size, int* ngrow); + + double** warpx_getParticleStructs(int speciesnumber, + int* num_tiles, int** particles_per_tile); + + double** warpx_getParticleArrays(int speciesnumber, int comp, + int* num_tiles, int** particles_per_tile); + + void warpx_ComputeDt (); + void warpx_MoveWindow (); + + void warpx_EvolveE (double dt); + void warpx_EvolveB (double dt); + void warpx_FillBoundaryE (); + void warpx_FillBoundaryB (); + void warpx_SyncCurrent (); + void warpx_UpdateAuxilaryData (); + void warpx_PushParticlesandDepose (double cur_time); + + int warpx_getistep (int lev); + void warpx_setistep (int lev, int ii); + double warpx_gett_new (int lev); + void warpx_sett_new (int lev, double time); + double warpx_getdt (int lev); + + int warpx_maxStep (); + double warpx_stopTime (); + + int warpx_checkInt (); + int warpx_plotInt (); + + void warpx_WriteCheckPointFile (); + void warpx_WritePlotFile (); + + int warpx_finestLevel (); + + void mypc_Redistribute (); + +#ifdef __cplusplus +} +#endif + +#endif diff --git a/Source/Utils/Make.package b/Source/Utils/Make.package index 5e4233a29..8cd124ea3 100644 --- a/Source/Utils/Make.package +++ b/Source/Utils/Make.package @@ -1,12 +1,10 @@ F90EXE_sources += utils_ES.F90 CEXE_sources += WarpXMovingWindow.cpp -CEXE_sources += WarpXWrappers.cpp CEXE_sources += WarpXConst.cpp CEXE_sources += WarpXTagging.cpp CEXE_sources += WarpXUtil.cpp CEXE_headers += WarpXConst.H CEXE_headers += WarpXUtil.H -CEXE_headers += WarpXWrappers.h INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Utils VPATH_LOCATIONS += $(WARPX_HOME)/Source/Utils diff --git a/Source/Utils/WarpXConst.cpp b/Source/Utils/WarpXConst.cpp index 4ca729872..bd3ebc3ef 100644 --- a/Source/Utils/WarpXConst.cpp +++ b/Source/Utils/WarpXConst.cpp @@ -8,7 +8,7 @@ #include #include #include -#include +#include std::string UserConstants::replaceStringValue(std::string math_expr){ std::string pattern, value_str; diff --git a/Source/Utils/WarpXWrappers.cpp b/Source/Utils/WarpXWrappers.cpp deleted file mode 100644 index 6ee06a7d5..000000000 --- a/Source/Utils/WarpXWrappers.cpp +++ /dev/null @@ -1,446 +0,0 @@ - -#include -#include - -#include -#include -#include -#include - -namespace -{ - double** getMultiFabPointers(const amrex::MultiFab& mf, int *num_boxes, int *ngrow, int **shapes) - { - *ngrow = mf.nGrow(); - *num_boxes = mf.local_size(); - *shapes = (int*) malloc(AMREX_SPACEDIM * (*num_boxes) * sizeof(int)); - double** data = (double**) malloc((*num_boxes) * sizeof(double*)); - - int i = 0; -#ifdef _OPENMP -#pragma omp parallel -#endif - for ( amrex::MFIter mfi(mf, false); mfi.isValid(); ++mfi, ++i ) { - data[i] = (double*) mf[mfi].dataPtr(); - for (int j = 0; j < AMREX_SPACEDIM; ++j) { - (*shapes)[AMREX_SPACEDIM*i+j] = mf[mfi].box().length(j); - } - } - return data; - } - int* getMultiFabLoVects(const amrex::MultiFab& mf, int *num_boxes, int *ngrow) - { - *ngrow = mf.nGrow(); - *num_boxes = mf.local_size(); - int *loVects = (int*) malloc((*num_boxes)*AMREX_SPACEDIM * sizeof(int)); - - int i = 0; - for ( amrex::MFIter mfi(mf, false); mfi.isValid(); ++mfi, ++i ) { - const int* loVect = mf[mfi].loVect(); - for (int j = 0; j < AMREX_SPACEDIM; ++j) { - loVects[AMREX_SPACEDIM*i+j] = loVect[j]; - } - } - return loVects; - } -} - -extern "C" -{ - - int warpx_nSpecies() - { - auto & mypc = WarpX::GetInstance().GetPartContainer(); - return mypc.nSpecies(); - } - - bool warpx_use_fdtd_nci_corr() - { - return WarpX::use_fdtd_nci_corr; - } - - int warpx_l_lower_order_in_v() - { - return WarpX::l_lower_order_in_v; - } - - int warpx_nComps() - { - return PIdx::nattribs; - } - - int warpx_SpaceDim() - { - return AMREX_SPACEDIM; - } - - void amrex_init (int argc, char* argv[]) - { - amrex::Initialize(argc,argv); - } - -#ifdef BL_USE_MPI - void amrex_init_with_inited_mpi (int argc, char* argv[], MPI_Comm mpicomm) - { - amrex::Initialize(argc,argv,true,mpicomm); - } -#endif - - void amrex_finalize (int finalize_mpi) - { - amrex::Finalize(finalize_mpi); - } - - void warpx_init () - { - WarpX& warpx = WarpX::GetInstance(); - warpx.InitData(); - if (warpx_py_afterinit) warpx_py_afterinit(); - if (warpx_py_particleloader) warpx_py_particleloader(); - } - - void warpx_finalize () - { - WarpX::ResetInstance(); - } - - void warpx_set_callback_py_afterinit (WARPX_CALLBACK_PY_FUNC_0 callback) - { - warpx_py_afterinit = callback; - } - void warpx_set_callback_py_beforeEsolve (WARPX_CALLBACK_PY_FUNC_0 callback) - { - warpx_py_beforeEsolve = callback; - } - void warpx_set_callback_py_afterEsolve (WARPX_CALLBACK_PY_FUNC_0 callback) - { - warpx_py_afterEsolve = callback; - } - void warpx_set_callback_py_beforedeposition (WARPX_CALLBACK_PY_FUNC_0 callback) - { - warpx_py_beforedeposition = callback; - } - void warpx_set_callback_py_afterdeposition (WARPX_CALLBACK_PY_FUNC_0 callback) - { - warpx_py_afterdeposition = callback; - } - void warpx_set_callback_py_particlescraper (WARPX_CALLBACK_PY_FUNC_0 callback) - { - warpx_py_particlescraper = callback; - } - void warpx_set_callback_py_particleloader (WARPX_CALLBACK_PY_FUNC_0 callback) - { - warpx_py_particleloader = callback; - } - void warpx_set_callback_py_beforestep (WARPX_CALLBACK_PY_FUNC_0 callback) - { - warpx_py_beforestep = callback; - } - void warpx_set_callback_py_afterstep (WARPX_CALLBACK_PY_FUNC_0 callback) - { - warpx_py_afterstep = callback; - } - void warpx_set_callback_py_afterrestart (WARPX_CALLBACK_PY_FUNC_0 callback) - { - warpx_py_afterrestart = callback; - } - void warpx_set_callback_py_particleinjection (WARPX_CALLBACK_PY_FUNC_0 callback) - { - warpx_py_particleinjection = callback; - } - void warpx_set_callback_py_appliedfields (WARPX_CALLBACK_PY_FUNC_0 callback) - { - warpx_py_appliedfields = callback; - } - - void warpx_evolve (int numsteps) - { - WarpX& warpx = WarpX::GetInstance(); - warpx.Evolve(numsteps); - } - - void warpx_addNParticles(int speciesnumber, int lenx, - double* x, double* y, double* z, - double* vx, double* vy, double* vz, - int nattr, double* attr, int uniqueparticles) - { - auto & mypc = WarpX::GetInstance().GetPartContainer(); - auto & myspc = mypc.GetParticleContainer(speciesnumber); - const int lev = 0; - myspc.AddNParticles(lev, lenx, x, y, z, vx, vy, vz, nattr, attr, uniqueparticles); - } - - double warpx_getProbLo(int dir) - { - WarpX& warpx = WarpX::GetInstance(); - const amrex::Geometry& geom = warpx.Geom(0); - return geom.ProbLo(dir); - } - - double warpx_getProbHi(int dir) - { - WarpX& warpx = WarpX::GetInstance(); - const amrex::Geometry& geom = warpx.Geom(0); - return geom.ProbHi(dir); - } - - long warpx_getNumParticles(int speciesnumber) { - auto & mypc = WarpX::GetInstance().GetPartContainer(); - auto & myspc = mypc.GetParticleContainer(speciesnumber); - return myspc.TotalNumberOfParticles(); - } - - double** warpx_getEfield(int lev, int direction, - int *return_size, int *ngrow, int **shapes) { - auto & mf = WarpX::GetInstance().getEfield(lev, direction); - return getMultiFabPointers(mf, return_size, ngrow, shapes); - } - - int* warpx_getEfieldLoVects(int lev, int direction, - int *return_size, int *ngrow) { - auto & mf = WarpX::GetInstance().getEfield(lev, direction); - return getMultiFabLoVects(mf, return_size, ngrow); - } - - double** warpx_getEfieldCP(int lev, int direction, - int *return_size, int *ngrow, int **shapes) { - auto & mf = WarpX::GetInstance().getEfield_cp(lev, direction); - return getMultiFabPointers(mf, return_size, ngrow, shapes); - } - - int* warpx_getEfieldCPLoVects(int lev, int direction, - int *return_size, int *ngrow) { - auto & mf = WarpX::GetInstance().getEfield_cp(lev, direction); - return getMultiFabLoVects(mf, return_size, ngrow); - } - - double** warpx_getEfieldFP(int lev, int direction, - int *return_size, int *ngrow, int **shapes) { - auto & mf = WarpX::GetInstance().getEfield_fp(lev, direction); - return getMultiFabPointers(mf, return_size, ngrow, shapes); - } - - int* warpx_getEfieldFPLoVects(int lev, int direction, - int *return_size, int *ngrow) { - auto & mf = WarpX::GetInstance().getEfield_fp(lev, direction); - return getMultiFabLoVects(mf, return_size, ngrow); - } - - double** warpx_getBfield(int lev, int direction, - int *return_size, int *ngrow, int **shapes) { - auto & mf = WarpX::GetInstance().getBfield(lev, direction); - return getMultiFabPointers(mf, return_size, ngrow, shapes); - } - - int* warpx_getBfieldLoVects(int lev, int direction, - int *return_size, int *ngrow) { - auto & mf = WarpX::GetInstance().getBfield(lev, direction); - return getMultiFabLoVects(mf, return_size, ngrow); - } - - double** warpx_getBfieldCP(int lev, int direction, - int *return_size, int *ngrow, int **shapes) { - auto & mf = WarpX::GetInstance().getBfield_cp(lev, direction); - return getMultiFabPointers(mf, return_size, ngrow, shapes); - } - - int* warpx_getBfieldCPLoVects(int lev, int direction, - int *return_size, int *ngrow) { - auto & mf = WarpX::GetInstance().getBfield_cp(lev, direction); - return getMultiFabLoVects(mf, return_size, ngrow); - } - - double** warpx_getBfieldFP(int lev, int direction, - int *return_size, int *ngrow, int **shapes) { - auto & mf = WarpX::GetInstance().getBfield_fp(lev, direction); - return getMultiFabPointers(mf, return_size, ngrow, shapes); - } - - int* warpx_getBfieldFPLoVects(int lev, int direction, - int *return_size, int *ngrow) { - auto & mf = WarpX::GetInstance().getBfield_fp(lev, direction); - return getMultiFabLoVects(mf, return_size, ngrow); - } - - double** warpx_getCurrentDensity(int lev, int direction, - int *return_size, int *ngrow, int **shapes) { - auto & mf = WarpX::GetInstance().getcurrent(lev, direction); - return getMultiFabPointers(mf, return_size, ngrow, shapes); - } - - int* warpx_getCurrentDensityLoVects(int lev, int direction, - int *return_size, int *ngrow) { - auto & mf = WarpX::GetInstance().getcurrent(lev, direction); - return getMultiFabLoVects(mf, return_size, ngrow); - } - - double** warpx_getCurrentDensityCP(int lev, int direction, - int *return_size, int *ngrow, int **shapes) { - auto & mf = WarpX::GetInstance().getcurrent_cp(lev, direction); - return getMultiFabPointers(mf, return_size, ngrow, shapes); - } - - int* warpx_getCurrentDensityCPLoVects(int lev, int direction, - int *return_size, int *ngrow) { - auto & mf = WarpX::GetInstance().getcurrent_cp(lev, direction); - return getMultiFabLoVects(mf, return_size, ngrow); - } - - double** warpx_getCurrentDensityFP(int lev, int direction, - int *return_size, int *ngrow, int **shapes) { - auto & mf = WarpX::GetInstance().getcurrent_fp(lev, direction); - return getMultiFabPointers(mf, return_size, ngrow, shapes); - } - - int* warpx_getCurrentDensityFPLoVects(int lev, int direction, - int *return_size, int *ngrow) { - auto & mf = WarpX::GetInstance().getcurrent_fp(lev, direction); - return getMultiFabLoVects(mf, return_size, ngrow); - } - - double** warpx_getParticleStructs(int speciesnumber, - int* num_tiles, int** particles_per_tile) { - auto & mypc = WarpX::GetInstance().GetPartContainer(); - auto & myspc = mypc.GetParticleContainer(speciesnumber); - - const int level = 0; - - int i = 0; - for (WarpXParIter pti(myspc, level); pti.isValid(); ++pti, ++i) {} - - // *num_tiles = myspc.numLocalTilesAtLevel(level); - *num_tiles = i; - *particles_per_tile = (int*) malloc(*num_tiles*sizeof(int)); - - double** data = (double**) malloc(*num_tiles*sizeof(typename WarpXParticleContainer::ParticleType*)); - i = 0; - for (WarpXParIter pti(myspc, level); pti.isValid(); ++pti, ++i) { - auto& aos = pti.GetArrayOfStructs(); - data[i] = (double*) aos.data(); - (*particles_per_tile)[i] = pti.numParticles(); - } - return data; - } - - double** warpx_getParticleArrays(int speciesnumber, int comp, - int* num_tiles, int** particles_per_tile) { - auto & mypc = WarpX::GetInstance().GetPartContainer(); - auto & myspc = mypc.GetParticleContainer(speciesnumber); - - const int level = 0; - - int i = 0; - for (WarpXParIter pti(myspc, level); pti.isValid(); ++pti, ++i) {} - - // *num_tiles = myspc.numLocalTilesAtLevel(level); - *num_tiles = i; - *particles_per_tile = (int*) malloc(*num_tiles*sizeof(int)); - - double** data = (double**) malloc(*num_tiles*sizeof(double*)); - i = 0; - for (WarpXParIter pti(myspc, level); pti.isValid(); ++pti, ++i) { - auto& soa = pti.GetStructOfArrays(); - data[i] = (double*) soa.GetRealData(comp).dataPtr(); - (*particles_per_tile)[i] = pti.numParticles(); - } - return data; - } - - void warpx_ComputeDt () { - WarpX& warpx = WarpX::GetInstance(); - warpx.ComputeDt (); - } - void warpx_MoveWindow () { - WarpX& warpx = WarpX::GetInstance(); - warpx.MoveWindow (true); - } - - void warpx_EvolveE (double dt) { - WarpX& warpx = WarpX::GetInstance(); - warpx.EvolveE (dt); - } - void warpx_EvolveB (double dt) { - WarpX& warpx = WarpX::GetInstance(); - warpx.EvolveB (dt); - } - void warpx_FillBoundaryE () { - WarpX& warpx = WarpX::GetInstance(); - warpx.FillBoundaryE (); - } - void warpx_FillBoundaryB () { - WarpX& warpx = WarpX::GetInstance(); - warpx.FillBoundaryB (); - } - void warpx_SyncCurrent () { - WarpX& warpx = WarpX::GetInstance(); - warpx.SyncCurrent (); - } - void warpx_UpdateAuxilaryData () { - WarpX& warpx = WarpX::GetInstance(); - warpx.UpdateAuxilaryData (); - } - void warpx_PushParticlesandDepose (double cur_time) { - WarpX& warpx = WarpX::GetInstance(); - warpx.PushParticlesandDepose (cur_time); - } - - int warpx_getistep (int lev) { - WarpX& warpx = WarpX::GetInstance(); - return warpx.getistep (lev); - } - void warpx_setistep (int lev, int ii) { - WarpX& warpx = WarpX::GetInstance(); - warpx.setistep (lev, ii); - } - double warpx_gett_new (int lev) { - WarpX& warpx = WarpX::GetInstance(); - return warpx.gett_new (lev); - } - void warpx_sett_new (int lev, double time) { - WarpX& warpx = WarpX::GetInstance(); - warpx.sett_new (lev, time); - } - double warpx_getdt (int lev) { - WarpX& warpx = WarpX::GetInstance(); - return warpx.getdt (lev); - } - - int warpx_maxStep () { - WarpX& warpx = WarpX::GetInstance(); - return warpx.maxStep (); - } - double warpx_stopTime () { - WarpX& warpx = WarpX::GetInstance(); - return warpx.stopTime (); - } - - int warpx_checkInt () { - WarpX& warpx = WarpX::GetInstance(); - return warpx.checkInt (); - } - int warpx_plotInt () { - WarpX& warpx = WarpX::GetInstance(); - return warpx.plotInt (); - } - - void warpx_WriteCheckPointFile () { - WarpX& warpx = WarpX::GetInstance(); - warpx.WriteCheckPointFile (); - } - void warpx_WritePlotFile () { - WarpX& warpx = WarpX::GetInstance(); - warpx.WritePlotFile (); - } - - int warpx_finestLevel () { - WarpX& warpx = WarpX::GetInstance(); - return warpx.finestLevel (); - } - - void mypc_Redistribute () { - auto & mypc = WarpX::GetInstance().GetPartContainer(); - mypc.Redistribute(); - } - -} - diff --git a/Source/Utils/WarpXWrappers.h b/Source/Utils/WarpXWrappers.h deleted file mode 100644 index 07d6f80f7..000000000 --- a/Source/Utils/WarpXWrappers.h +++ /dev/null @@ -1,120 +0,0 @@ -#ifndef WARPX_WRAPPERS_H_ -#define WARPX_WRAPPERS_H_ - -#ifdef BL_USE_MPI -#include -#endif - -#ifdef __cplusplus -extern "C" { -#endif - - int warpx_nSpecies(); - - bool warpx_use_fdtd_nci_corr(); - - int warpx_l_lower_order_in_v(); - - int warpx_nComps(); - - int warpx_SpaceDim(); - - void amrex_init (int argc, char* argv[]); - -#ifdef BL_USE_MPI - void amrex_init_with_inited_mpi (int argc, char* argv[], MPI_Comm mpicomm); -#endif - - void amrex_finalize (int finalize_mpi); - - void warpx_init (); - - void warpx_finalize (); - - typedef void(*WARPX_CALLBACK_PY_FUNC_0)(); - - void warpx_set_callback_py_afterinit (WARPX_CALLBACK_PY_FUNC_0); - void warpx_set_callback_py_beforeEsolve (WARPX_CALLBACK_PY_FUNC_0); - void warpx_set_callback_py_afterEsolve (WARPX_CALLBACK_PY_FUNC_0); - void warpx_set_callback_py_beforedeposition (WARPX_CALLBACK_PY_FUNC_0); - void warpx_set_callback_py_afterdeposition (WARPX_CALLBACK_PY_FUNC_0); - void warpx_set_callback_py_particlescraper (WARPX_CALLBACK_PY_FUNC_0); - void warpx_set_callback_py_particleloader (WARPX_CALLBACK_PY_FUNC_0); - void warpx_set_callback_py_beforestep (WARPX_CALLBACK_PY_FUNC_0); - void warpx_set_callback_py_afterstep (WARPX_CALLBACK_PY_FUNC_0); - void warpx_set_callback_py_afterrestart (WARPX_CALLBACK_PY_FUNC_0); - void warpx_set_callback_py_particleinjection (WARPX_CALLBACK_PY_FUNC_0); - void warpx_set_callback_py_appliedfields (WARPX_CALLBACK_PY_FUNC_0); - - void warpx_evolve (int numsteps); // -1 means the inputs parameter will be used. - - void warpx_addNParticles(int speciesnumber, int lenx, - double* x, double* y, double* z, - double* vx, double* vy, double* vz, - int nattr, double* attr, int uniqueparticles); - - double warpx_getProbLo(int dir); - - double warpx_getProbHi(int dir); - - long warpx_getNumParticles(int speciesnumber); - - double** warpx_getEfield(int lev, int direction, - int *return_size, int* ngrow, int **shapes); - - int* warpx_getEfieldLoVects(int lev, int direction, - int *return_size, int* ngrow); - - double** warpx_getBfield(int lev, int direction, - int *return_size, int* ngrow, int **shapes); - - int* warpx_getBfieldLoVects(int lev, int direction, - int *return_size, int* ngrow); - - double** warpx_getCurrentDensity(int lev, int direction, - int *return_size, int* ngrow, int **shapes); - - int* warpx_getCurrentDensityLoVects(int lev, int direction, - int *return_size, int* ngrow); - - double** warpx_getParticleStructs(int speciesnumber, - int* num_tiles, int** particles_per_tile); - - double** warpx_getParticleArrays(int speciesnumber, int comp, - int* num_tiles, int** particles_per_tile); - - void warpx_ComputeDt (); - void warpx_MoveWindow (); - - void warpx_EvolveE (double dt); - void warpx_EvolveB (double dt); - void warpx_FillBoundaryE (); - void warpx_FillBoundaryB (); - void warpx_SyncCurrent (); - void warpx_UpdateAuxilaryData (); - void warpx_PushParticlesandDepose (double cur_time); - - int warpx_getistep (int lev); - void warpx_setistep (int lev, int ii); - double warpx_gett_new (int lev); - void warpx_sett_new (int lev, double time); - double warpx_getdt (int lev); - - int warpx_maxStep (); - double warpx_stopTime (); - - int warpx_checkInt (); - int warpx_plotInt (); - - void warpx_WriteCheckPointFile (); - void warpx_WritePlotFile (); - - int warpx_finestLevel (); - - void mypc_Redistribute (); - -#ifdef __cplusplus -} -#endif - -#endif -- cgit v1.2.3