diff options
Diffstat (limited to 'Source')
17 files changed, 515 insertions, 583 deletions
diff --git a/Source/Particles/ElementaryProcess/Ionization.H b/Source/Particles/ElementaryProcess/Ionization.H new file mode 100644 index 000000000..89e957fc6 --- /dev/null +++ b/Source/Particles/ElementaryProcess/Ionization.H @@ -0,0 +1,73 @@ +#ifndef IONIZATION_H_ +#define IONIZATION_H_ + +#include "WarpXConst.H" +#include "WarpXParticleContainer.H" + +struct IonizationFilterFunc +{ + const amrex::Real* const AMREX_RESTRICT m_ionization_energies; + const amrex::Real* const AMREX_RESTRICT m_adk_prefactor; + const amrex::Real* const AMREX_RESTRICT m_adk_exp_prefactor; + const amrex::Real* const AMREX_RESTRICT m_adk_power; + + int comp; + int m_atomic_number; + + template <typename PData> + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + bool operator() (const PData& ptd, int i) const noexcept + { + int ion_lev = ptd.m_runtime_idata[comp][i]; + if (ion_lev < m_atomic_number) + { + amrex::Real c = PhysConst::c; + amrex::Real c2_inv = 1./c/c; + + // Compute electric field amplitude in the particle's frame of + // reference (particularly important when in boosted frame). + amrex::ParticleReal ux = ptd.m_rdata[PIdx::ux][i]; + amrex::ParticleReal uy = ptd.m_rdata[PIdx::uy][i]; + amrex::ParticleReal uz = ptd.m_rdata[PIdx::uz][i]; + amrex::ParticleReal ex = ptd.m_rdata[PIdx::Ex][i]; + amrex::ParticleReal ey = ptd.m_rdata[PIdx::Ey][i]; + amrex::ParticleReal ez = ptd.m_rdata[PIdx::Ez][i]; + amrex::ParticleReal bx = ptd.m_rdata[PIdx::Bx][i]; + amrex::ParticleReal by = ptd.m_rdata[PIdx::By][i]; + amrex::ParticleReal bz = ptd.m_rdata[PIdx::Bz][i]; + + amrex::Real ga = std::sqrt(1. + (ux*ux + uy*uy + uz*uz) * c2_inv); + amrex::Real E = std::sqrt( + - ( ux*ex + uy*ey + uz*ez ) * ( ux*ex + uy*ey + uz*ez ) * c2_inv + + ( ga *ex + uy*bz - uz*by ) * ( ga *ex + uy*bz - uz*by ) + + ( ga *ey + uz*bx - ux*bz ) * ( ga *ey + uz*bx - ux*bz ) + + ( ga *ez + ux*by - uy*bx ) * ( ga *ez + ux*by - uy*bx ) + ); + + // Compute probability of ionization p + amrex::Real w_dtau = 1./ ga * m_adk_prefactor[ion_lev] * + std::pow(E, m_adk_power[ion_lev]) * + std::exp( m_adk_exp_prefactor[ion_lev]/E ); + amrex::Real p = 1. - std::exp( - w_dtau ); + + amrex::Real random_draw = amrex::Random(); + if (random_draw < p) + { + return true; + } + } + return false; + } +}; + +struct IonizationTransformFunc +{ + template <typename DstData, typename SrcData> + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + void operator() (DstData& dst, SrcData& src, int i_src, int i_dst) const noexcept + { + src.m_runtime_idata[0][i_src] += 1; + } +}; + +#endif diff --git a/Source/Particles/ElementaryProcess/Make.package b/Source/Particles/ElementaryProcess/Make.package new file mode 100644 index 000000000..8caa7913a --- /dev/null +++ b/Source/Particles/ElementaryProcess/Make.package @@ -0,0 +1,4 @@ +CEXE_headers += Ionization.H + +INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Particles/ElementaryProcess/ +VPATH_LOCATIONS += $(WARPX_HOME)/Source/Particles/ElementaryProcess/ diff --git a/Source/Particles/Make.package b/Source/Particles/Make.package index 8e336ba69..3bd4829a0 100644 --- a/Source/Particles/Make.package +++ b/Source/Particles/Make.package @@ -19,6 +19,7 @@ include $(WARPX_HOME)/Source/Particles/Deposition/Make.package include $(WARPX_HOME)/Source/Particles/Gather/Make.package include $(WARPX_HOME)/Source/Particles/Sorting/Make.package include $(WARPX_HOME)/Source/Particles/ParticleCreation/Make.package +include $(WARPX_HOME)/Source/Particles/ElementaryProcess/Make.package include $(WARPX_HOME)/Source/Particles/Collision/Make.package INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Particles diff --git a/Source/Particles/MultiParticleContainer.H b/Source/Particles/MultiParticleContainer.H index 65c13e39b..767abbeeb 100644 --- a/Source/Particles/MultiParticleContainer.H +++ b/Source/Particles/MultiParticleContainer.H @@ -11,14 +11,14 @@ #ifndef WARPX_ParticleContainer_H_ #define WARPX_ParticleContainer_H_ -#include "ElementaryProcess.H" - #include <WarpXParticleContainer.H> #include <PhysicalParticleContainer.H> #include <RigidInjectedParticleContainer.H> #include <PhotonParticleContainer.H> #include <LaserParticleContainer.H> #include <WarpXParserWrapper.H> +#include <SmartCopy.H> +#include <FilterCopyTransform.H> #include <AMReX_Particles.H> #ifdef WARPX_QED @@ -217,8 +217,6 @@ public: PhysicalParticleContainer& GetPCtmp () { return *pc_tmp; } - IonizationProcess ionization_process; - std::string m_B_ext_particle_s = "default"; std::string m_E_ext_particle_s = "default"; // External fields added to particle fields. @@ -255,6 +253,9 @@ protected: std::vector<PCTypes> species_types; + amrex::MFItInfo getMFItInfo (const WarpXParticleContainer& pc_src, + const WarpXParticleContainer& pc_dst) const noexcept; + #ifdef WARPX_QED // The QED engines std::shared_ptr<BreitWheelerEngine> m_shr_p_bw_engine; diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index f9a0d230b..db47c7fa2 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -63,7 +63,6 @@ MultiParticleContainer::MultiParticleContainer (AmrCore* amr_core) nspecies_back_transformed_diagnostics += 1; } } - ionization_process = IonizationProcess(); // collision allcollisions.resize(ncollisions); @@ -642,77 +641,51 @@ void MultiParticleContainer::doFieldIonization () { BL_PROFILE("MPC::doFieldIonization"); + // Loop over all species. // Ionized particles in pc_source create particles in pc_product - for (auto& pc_source : allcontainers){ - - // Skip if not ionizable + for (auto& pc_source : allcontainers) + { if (!pc_source->do_field_ionization){ continue; } - - // Get product species + auto& pc_product = allcontainers[pc_source->ionization_product]; - for (int lev = 0; lev <= pc_source->finestLevel(); ++lev){ + SmartCopyFactory copy_factory(*pc_source, *pc_product); + auto phys_pc_ptr = static_cast<PhysicalParticleContainer*>(pc_source.get()); - // When using runtime components, AMReX requires to touch all tiles - // in serial and create particles tiles with runtime components if - // they do not exist (or if they were defined by default, i.e., - // without runtime component). -#ifdef _OPENMP - // Touch all tiles of source species in serial if runtime attribs - for (MFIter mfi = pc_source->MakeMFIter(lev); mfi.isValid(); ++mfi) { - const int grid_id = mfi.index(); - const int tile_id = mfi.LocalTileIndex(); - pc_source->GetParticles(lev)[std::make_pair(grid_id,tile_id)]; - if ( (pc_source->NumRuntimeRealComps()>0) || (pc_source->NumRuntimeIntComps()>0) ) { - pc_source->DefineAndReturnParticleTile(lev, grid_id, tile_id); - } - } -#endif - // Touch all tiles of product species in serial - for (MFIter mfi = pc_source->MakeMFIter(lev); mfi.isValid(); ++mfi) { - const int grid_id = mfi.index(); - const int tile_id = mfi.LocalTileIndex(); - pc_product->GetParticles(lev)[std::make_pair(grid_id,tile_id)]; - pc_product->DefineAndReturnParticleTile(lev, grid_id, tile_id); - } + auto Filter = phys_pc_ptr->getIonizationFunc(); + auto Copy = copy_factory.getSmartCopy(); + auto Transform = IonizationTransformFunc(); - // Enable tiling - MFItInfo info; - if (pc_source->do_tiling && Gpu::notInLaunchRegion()) { - AMREX_ALWAYS_ASSERT_WITH_MESSAGE( - pc_product->do_tiling, - "For ionization, either all or none of the " - "particle species must use tiling."); - info.EnableTiling(pc_source->tile_size); - } + pc_source ->defineAllParticleTiles(); + pc_product->defineAllParticleTiles(); + + for (int lev = 0; lev <= pc_source->finestLevel(); ++lev) + { + auto info = getMFItInfo(*pc_source, *pc_product); #ifdef _OPENMP - info.SetDynamic(true); -#pragma omp parallel -#endif - // Loop over all grids (if not tiling) or grids and tiles (if tiling) +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif for (MFIter mfi = pc_source->MakeMFIter(lev, info); mfi.isValid(); ++mfi) { - // Ionization mask: one element per source particles. - // 0 if not ionized, 1 if ionized. - amrex::Gpu::ManagedDeviceVector<int> is_ionized; - pc_source->buildIonizationMask(mfi, lev, is_ionized); - // Create particles in pc_product - int do_boost = WarpX::do_back_transformed_diagnostics - && pc_product->doBackTransformedDiagnostics(); - amrex::Gpu::ManagedDeviceVector<int> v_do_back_transformed_product{do_boost}; - const amrex::Vector<WarpXParticleContainer*> v_pc_product {pc_product.get()}; - // Copy source to product particles, and increase ionization - // level of source particle - ionization_process.createParticles(lev, mfi, pc_source, v_pc_product, - is_ionized, v_do_back_transformed_product); - // Synchronize to prevent the destruction of temporary arrays (at the - // end of the function call) before the kernel executes. - Gpu::streamSynchronize(); + auto& src_tile = pc_source ->ParticlesAt(lev, mfi); + auto& dst_tile = pc_product->ParticlesAt(lev, mfi); + + auto np_src = src_tile.numParticles(); + auto np_dst = dst_tile.numParticles(); + + dst_tile.resize(np_dst + np_src); + + auto num_added = filterCopyTransformParticles(dst_tile, src_tile, np_dst, + Filter, Copy, Transform); + + dst_tile.resize(np_dst + num_added); + + setNewParticleIDs(dst_tile, np_dst, num_added); } - } // lev - } // pc_source + } + } } void @@ -735,7 +708,7 @@ MultiParticleContainer::doCoulombCollisions () // Loop over all grids/tiles at this level #ifdef _OPENMP info.SetDynamic(true); - #pragma omp parallel if (Gpu::notInLaunchRegion()) +#pragma omp parallel if (Gpu::notInLaunchRegion()) #endif for (MFIter mfi = species1->MakeMFIter(lev, info); mfi.isValid(); ++mfi){ @@ -749,6 +722,25 @@ MultiParticleContainer::doCoulombCollisions () } } +MFItInfo MultiParticleContainer::getMFItInfo (const WarpXParticleContainer& pc_src, + const WarpXParticleContainer& pc_dst) const noexcept +{ + MFItInfo info; + + if (pc_src.do_tiling && Gpu::notInLaunchRegion()) { + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(pc_dst.do_tiling, + "For ionization, either all or none of the " + "particle species must use tiling."); + info.EnableTiling(pc_src.tile_size); + } + +#ifdef _OPENMP + info.SetDynamic(true); +#endif + + return info; +} + #ifdef WARPX_QED void MultiParticleContainer::InitQED () { diff --git a/Source/Particles/ParticleCreation/CopyParticle.H b/Source/Particles/ParticleCreation/CopyParticle.H deleted file mode 100644 index 8b2770891..000000000 --- a/Source/Particles/ParticleCreation/CopyParticle.H +++ /dev/null @@ -1,96 +0,0 @@ -/* Copyright 2019 Axel Huebl, Maxence Thevenet - * - * This file is part of WarpX. - * - * License: BSD-3-Clause-LBNL - */ -#ifndef COPYPARTICLE_H_ -#define COPYPARTICLE_H_ - -#include "WarpXParticleContainer.H" - -/** - * \brief Functor to copy one particle - * - * This is meant to be a very small class captured by value in kernel launches, - * that can be initialized on the host and copied to the device at each - * iteration. It should be general enough to be used by all processes. - * - * Pointers to SoA data are saved when constructor is called. - * AoS data is passed as argument of operator (). - */ -class CopyParticle -{ - -public: - - // ID of MPI rank - int m_cpuid; - // If true, will copy old attribs for back-transforme diagnostics - bool m_do_back_transformed_product; - // Source old (runtime) attribs for back-transformed diagnostics - amrex::GpuArray<amrex::ParticleReal*,3> m_runtime_uold_source; - // Source attribs - amrex::GpuArray<amrex::ParticleReal*,PIdx::nattribs> m_attribs_source; - // Product attribs - amrex::GpuArray<amrex::ParticleReal*,PIdx::nattribs> m_attribs_product; - // Product runtime attribs - amrex::GpuArray<amrex::ParticleReal*,6> m_runtime_attribs_product; - - // Simple constructor - AMREX_GPU_HOST_DEVICE - CopyParticle( - const int cpuid, const int do_back_transformed_product, - const amrex::GpuArray<amrex::ParticleReal*,3> runtime_uold_source, - const amrex::GpuArray<amrex::ParticleReal*,PIdx::nattribs> attribs_source, - const amrex::GpuArray<amrex::ParticleReal*,PIdx::nattribs> attribs_product, - const amrex::GpuArray<amrex::ParticleReal*,6> runtime_attribs_product) - : - m_cpuid(cpuid), - m_do_back_transformed_product(do_back_transformed_product), - m_runtime_uold_source(runtime_uold_source), - m_attribs_source(attribs_source), - m_attribs_product(attribs_product), - m_runtime_attribs_product(runtime_attribs_product){} - - /** - * \brief Overload operator () to do the copy of 1 particle - * - * \param is index of source particle - * \param ip index of product particle - * \param pid_product ID of first product particle - * \param p_source Struct with data for 1 source particle (position etc.) - * \param p_product Struct with data for 1 product particle (position etc.) - */ - AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE - void operator () (int is, int ip, int pid_product, - WarpXParticleContainer::ParticleType& p_source, - WarpXParticleContainer::ParticleType& p_product) const noexcept - { - // Copy particle from source to product: AoS - p_product.id() = pid_product + ip; - p_product.cpu() = m_cpuid; - p_product.pos(0) = p_source.pos(0); - p_product.pos(1) = p_source.pos(1); -#if (AMREX_SPACEDIM == 3) - p_product.pos(2) = p_source.pos(2); -#endif - // Copy particle from source to product: SoA - for (int ia = 0; ia < PIdx::nattribs; ++ia) { - m_attribs_product[ia][ip] = m_attribs_source[ia][is]; - } - // Update xold etc. if boosted frame diagnostics required - // for product species. Fill runtime attribs with a copy of - // current properties (xold = x etc.). - if (m_do_back_transformed_product) { - m_runtime_attribs_product[0][ip] = p_source.pos(0); - m_runtime_attribs_product[1][ip] = p_source.pos(1); - m_runtime_attribs_product[2][ip] = p_source.pos(2); - m_runtime_attribs_product[3][ip] = m_runtime_uold_source[0][ip]; - m_runtime_attribs_product[4][ip] = m_runtime_uold_source[1][ip]; - m_runtime_attribs_product[5][ip] = m_runtime_uold_source[2][ip]; - } - } -}; - -#endif // COPYPARTICLE_H_ diff --git a/Source/Particles/ParticleCreation/DefaultInitialization.H b/Source/Particles/ParticleCreation/DefaultInitialization.H new file mode 100644 index 000000000..f2f0dd3d8 --- /dev/null +++ b/Source/Particles/ParticleCreation/DefaultInitialization.H @@ -0,0 +1,86 @@ +#ifndef DEFAULTINITIALIZATION_H_ +#define DEFAULTINITIALIZATION_H_ + +#include "AMReX_REAL.H" +#include "AMReX_GpuContainers.H" + +#include <map> +#include <string> + +/** + * \brief This set of initialization policies describes what happens + * when we need to create a new particle due to an elementary process. + * For example, when an ionization event creates an electron, these + * policies control the initial values of the electron's components. + * These can always be over-written later. + * + * The specific meanings are as follows: + * Zero - set the component to zero + * One - set the component to one + * SignalingNan - set the component to a signalling nan. If you run + * with amrex.fpe_trap_invalid=1 AND built the code + * with gcc, this will throw an exception the first time + * it is used. + * RandomTau - a special flag for the "tau" component used by certain + * QED processes, which gets a random initial value. + * + */ +enum struct InitializationPolicy {Zero=0, One, SignalingNan, RandomTau}; + +/** + * \brief This map sets the initialization policy for each particle component + * used in WarpX. + */ +static std::map<std::string, InitializationPolicy> initialization_policies = { + {"w", InitializationPolicy::Zero }, + {"ux", InitializationPolicy::Zero }, + {"uy", InitializationPolicy::Zero }, + {"uz", InitializationPolicy::Zero }, + {"Ex", InitializationPolicy::Zero }, + {"Ey", InitializationPolicy::Zero }, + {"Ez", InitializationPolicy::Zero }, + {"Bx", InitializationPolicy::Zero }, + {"By", InitializationPolicy::Zero }, + {"Bz", InitializationPolicy::Zero }, +#ifdef WARPX_DIM_RZ + {"theta", InitializationPolicy::Zero}, +#endif + {"tau", InitializationPolicy::RandomTau} +}; + +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +amrex::ParticleReal initializeRealValue (const InitializationPolicy policy) noexcept +{ + switch (policy) { + case InitializationPolicy::Zero : return 0.0; + case InitializationPolicy::One : return 1.0; + case InitializationPolicy::SignalingNan : { + return 1.0; + } + case InitializationPolicy::RandomTau : { + return 1.0; + } + default : { + amrex::Abort("Initialization Policy not recognized"); + return 1.0; + } + } +} + +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +int initializeIntValue (const InitializationPolicy policy) noexcept +{ + switch (policy) { + case InitializationPolicy::Zero : return 0; + case InitializationPolicy::One : return 1; + case InitializationPolicy::SignalingNan : { + return 1; + } + default : { + amrex::Abort("Initialization Policy not recognized"); + return 1; + } + } +} + +#endif diff --git a/Source/Particles/ParticleCreation/ElementaryProcess.H b/Source/Particles/ParticleCreation/ElementaryProcess.H deleted file mode 100644 index 3fe2240cc..000000000 --- a/Source/Particles/ParticleCreation/ElementaryProcess.H +++ /dev/null @@ -1,269 +0,0 @@ -/* Copyright 2019 Axel Huebl, Maxence Thevenet, Weiqun Zhang - * - * - * This file is part of WarpX. - * - * License: BSD-3-Clause-LBNL - */ -#ifndef ELEMENTARYPROCESS_H_ -#define ELEMENTARYPROCESS_H_ - -#include "WarpXParticleContainer.H" -#include "CopyParticle.H" -#include "TransformParticle.H" - -/** - * \brief Base class for particle creation processes - * - * Particles in a source particle container are copied to product particle - * containers if they are flagged. Both source and product particles can be - * modified. - * - * - initialize_copy initializes a general copy functor - * - createParticles formats the data to prepare for the copy, then - * calls copyAndTransformParticles - * - copyAndTransformParticles loops over particles, performs the copy and - * transform source and product particles if needed. - * - * The class is templated on the process type. This gives flexibility - * on source and product particle transformations. - */ -template<ElementaryProcessType ProcessT> -class ElementaryProcess -{ -public: - - /** - * \brief initialize and return functor for particle copy - * - * \param cpuid id of MPI rank - * \param do_back_transformed_product whether to copy old attribs - * \param runtime_uold_source momentum of source particles - * \param attribs_source attribs of source particles - * \param attribs_product attribs of product particles - * \param runtime_attribs_product runtime attribs for product, to store old attribs - */ - CopyParticle initialize_copy( - const int cpuid, const int do_back_transformed_product, - const amrex::GpuArray<amrex::ParticleReal*,3> runtime_uold_source, - const amrex::GpuArray<amrex::ParticleReal*,PIdx::nattribs> attribs_source, - const amrex::GpuArray<amrex::ParticleReal*,PIdx::nattribs> attribs_product, - const amrex::GpuArray<amrex::ParticleReal*,6> runtime_attribs_product) const noexcept - { - return CopyParticle( - cpuid, do_back_transformed_product, runtime_uold_source, - attribs_source, attribs_product, runtime_attribs_product - ); - }; - - /** - * \brief create particles in product particle containers - * - * For particle i in mfi, if is_flagged[i]=1, copy particle - * particle i from container pc_source into each container in v_pc_product - * - * \param lev MR level - * \param mfi MFIter where source particles are located - * \param pc_source source particle container - * \param v_pc_product vector of product particle containers - * \param is_flagged particles for which is_flagged is 1 are copied - * \param v_do_back_transformed_product vector: whether to copy old attribs for - * each product container - */ - void createParticles ( - int lev, const amrex::MFIter& mfi, - std::unique_ptr< WarpXParticleContainer>& pc_source, - amrex::Vector<WarpXParticleContainer*> v_pc_product, - amrex::Gpu::ManagedDeviceVector<int>& is_flagged, - amrex::Gpu::ManagedDeviceVector<int> v_do_back_transformed_product) - { - - BL_PROFILE("createParticles"); - const int grid_id = mfi.index(); - const int tile_id = mfi.LocalTileIndex(); - - // Get source particle data - auto& ptile_source = pc_source->GetParticles(lev)[std::make_pair(grid_id,tile_id)]; - const int np_source = ptile_source.GetArrayOfStructs().size(); - if (np_source == 0) return; - // --- source AoS particle data - WarpXParticleContainer::ParticleType* particles_source = ptile_source.GetArrayOfStructs()().data(); - // --- source SoA particle data - auto& soa_source = ptile_source.GetStructOfArrays(); - amrex::GpuArray<amrex::ParticleReal*,PIdx::nattribs> attribs_source; - for (int ia = 0; ia < PIdx::nattribs; ++ia) { - attribs_source[ia] = soa_source.GetRealData(ia).data(); - } - // --- source runtime attribs - amrex::GpuArray<amrex::ParticleReal*,3> runtime_uold_source; - // Prepare arrays for boosted frame diagnostics. - runtime_uold_source[0] = soa_source.GetRealData(PIdx::ux).data(); - runtime_uold_source[1] = soa_source.GetRealData(PIdx::uy).data(); - runtime_uold_source[2] = soa_source.GetRealData(PIdx::uz).data(); - - // --- source runtime integer attribs - amrex::GpuArray<int*,1> runtime_iattribs_source; - std::map<std::string, int> icomps_source = pc_source->getParticleiComps(); - runtime_iattribs_source[0] = soa_source.GetIntData(icomps_source["ionization_level"]).data(); - - int nproducts = v_pc_product.size(); - AMREX_ALWAYS_ASSERT_WITH_MESSAGE( - v_do_back_transformed_product.size() == nproducts, - "v_pc_product and v_do_back_transformed_product must have the same size."); - - // Indices of product particle for each source particle. - // i_product[i]-1 is the location in product tile of product particle - // from source particle i. - amrex::Gpu::ManagedDeviceVector<int> i_product; - i_product.resize(np_source); - // 0<i<np_source - // 0<i_product<np_flagged - // Strictly speaking, i_product should be an exclusive_scan of - // is_flagged. However, for indices where is_flagged is 1, the - // inclusive scan gives the same result with an offset of 1. - // The advantage of inclusive_scan is that the sum of is_flagged - // is in the last element, so no other reduction is required to get - // number of particles. - // Gpu::inclusive_scan runs on the current GPU stream, and synchronizes - // with the CPU, so that the next line (executed by the CPU) has the - // updated values of i_product - amrex::Gpu::inclusive_scan(is_flagged.begin(), is_flagged.end(), i_product.begin()); - int np_flagged = i_product[np_source-1]; - if (np_flagged == 0) return; - - amrex::Gpu::ManagedDeviceVector<CopyParticle> v_copy_functor; - v_copy_functor.reserve(nproducts); - amrex::Gpu::ManagedDeviceVector<int> v_pid_product(nproducts); - amrex::Gpu::ManagedDeviceVector<WarpXParticleContainer::ParticleType*> v_particles_product(nproducts); - for (int iproduct=0; iproduct<nproducts; iproduct++){ - WarpXParticleContainer*& pc_product = v_pc_product[iproduct]; - // Get product particle data - auto& ptile_product = pc_product->GetParticles(lev)[std::make_pair(grid_id,tile_id)]; - // old and new (i.e., including flagged particles) number of particles - // for product species - const int np_product_old = ptile_product.GetArrayOfStructs().size(); - const int np_product_new = np_product_old + np_flagged; - // Allocate extra space in product species for flagged particles. - ptile_product.resize(np_product_new); - // --- product AoS particle data - // WarpXParticleContainer::ParticleType* particles_product = ptile_product.GetArrayOfStructs()().data() + np_product_old; - v_particles_product[iproduct] = ptile_product.GetArrayOfStructs()().data() + np_product_old; - // First element is the first newly-created product particle - // --- product SoA particle data - auto& soa_product = ptile_product.GetStructOfArrays(); - amrex::GpuArray<amrex::ParticleReal*,PIdx::nattribs> attribs_product; - for (int ia = 0; ia < PIdx::nattribs; ++ia) { - // First element is the first newly-created product particle - attribs_product[ia] = soa_product.GetRealData(ia).data() + np_product_old; - } - // --- product runtime attribs - amrex::GpuArray<amrex::ParticleReal*,6> runtime_attribs_product; - const int do_boost = v_do_back_transformed_product[iproduct]; - if (do_boost) { - std::map<std::string, int> comps_product = pc_product->getParticleComps(); - runtime_attribs_product[0] = soa_product.GetRealData(comps_product[ "xold"]).data() + np_product_old; - runtime_attribs_product[1] = soa_product.GetRealData(comps_product[ "yold"]).data() + np_product_old; - runtime_attribs_product[2] = soa_product.GetRealData(comps_product[ "zold"]).data() + np_product_old; - runtime_attribs_product[3] = soa_product.GetRealData(comps_product["uxold"]).data() + np_product_old; - runtime_attribs_product[4] = soa_product.GetRealData(comps_product["uyold"]).data() + np_product_old; - runtime_attribs_product[5] = soa_product.GetRealData(comps_product["uzold"]).data() + np_product_old; - } - - // --- product runtime integer attribs - int pid_product; -#ifdef _OPENMP -#pragma omp critical (doParticleCreation_nextid) -#endif - { - // ID of first newly-created product particle - pid_product = pc_product->NextID(); - // Update NextID to include particles created in this function - pc_product->setNextID(pid_product+np_flagged); - } - const int cpuid = amrex::ParallelDescriptor::MyProc(); - - // Create instance of copy functor, and add it to the vector - v_copy_functor.push_back( initialize_copy( - cpuid, v_do_back_transformed_product[iproduct], - runtime_uold_source, - attribs_source, - attribs_product, - runtime_attribs_product - ) ); - v_pid_product[iproduct] = pid_product; - - } - // Loop over source particles and create product particles - copyAndTransformParticles(is_flagged, i_product, np_source, v_pid_product, - v_particles_product, particles_source, v_copy_functor, - runtime_iattribs_source); - // Synchronize to prevent the destruction of temporary arrays (at the - // end of the function call) before the kernel executes. - amrex::Gpu::streamSynchronize(); - } - - /** - * \brief Loop over source particles and create product particles - * - * \param is_flagged particles for which is_flagged is 1 are copied - * \param i_product relative indices of product particle. This is the same - * for all product particle containers. - * \param np_source number of particles in source container - * \param v_pid_product for each product particle container: ID of the - * first product particle. Others IDs are incremented from this one - * \param v_particles_product for each product particle container: pointer - * to the AoS particle data - * \param particles_source pointter to the AoS source particle data - * \param v_copy_functor vector of copy functors, one for each product particle container - * \param runtime_iattribs_source pointer to source runtime integer attribs - */ - void copyAndTransformParticles( - amrex::Gpu::ManagedDeviceVector<int>& is_flagged, - amrex::Gpu::ManagedDeviceVector<int>& i_product, - int np_source, - amrex::Gpu::ManagedDeviceVector<int> v_pid_product, - amrex::Gpu::ManagedDeviceVector<WarpXParticleContainer::ParticleType*> v_particles_product, - WarpXParticleContainer::ParticleType* particles_source, - const amrex::Gpu::ManagedDeviceVector<CopyParticle>& v_copy_functor, - amrex::GpuArray<int*,1> runtime_iattribs_source) - { - int const * const AMREX_RESTRICT p_is_flagged = is_flagged.dataPtr(); - int const * const AMREX_RESTRICT p_i_product = i_product.dataPtr(); - CopyParticle const * const AMREX_RESTRICT p_copy_functor = v_copy_functor.dataPtr(); - WarpXParticleContainer::ParticleType * * p_particles_product = v_particles_product.data(); - int const * const p_pid_product = v_pid_product.data(); - - // Loop over all source particles. If is_flagged, copy particle data - // to corresponding product particle. - int nproducts = v_particles_product.size(); - amrex::For( - np_source, [=] AMREX_GPU_DEVICE (int is) noexcept - { - if(p_is_flagged[is]){ - // offset of 1 due to inclusive scan - int ip = p_i_product[is]-1; - WarpXParticleContainer::ParticleType& p_source = particles_source[is]; - for (int iproduct=0; iproduct<nproducts; iproduct++){ - // is: index of flagged particle in source species - // ip: index of corresponding new particle in product species - WarpXParticleContainer::ParticleType* particles_product = p_particles_product[iproduct]; - WarpXParticleContainer::ParticleType& p_product = particles_product[ip]; - p_copy_functor[iproduct](is, ip, p_pid_product[iproduct], p_source, p_product); - } - // Transform source particle - transformSourceParticle<ProcessT>(is, p_source, runtime_iattribs_source); - // Transform product particles. The loop over product - // species is done inside the function to allow for - // more flexibility. - transformProductParticle<ProcessT>(ip, p_particles_product); - } - } - ); - } -}; - -// Derived class for ionization process -class IonizationProcess: public ElementaryProcess<ElementaryProcessType::Ionization> -{}; - -#endif // ELEMENTARYPROCESS_H_ diff --git a/Source/Particles/ParticleCreation/FilterCopyTransform.H b/Source/Particles/ParticleCreation/FilterCopyTransform.H new file mode 100644 index 000000000..6e47dd732 --- /dev/null +++ b/Source/Particles/ParticleCreation/FilterCopyTransform.H @@ -0,0 +1,62 @@ +#ifndef FILTER_COPY_TRANSFORM_H_ +#define FILTER_COPY_TRANSFORM_H_ + +#include <AMReX_GpuContainers.H> +#include <AMReX_TypeTraits.H> + +template <typename DstTile, typename SrcTile, typename Index, + typename TransFunc, typename CopyFunc, + amrex::EnableIf_t<std::is_integral<Index>::value, int> foo = 0> +Index filterCopyTransformParticles (DstTile& dst, SrcTile& src, Index* mask, Index dst_index, + CopyFunc&& copy, TransFunc&& transform) noexcept +{ + using namespace amrex; + auto np = src.numParticles(); + Gpu::DeviceVector<Index> offsets(np); + Gpu::exclusive_scan(mask, mask+np, offsets.begin()); + + Index last_mask, last_offset; + Gpu::copyAsync(Gpu::deviceToHost, mask+np-1, mask + np, &last_mask); + Gpu::copyAsync(Gpu::deviceToHost, offsets.data()+np-1, offsets.data()+np, &last_offset); + + auto p_offsets = offsets.dataPtr(); + + const auto src_data = src.getParticleTileData(); + auto dst_data = dst.getParticleTileData(); + + AMREX_HOST_DEVICE_FOR_1D( np, i, + { + if (mask[i]) + { + copy(dst_data, src_data, i, p_offsets[i] + dst_index); + transform(dst_data, src_data, i, p_offsets[i] + dst_index); + } + }); + + Gpu::streamSynchronize(); + return last_mask + last_offset; +} + +template <typename DstTile, typename SrcTile, typename Index, + typename PredFunc, typename TransFunc, typename CopyFunc> +Index filterCopyTransformParticles (DstTile& dst, SrcTile& src, Index dst_index, + PredFunc&& filter, CopyFunc&& copy, TransFunc&& transform) noexcept +{ + using namespace amrex; + auto np = src.numParticles(); + Gpu::DeviceVector<Index> mask(np); + + auto p_mask = mask.dataPtr(); + const auto src_data = src.getParticleTileData(); + + AMREX_HOST_DEVICE_FOR_1D(np, i, + { + p_mask[i] = filter(src_data, i); + }); + + return filterCopyTransformParticles(dst, src, mask.dataPtr(), dst_index, + std::forward<CopyFunc>(copy), + std::forward<TransFunc>(transform)); +} + +#endif diff --git a/Source/Particles/ParticleCreation/Make.package b/Source/Particles/ParticleCreation/Make.package index 6e32f4a77..32c0add64 100644 --- a/Source/Particles/ParticleCreation/Make.package +++ b/Source/Particles/ParticleCreation/Make.package @@ -1,6 +1,7 @@ -CEXE_headers += ElementaryProcess.H -CEXE_headers += CopyParticle.H -CEXE_headers += TransformParticle.H +CEXE_headers += DefaultInitialization.H +CEXE_headers += FilterCopyTransform.H +CEXE_headers += SmartCopy.H +CEXE_sources += SmartCopy.cpp INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Particles/ParticleCreation/ VPATH_LOCATIONS += $(WARPX_HOME)/Source/Particles/ParticleCreation/ diff --git a/Source/Particles/ParticleCreation/SmartCopy.H b/Source/Particles/ParticleCreation/SmartCopy.H new file mode 100644 index 000000000..e72a2f299 --- /dev/null +++ b/Source/Particles/ParticleCreation/SmartCopy.H @@ -0,0 +1,146 @@ +#ifndef SMART_COPY_H_ +#define SMART_COPY_H_ + +#include <DefaultInitialization.H> + +#include <AMReX_GpuContainers.H> +#include <AMReX_ParallelDescriptor.H> + +#include <map> +#include <string> + +using NameMap = std::map<std::string, int>; +using PolicyVec = amrex::Gpu::DeviceVector<InitializationPolicy>; + +struct SmartCopyTag +{ + std::vector<std::string> common_names; + amrex::Gpu::DeviceVector<int> src_comps; + amrex::Gpu::DeviceVector<int> dst_comps; + + int size () const noexcept { return common_names.size(); } +}; + +PolicyVec getPolicies (const NameMap& names) noexcept; + +SmartCopyTag getSmartCopyTag (const NameMap& src, const NameMap& dst) noexcept; + +template <typename PTile> +void setNewParticleIDs (PTile& ptile, int old_size, int num_added) +{ + int pid; +#ifdef _OPENMP +#pragma omp critical (ionization_nextid) +#endif + { + pid = PTile::ParticleType::NextID(); + PTile::ParticleType::NextID(pid + num_added); + } + + const int cpuid = amrex::ParallelDescriptor::MyProc(); + auto pp = ptile.GetArrayOfStructs()().data() + old_size; + amrex::For(num_added, [=] AMREX_GPU_DEVICE (int ip) noexcept + { + auto& p = pp[ip]; + p.id() = pid+ip; + p.cpu() = cpuid; + }); +} + +struct SmartCopy +{ + int m_num_copy_real; + const int* m_src_comps_r; + const int* m_dst_comps_r; + + int m_num_copy_int; + const int* m_src_comps_i; + const int* m_dst_comps_i; + + const InitializationPolicy* m_policy_real; + const InitializationPolicy* m_policy_int; + + template <typename DstData, typename SrcData> + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + void operator() (DstData& dst, const SrcData& src, int i_src, int i_dst) const noexcept + { + // the particle struct is always copied over + dst.m_aos[i_dst] = src.m_aos[i_src]; + + // initialize the real components + for (int j = 0; j < DstData::NAR; ++j) + dst.m_rdata[j][i_dst] = initializeRealValue(m_policy_real[j]); + for (int j = 0; j < dst.m_num_runtime_real; ++j) + dst.m_rdata[j][i_dst] = initializeRealValue(m_policy_real[j+DstData::NAR]); + + // initialize the int components + for (int j = 0; j < DstData::NAI; ++j) + dst.m_idata[j][i_dst] = initializeIntValue(m_policy_int[j]); + for (int j = 0; j < dst.m_num_runtime_int; ++j) + dst.m_idata[j][i_dst] = initializeIntValue(m_policy_int[j+DstData::NAI]); + + // copy the shared real components + for (int j = 0; j < m_num_copy_real; ++j) + { + int src_comp = (m_src_comps_r[j] < SrcData::NAR) ? m_src_comps_r[j] : m_src_comps_r[j] - SrcData::NAR; + int dst_comp = (m_dst_comps_r[j] < DstData::NAR) ? m_dst_comps_r[j] : m_dst_comps_r[j] - DstData::NAR; + + amrex::ParticleReal* AMREX_RESTRICT dst_data = (m_dst_comps_r[j] < DstData::NAR) ? dst.m_rdata[dst_comp] : dst.m_runtime_rdata[dst_comp]; + const amrex::ParticleReal* AMREX_RESTRICT src_data = (m_src_comps_r[j] < SrcData::NAR) ? src.m_rdata[src_comp] : src.m_runtime_rdata[src_comp]; + + dst_data[i_dst] = src_data[i_src]; + } + + // copy the shared int components + for (int j = 0; j < m_num_copy_int; ++j) + { + int src_comp = (m_src_comps_i[j] < SrcData::NAI) ? m_src_comps_i[j] : m_src_comps_i[j] - SrcData::NAI; + int dst_comp = (m_dst_comps_i[j] < DstData::NAI) ? m_dst_comps_i[j] : m_dst_comps_i[j] - DstData::NAI; + + int* AMREX_RESTRICT dst_data = (m_dst_comps_i[j] < DstData::NAI) ? dst.m_idata[dst_comp] : dst.m_runtime_idata[dst_comp]; + const int* AMREX_RESTRICT src_data = (m_src_comps_i[j] < SrcData::NAI) ? src.m_idata[src_comp] : src.m_runtime_idata[src_comp]; + + dst_data[i_dst] = src_data[i_src]; + } + } +}; + +class SmartCopyFactory +{ + SmartCopyTag m_tag_real; + SmartCopyTag m_tag_int; + PolicyVec m_policy_real; + PolicyVec m_policy_int; + bool m_defined; + +public: + template <class SrcPC, class DstPC> + SmartCopyFactory (const SrcPC& src, const DstPC& dst) noexcept + : m_defined(false) + { + m_tag_real = getSmartCopyTag(src.getParticleComps(), dst.getParticleComps()); + m_tag_int = getSmartCopyTag(src.getParticleiComps(), dst.getParticleiComps()); + + m_policy_real = getPolicies(dst.getParticleComps()); + m_policy_int = getPolicies(dst.getParticleiComps()); + + m_defined = true; + } + + SmartCopy getSmartCopy () const noexcept + { + AMREX_ASSERT(m_defined); + return SmartCopy{m_tag_real.size(), + m_tag_real.src_comps.dataPtr(), + m_tag_real.dst_comps.dataPtr(), + m_tag_int.size(), + m_tag_int. src_comps.dataPtr(), + m_tag_int. dst_comps.dataPtr(), + m_policy_real.dataPtr(), + m_policy_int.dataPtr()}; + } + + bool isDefined () const noexcept { return m_defined; } +}; + +#endif diff --git a/Source/Particles/ParticleCreation/SmartCopy.cpp b/Source/Particles/ParticleCreation/SmartCopy.cpp new file mode 100644 index 000000000..5ceb625a8 --- /dev/null +++ b/Source/Particles/ParticleCreation/SmartCopy.cpp @@ -0,0 +1,45 @@ +#include "SmartCopy.H" + +PolicyVec getPolicies (const NameMap& names) noexcept +{ + PolicyVec policies; + policies.resize(names.size()); + for (const auto& kv : names) + { + policies[kv.second] = initialization_policies[kv.first]; + } + return policies; +} + +SmartCopyTag getSmartCopyTag (const NameMap& src, const NameMap& dst) noexcept +{ + SmartCopyTag tag; + + // we use the fact that maps are sorted + auto i_src = src.begin(); + auto i_dst = dst.begin(); + while ( (i_src != src.end()) and (i_dst != dst.end()) ) + { + if (i_src->first < i_dst->first) + { + // names are not the same and src is lower + ++i_src; + } + else if (i_src->first > i_dst->first) + { + // names are not the same and dst is lower + ++i_dst; + } + else + { + // name is in both... + tag.common_names.push_back(i_src->first); + tag.src_comps.push_back(i_src->second); + tag.dst_comps.push_back(i_dst->second); + ++i_src; + ++i_dst; + } + } + + return tag; +} diff --git a/Source/Particles/ParticleCreation/TransformParticle.H b/Source/Particles/ParticleCreation/TransformParticle.H deleted file mode 100644 index eb5820e32..000000000 --- a/Source/Particles/ParticleCreation/TransformParticle.H +++ /dev/null @@ -1,50 +0,0 @@ -/* Copyright 2019 Axel Huebl, Maxence Thevenet - * - * This file is part of WarpX. - * - * License: BSD-3-Clause-LBNL - */ -#ifndef TRANSFORMPARTICLE_H_ -#define TRANSFORMPARTICLE_H_ - -#include "WarpXParticleContainer.H" - -enum struct ElementaryProcessType { Ionization }; - -/** - * \brief small modifications on source particle - * \param i index of particle - * \param particle particle struct - * \param runtime_iattribs integer attribs - */ -template < ElementaryProcessType ProcessT > -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void transformSourceParticle( - int i, - WarpXParticleContainer::ParticleType& particle, - amrex::GpuArray<int*,1> runtime_iattribs){} - -/** - * \brief small modifications on product particle - * \param i index of particle - * \param v_particle pointer to vector of particles - */ -template < ElementaryProcessType ProcessT > -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void transformProductParticle( - int i, - WarpXParticleContainer::ParticleType* * v_particle){} - -// For ionization, increase ionization level of source particle -template <> -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void transformSourceParticle < ElementaryProcessType::Ionization > ( - int i, - WarpXParticleContainer::ParticleType& particle, - amrex::GpuArray<int*,1> runtime_iattribs) -{ - // increment particle's ionization level - runtime_iattribs[0][i] += 1; -} - -#endif // TRANSFORMPARTICLE_H_ diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index 942f7b7a0..421dcd842 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -13,6 +13,7 @@ #include <PlasmaInjector.H> #include <WarpXParticleContainer.H> #include <NCIGodfreyFilter.H> +#include <Ionization.H> #include <AMReX_IArrayBox.H> @@ -187,8 +188,7 @@ public: void SplitParticles(int lev); - virtual void buildIonizationMask (const amrex::MFIter& mfi, const int lev, - amrex::Gpu::ManagedDeviceVector<int>& ionization_mask) override; + IonizationFilterFunc getIonizationFunc (); // Inject particles in Box 'part_box' virtual void AddParticles (int lev); diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 588f8f0ae..6572657ff 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -375,19 +375,7 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) scale_fac = dx[0]*dx[1]/num_ppc; #endif -#ifdef _OPENMP - // First touch all tiles in the map in serial - for (MFIter mfi = MakeMFIter(lev); mfi.isValid(); ++mfi) { - auto index = std::make_pair(mfi.index(), mfi.LocalTileIndex()); - GetParticles(lev)[index]; - tmp_particle_data.resize(finestLevel()+1); - // Create map entry if not there - tmp_particle_data[lev][index]; - if ( (NumRuntimeRealComps()>0) || (NumRuntimeIntComps()>0) ) { - DefineAndReturnParticleTile(lev, mfi.index(), mfi.LocalTileIndex()); - } - } -#endif + defineAllParticleTiles(); MultiFab* cost = WarpX::getCosts(lev); @@ -2338,86 +2326,17 @@ void PhysicalParticleContainer::InitIonizationModule () } } -/* \brief create mask of ionized particles (1 if ionized, 0 otherwise) - * - * \param mfi: tile or grid - * \param lev: MR level - * \param ionization_mask: Array with as many elements as particles in mfi. - * This function initialized the array, and set each element to 1 or 0 - * depending on whether the particle is ionized or not. - */ -void -PhysicalParticleContainer::buildIonizationMask (const amrex::MFIter& mfi, const int lev, - amrex::Gpu::ManagedDeviceVector<int>& ionization_mask) +IonizationFilterFunc +PhysicalParticleContainer::getIonizationFunc () { - BL_PROFILE("PPC::buildIonizationMask"); - // Get pointers to ionization data from pc_source - const Real * const AMREX_RESTRICT p_ionization_energies = ionization_energies.dataPtr(); - const Real * const AMREX_RESTRICT p_adk_prefactor = adk_prefactor.dataPtr(); - const Real * const AMREX_RESTRICT p_adk_exp_prefactor = adk_exp_prefactor.dataPtr(); - const Real * const AMREX_RESTRICT p_adk_power = adk_power.dataPtr(); - - // Current tile info - const int grid_id = mfi.index(); - const int tile_id = mfi.LocalTileIndex(); - - // Get GPU-friendly arrays of particle data - auto& ptile = GetParticles(lev)[std::make_pair(grid_id,tile_id)]; - // Only need attribs (i.e., SoA data) - auto& soa = ptile.GetStructOfArrays(); - const int np = ptile.GetArrayOfStructs().size(); - - // If no particle, nothing to do. - if (np == 0) return; - // Otherwise, resize ionization_mask, and get poiters to attribs arrays. - ionization_mask.resize(np); - int * const AMREX_RESTRICT p_ionization_mask = ionization_mask.data(); - const ParticleReal * const AMREX_RESTRICT ux = soa.GetRealData(PIdx::ux).data(); - const ParticleReal * const AMREX_RESTRICT uy = soa.GetRealData(PIdx::uy).data(); - const ParticleReal * const AMREX_RESTRICT uz = soa.GetRealData(PIdx::uz).data(); - const ParticleReal * const AMREX_RESTRICT ex = soa.GetRealData(PIdx::Ex).data(); - const ParticleReal * const AMREX_RESTRICT ey = soa.GetRealData(PIdx::Ey).data(); - const ParticleReal * const AMREX_RESTRICT ez = soa.GetRealData(PIdx::Ez).data(); - const ParticleReal * const AMREX_RESTRICT bx = soa.GetRealData(PIdx::Bx).data(); - const ParticleReal * const AMREX_RESTRICT by = soa.GetRealData(PIdx::By).data(); - const ParticleReal * const AMREX_RESTRICT bz = soa.GetRealData(PIdx::Bz).data(); - int* ion_lev = soa.GetIntData(particle_icomps["ionization_level"]).data(); - - Real c = PhysConst::c; - Real c2_inv = 1./c/c; - int atomic_number = ion_atomic_number; - - // Loop over all particles in grid/tile. If ionized, set mask to 1 - // and increment ionization level. - ParallelFor( - np, - [=] AMREX_GPU_DEVICE (long i) { - // Get index of ionization_level - p_ionization_mask[i] = 0; - if ( ion_lev[i]<atomic_number ){ - Real random_draw = amrex::Random(); - // Compute electric field amplitude in the particle's frame of - // reference (particularly important when in boosted frame). - Real ga = std::sqrt(1. + (ux[i]*ux[i] + uy[i]*uy[i] + uz[i]*uz[i]) * c2_inv); - Real E = std::sqrt( - - ( ux[i]*ex[i] + uy[i]*ey[i] + uz[i]*ez[i] ) * ( ux[i]*ex[i] + uy[i]*ey[i] + uz[i]*ez[i] ) * c2_inv - + ( ga *ex[i] + uy[i]*bz[i] - uz[i]*by[i] ) * ( ga *ex[i] + uy[i]*bz[i] - uz[i]*by[i] ) - + ( ga *ey[i] + uz[i]*bx[i] - ux[i]*bz[i] ) * ( ga *ey[i] + uz[i]*bx[i] - ux[i]*bz[i] ) - + ( ga *ez[i] + ux[i]*by[i] - uy[i]*bx[i] ) * ( ga *ez[i] + ux[i]*by[i] - uy[i]*bx[i] ) - ); - // Compute probability of ionization p - Real w_dtau = 1./ ga * p_adk_prefactor[ion_lev[i]] * - std::pow(E,p_adk_power[ion_lev[i]]) * - std::exp( p_adk_exp_prefactor[ion_lev[i]]/E ); - Real p = 1. - std::exp( - w_dtau ); - - if (random_draw < p){ - // update mask - p_ionization_mask[i] = 1; - } - } - } - ); + BL_PROFILE("PPC::getIonizationFunc"); + + return IonizationFilterFunc{ionization_energies.dataPtr(), + adk_prefactor.dataPtr(), + adk_exp_prefactor.dataPtr(), + adk_power.dataPtr(), + particle_icomps["ionization_level"], + ion_atomic_number}; } //This function return true if the PhysicalParticleContainer contains electrons diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index 845977b00..b1def39d3 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -302,12 +302,8 @@ public: int doBackTransformedDiagnostics () const { return do_back_transformed_diagnostics; } - virtual void buildIonizationMask (const amrex::MFIter& mfi, const int lev, - amrex::Gpu::ManagedDeviceVector<int>& ionization_mask) - {}; - - std::map<std::string, int> getParticleComps () { return particle_comps;} - std::map<std::string, int> getParticleiComps () { return particle_icomps;} + std::map<std::string, int> getParticleComps () const noexcept { return particle_comps;} + std::map<std::string, int> getParticleiComps () const noexcept { return particle_icomps;} //amrex::Real getCharge () {return charge;} amrex::ParticleReal getCharge () const {return charge;} @@ -386,6 +382,8 @@ protected: amrex::Vector<std::map<PairIndex, std::array<DataContainer, TmpIdx::nattribs> > > tmp_particle_data; + void defineAllParticleTiles () noexcept; + private: virtual void particlePostLocate(ParticleType& p, const amrex::ParticleLocData& pld, const int lev) override; diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 59699d710..cd5bcaaf9 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -748,6 +748,25 @@ WarpXParticleContainer::PushX (int lev, amrex::Real dt) } } +// When using runtime components, AMReX requires to touch all tiles +// in serial and create particles tiles with runtime components if +// they do not exist (or if they were defined by default, i.e., +// without runtime component). +void WarpXParticleContainer::defineAllParticleTiles () noexcept +{ + tmp_particle_data.resize(finestLevel()+1); + for (int lev = 0; lev <= finestLevel(); ++lev) + { + for (auto mfi = MakeMFIter(lev); mfi.isValid(); ++mfi) + { + const int grid_id = mfi.index(); + const int tile_id = mfi.LocalTileIndex(); + tmp_particle_data[lev][std::make_pair(grid_id,tile_id)]; + DefineAndReturnParticleTile(lev, grid_id, tile_id); + } + } +} + // This function is called in Redistribute, just after locate void WarpXParticleContainer::particlePostLocate(ParticleType& p, |