diff options
Diffstat (limited to 'Source')
21 files changed, 1123 insertions, 588 deletions
diff --git a/Source/Diagnostics/ReducedDiags/BeamRelevant.H b/Source/Diagnostics/ReducedDiags/BeamRelevant.H new file mode 100644 index 000000000..fa733ae66 --- /dev/null +++ b/Source/Diagnostics/ReducedDiags/BeamRelevant.H @@ -0,0 +1,35 @@ +/* Copyright 2019-2020 Yinjian Zhao + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ + +#ifndef WARPX_DIAGNOSTICS_REDUCEDDIAGS_BEAMRELEVANT_H_ +#define WARPX_DIAGNOSTICS_REDUCEDDIAGS_BEAMRELEVANT_H_ + +#include "ReducedDiags.H" +#include <fstream> + +/** + * This class contains diagnostics that are relevant to beam. + */ +class BeamRelevant : public ReducedDiags +{ +public: + + /** constructor + * @param[in] rd_name reduced diags names */ + BeamRelevant(std::string rd_name); + + /// name of beam species + std::string m_beam_name; + + /** This funciton computes beam relevant quantites. + * \param [in] step current time step + */ + virtual void ComputeDiags(int step) override final; + +}; + +#endif diff --git a/Source/Diagnostics/ReducedDiags/BeamRelevant.cpp b/Source/Diagnostics/ReducedDiags/BeamRelevant.cpp new file mode 100644 index 000000000..b82d24645 --- /dev/null +++ b/Source/Diagnostics/ReducedDiags/BeamRelevant.cpp @@ -0,0 +1,389 @@ +/* Copyright 2019-2020 Yinjian Zhao + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ + +#include "BeamRelevant.H" +#include "WarpX.H" +#include "WarpXConst.H" +#include "AMReX_REAL.H" +#include "AMReX_ParticleReduce.H" +#include <iostream> +#include <cmath> +#include <limits> + +using namespace amrex; + +// constructor +BeamRelevant::BeamRelevant (std::string rd_name) +: ReducedDiags{rd_name} +{ + +#if (defined WARPX_DIM_RZ) + // RZ coordinate is not working + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(false, + "BeamRelevant reduced diagnostics does not work for RZ coordinate."); +#endif + + // read beam name + ParmParse pp(rd_name); + pp.get("species",m_beam_name); + + // resize data array +#if (AMREX_SPACEDIM == 3) + // 0, 1, 2: mean x,y,z + // 3, 4, 5: mean px,py,pz + // 6: gamma + // 7, 8, 9: rms x,y,z + // 10,11,12: rms px,py,pz + // 13: rms gamma + // 14,15,16: emittance x,y,z + m_data.resize(17,0.0); +#elif (AMREX_SPACEDIM == 2) + // 0, 1: mean x,z + // 2, 3, 4: mean px,py,pz + // 5: gamma + // 6, 7: rms x,z + // 8, 9,10: rms px,py,pz + // 11: rms gamma + // 12,13: emittance x,z + m_data.resize(14,0.0); +#endif + + if (ParallelDescriptor::IOProcessor()) + { + if ( m_IsNotRestart ) + { + // open file + std::ofstream ofs; + ofs.open(m_path + m_rd_name + "." + m_extension, + std::ofstream::out | std::ofstream::app); + // write header row +#if (AMREX_SPACEDIM == 3) + ofs << "#"; + ofs << "[1]step"; ofs << m_sep; + ofs << "[2]time(s)"; ofs << m_sep; + ofs << "[3]x_mean(m)"; ofs << m_sep; + ofs << "[4]y_mean(m)"; ofs << m_sep; + ofs << "[5]z_mean(m)"; ofs << m_sep; + ofs << "[6]px_mean(kg m/s)"; ofs << m_sep; + ofs << "[7]py_mean(kg m/s)"; ofs << m_sep; + ofs << "[8]pz_mean(kg m/s)"; ofs << m_sep; + ofs << "[9]gamma_mean"; ofs << m_sep; + ofs << "[10]x_rms(m)"; ofs << m_sep; + ofs << "[11]y_rms(m)"; ofs << m_sep; + ofs << "[12]z_rms(m)"; ofs << m_sep; + ofs << "[13]px_rms(kg m/s)"; ofs << m_sep; + ofs << "[14]py_rms(kg m/s)"; ofs << m_sep; + ofs << "[15]pz_rms(kg m/s)"; ofs << m_sep; + ofs << "[16]gamma_rms"; ofs << m_sep; + ofs << "[17]emittance_x(m)"; ofs << m_sep; + ofs << "[18]emittance_y(m)"; ofs << m_sep; + ofs << "[19]emittance_z(m)"; ofs << std::endl; +#elif (AMREX_SPACEDIM == 2) + ofs << "#"; + ofs << "[1]step"; ofs << m_sep; + ofs << "[2]time(s)"; ofs << m_sep; + ofs << "[3]x_mean(m)"; ofs << m_sep; + ofs << "[4]z_mean(m)"; ofs << m_sep; + ofs << "[5]px_mean(kg m/s)"; ofs << m_sep; + ofs << "[6]py_mean(kg m/s)"; ofs << m_sep; + ofs << "[7]pz_mean(kg m/s)"; ofs << m_sep; + ofs << "[8]gamma_mean"; ofs << m_sep; + ofs << "[9]x_rms(m)"; ofs << m_sep; + ofs << "[10]z_rms(m)"; ofs << m_sep; + ofs << "[11]px_rms(kg m/s)"; ofs << m_sep; + ofs << "[12]py_rms(kg m/s)"; ofs << m_sep; + ofs << "[13]pz_rms(kg m/s)"; ofs << m_sep; + ofs << "[14]gamma_rms"; ofs << m_sep; + ofs << "[15]emittance_x(m)"; ofs << m_sep; + ofs << "[16]emittance_z(m)"; ofs << std::endl; +#endif + // close file + ofs.close(); + } + } + +} +// end constructor + +// function that compute beam relevant quantities +void BeamRelevant::ComputeDiags (int step) +{ + + // Judge if the diags should be done + if ( (step+1) % m_freq != 0 ) { return; } + + // get MultiParticleContainer class object + auto & mypc = WarpX::GetInstance().GetPartContainer(); + + // get number of species + int const nSpecies = mypc.nSpecies(); + + // get species names (std::vector<std::string>) + auto const species_names = mypc.GetSpeciesNames(); + + // inverse of speed of light squared + Real constexpr inv_c2 = 1.0 / (PhysConst::c * PhysConst::c); + + // If 2D-XZ, p.pos(1) is z, rather than p.pos(2). +#if (AMREX_SPACEDIM == 3) + int const index_z = 2; +#elif (AMREX_SPACEDIM == 2) + int const index_z = 1; +#endif + + // loop over species + for (int i_s = 0; i_s < nSpecies; ++i_s) + { + + // only beam species does + if (species_names[i_s] != m_beam_name) { continue; } + + // get WarpXParticleContainer class object + auto const & myspc = mypc.GetParticleContainer(i_s); + + // get mass (Real) + Real const m = myspc.getMass(); + + using PType = typename WarpXParticleContainer::SuperParticleType; + + // weight sum + Real w_sum = ReduceSum( myspc, + [=] AMREX_GPU_HOST_DEVICE (const PType& p) -> Real + { return p.rdata(PIdx::w); }); + + // reduced sum over mpi ranks + ParallelDescriptor::ReduceRealSum(w_sum); + + if (w_sum < std::numeric_limits<Real>::min() ) + { + for (int i = 0; i < m_data.size(); ++i) { m_data[i] = 0.0; } + return; + } + + // x mean + Real x_mean = ReduceSum( myspc, + [=] AMREX_GPU_HOST_DEVICE (const PType& p) -> Real + { return p.pos(0) * p.rdata(PIdx::w); }); + +#if (AMREX_SPACEDIM == 3) + // y mean + Real y_mean = ReduceSum( myspc, + [=] AMREX_GPU_HOST_DEVICE (const PType& p) -> Real + { return p.pos(1) * p.rdata(PIdx::w); }); +#endif + + // z mean + Real z_mean = ReduceSum( myspc, + [=] AMREX_GPU_HOST_DEVICE (const PType& p) -> Real + { return p.pos(index_z) * p.rdata(PIdx::w); }); + + // ux mean + Real ux_mean = ReduceSum( myspc, + [=] AMREX_GPU_HOST_DEVICE (const PType& p) -> Real + { return p.rdata(PIdx::ux) * p.rdata(PIdx::w); }); + + // uy mean + Real uy_mean = ReduceSum( myspc, + [=] AMREX_GPU_HOST_DEVICE (const PType& p) -> Real + { return p.rdata(PIdx::uy) * p.rdata(PIdx::w); }); + + // uz mean + Real uz_mean = ReduceSum( myspc, + [=] AMREX_GPU_HOST_DEVICE (const PType& p) -> Real + { return p.rdata(PIdx::uz) * p.rdata(PIdx::w); }); + + // gamma mean + Real gm_mean = ReduceSum( myspc, + [=] AMREX_GPU_HOST_DEVICE (const PType& p) -> Real + { + Real ux = p.rdata(PIdx::ux); + Real uy = p.rdata(PIdx::uy); + Real uz = p.rdata(PIdx::uz); + Real us = ux*ux + uy*uy + uz*uz; + return std::sqrt(1.0 + us*inv_c2) * p.rdata(PIdx::w); + }); + + // reduced sum over mpi ranks + ParallelDescriptor::ReduceRealSum(x_mean); x_mean /= w_sum; +#if (AMREX_SPACEDIM == 3) + ParallelDescriptor::ReduceRealSum(y_mean); y_mean /= w_sum; +#endif + ParallelDescriptor::ReduceRealSum(z_mean); z_mean /= w_sum; + ParallelDescriptor::ReduceRealSum(ux_mean); ux_mean /= w_sum; + ParallelDescriptor::ReduceRealSum(uy_mean); uy_mean /= w_sum; + ParallelDescriptor::ReduceRealSum(uz_mean); uz_mean /= w_sum; + ParallelDescriptor::ReduceRealSum(gm_mean); gm_mean /= w_sum; + + // x mean square + Real x_ms = ReduceSum( myspc, + [=] AMREX_GPU_HOST_DEVICE (const PType& p) -> Real + { + Real const a = (p.pos(0)-x_mean) * (p.pos(0)-x_mean); + return a * p.rdata(PIdx::w); + }); + +#if (AMREX_SPACEDIM == 3) + // y mean square + Real y_ms = ReduceSum( myspc, + [=] AMREX_GPU_HOST_DEVICE (const PType& p) -> Real + { + Real const a = (p.pos(1)-y_mean) * (p.pos(1)-y_mean); + return a * p.rdata(PIdx::w); + }); +#endif + + // z mean square + Real z_ms = ReduceSum( myspc, + [=] AMREX_GPU_HOST_DEVICE (const PType& p) -> Real + { + Real const a = (p.pos(index_z)-z_mean) * (p.pos(index_z)-z_mean); + return a * p.rdata(PIdx::w); + }); + + // ux mean square + Real ux_ms = ReduceSum( myspc, + [=] AMREX_GPU_HOST_DEVICE (const PType& p) -> Real + { + Real const a = (p.rdata(PIdx::ux)-ux_mean) * + (p.rdata(PIdx::ux)-ux_mean); + return a * p.rdata(PIdx::w); + }); + + // uy mean square + Real uy_ms = ReduceSum( myspc, + [=] AMREX_GPU_HOST_DEVICE (const PType& p) -> Real + { + Real const a = (p.rdata(PIdx::uy)-uy_mean) * + (p.rdata(PIdx::uy)-uy_mean); + return a * p.rdata(PIdx::w); + }); + + // uz mean square + Real uz_ms = ReduceSum( myspc, + [=] AMREX_GPU_HOST_DEVICE (const PType& p) -> Real + { + Real const a = (p.rdata(PIdx::uz)-uz_mean) * + (p.rdata(PIdx::uz)-uz_mean); + return a * p.rdata(PIdx::w); + }); + + // gamma mean square + Real gm_ms = ReduceSum( myspc, + [=] AMREX_GPU_HOST_DEVICE (const PType& p) -> Real + { + Real const ux = p.rdata(PIdx::ux); + Real const uy = p.rdata(PIdx::uy); + Real const uz = p.rdata(PIdx::uz); + Real const us = ux*ux + uy*uy + uz*uz; + Real const gm = std::sqrt(1.0 + us*inv_c2); + Real const a = (gm - gm_mean) * (gm - gm_mean); + return a * p.rdata(PIdx::w); + }); + + // x times ux + Real xux = ReduceSum( myspc, + [=] AMREX_GPU_HOST_DEVICE (const PType& p) -> Real + { + Real const a = (p.pos(0)-x_mean) * (p.rdata(PIdx::ux)-ux_mean); + return a * p.rdata(PIdx::w); + }); + +#if (AMREX_SPACEDIM == 3) + // y times uy + Real yuy = ReduceSum( myspc, + [=] AMREX_GPU_HOST_DEVICE (const PType& p) -> Real + { + Real const a = (p.pos(1)-y_mean) * (p.rdata(PIdx::uy)-uy_mean); + return a * p.rdata(PIdx::w); + }); +#endif + + // z times uz + Real zuz = ReduceSum( myspc, + [=] AMREX_GPU_HOST_DEVICE (const PType& p) -> Real + { + Real const a = (p.pos(index_z)-z_mean) * (p.rdata(PIdx::uz)-uz_mean); + return a * p.rdata(PIdx::w); + }); + + // reduced sum over mpi ranks + ParallelDescriptor::ReduceRealSum + ( x_ms, ParallelDescriptor::IOProcessorNumber()); + x_ms /= w_sum; +#if (AMREX_SPACEDIM == 3) + ParallelDescriptor::ReduceRealSum + ( y_ms, ParallelDescriptor::IOProcessorNumber()); + y_ms /= w_sum; +#endif + ParallelDescriptor::ReduceRealSum + ( z_ms, ParallelDescriptor::IOProcessorNumber()); + z_ms /= w_sum; + ParallelDescriptor::ReduceRealSum + (ux_ms, ParallelDescriptor::IOProcessorNumber()); + ux_ms /= w_sum; + ParallelDescriptor::ReduceRealSum + (uy_ms, ParallelDescriptor::IOProcessorNumber()); + uy_ms /= w_sum; + ParallelDescriptor::ReduceRealSum + (uz_ms, ParallelDescriptor::IOProcessorNumber()); + uz_ms /= w_sum; + ParallelDescriptor::ReduceRealSum + (gm_ms, ParallelDescriptor::IOProcessorNumber()); + gm_ms /= w_sum; + ParallelDescriptor::ReduceRealSum + ( xux, ParallelDescriptor::IOProcessorNumber()); + xux /= w_sum; +#if (AMREX_SPACEDIM == 3) + ParallelDescriptor::ReduceRealSum + ( yuy, ParallelDescriptor::IOProcessorNumber()); + yuy /= w_sum; +#endif + ParallelDescriptor::ReduceRealSum + ( zuz, ParallelDescriptor::IOProcessorNumber()); + zuz /= w_sum; + + // save data +#if (AMREX_SPACEDIM == 3) + m_data[0] = x_mean; + m_data[1] = y_mean; + m_data[2] = z_mean; + m_data[3] = ux_mean * m; + m_data[4] = uy_mean * m; + m_data[5] = uz_mean * m; + m_data[6] = gm_mean; + m_data[7] = std::sqrt(x_ms); + m_data[8] = std::sqrt(y_ms); + m_data[9] = std::sqrt(z_ms); + m_data[10] = std::sqrt(ux_ms * m); + m_data[11] = std::sqrt(uy_ms * m); + m_data[12] = std::sqrt(uz_ms * m); + m_data[13] = std::sqrt(gm_ms); + m_data[14] = std::sqrt(std::abs(x_ms*ux_ms-xux*xux)) / PhysConst::c; + m_data[15] = std::sqrt(std::abs(y_ms*uy_ms-yuy*yuy)) / PhysConst::c; + m_data[16] = std::sqrt(std::abs(z_ms*uz_ms-zuz*zuz)) / PhysConst::c; +#elif (AMREX_SPACEDIM == 2) + m_data[0] = x_mean; + m_data[1] = z_mean; + m_data[2] = ux_mean * m; + m_data[3] = uy_mean * m; + m_data[4] = uz_mean * m; + m_data[5] = gm_mean; + m_data[6] = std::sqrt(x_ms); + m_data[7] = std::sqrt(z_ms); + m_data[8] = std::sqrt(ux_ms * m); + m_data[9] = std::sqrt(uy_ms * m); + m_data[10] = std::sqrt(uz_ms * m); + m_data[11] = std::sqrt(gm_ms); + m_data[12] = std::sqrt(std::abs(x_ms*ux_ms-xux*xux)) / PhysConst::c; + m_data[13] = std::sqrt(std::abs(z_ms*uz_ms-zuz*zuz)) / PhysConst::c; +#endif + + } + // end loop over species + +} +// end void BeamRelevant::ComputeDiags diff --git a/Source/Diagnostics/ReducedDiags/Make.package b/Source/Diagnostics/ReducedDiags/Make.package index 37f76d3d5..4b1207593 100644 --- a/Source/Diagnostics/ReducedDiags/Make.package +++ b/Source/Diagnostics/ReducedDiags/Make.package @@ -10,5 +10,8 @@ CEXE_sources += ParticleEnergy.cpp CEXE_headers += FieldEnergy.H CEXE_sources += FieldEnergy.cpp +CEXE_headers += BeamRelevant.H +CEXE_sources += BeamRelevant.cpp + INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Diagnostics/ReducedDiags VPATH_LOCATIONS += $(WARPX_HOME)/Source/Diagnostics/ReducedDiags diff --git a/Source/Diagnostics/ReducedDiags/MultiReducedDiags.cpp b/Source/Diagnostics/ReducedDiags/MultiReducedDiags.cpp index 71a33924b..2544a90d9 100644 --- a/Source/Diagnostics/ReducedDiags/MultiReducedDiags.cpp +++ b/Source/Diagnostics/ReducedDiags/MultiReducedDiags.cpp @@ -5,6 +5,7 @@ * License: BSD-3-Clause-LBNL */ +#include "BeamRelevant.H" #include "ParticleEnergy.H" #include "FieldEnergy.H" #include "MultiReducedDiags.H" @@ -49,6 +50,11 @@ MultiReducedDiags::MultiReducedDiags () m_multi_rd[i_rd].reset ( new FieldEnergy(m_rd_names[i_rd])); } + else if (rd_type.compare("BeamRelevant") == 0) + { + m_multi_rd[i_rd].reset + ( new BeamRelevant(m_rd_names[i_rd])); + } else { Abort("No matching reduced diagnostics type found."); } // end if match diags diff --git a/Source/Particles/ElementaryProcess/Ionization.H b/Source/Particles/ElementaryProcess/Ionization.H new file mode 100644 index 000000000..279339d07 --- /dev/null +++ b/Source/Particles/ElementaryProcess/Ionization.H @@ -0,0 +1,80 @@ +/* Copyright 2019-2020 Andrew Myers, Axel Huebl, + * Maxence Thevenet + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ +#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 + { + const int ion_lev = ptd.m_runtime_idata[comp][i]; + if (ion_lev < m_atomic_number) + { + constexpr amrex::Real c = PhysConst::c; + constexpr 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..8d951263e 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..3744fd9b3 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); @@ -489,11 +488,11 @@ MultiParticleContainer::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<WarpXParticleContainer::DiagnosticParticleData>& parts) const +::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<WarpXParticleContainer::DiagnosticParticleData>& parts) const { BL_PROFILE("MultiParticleContainer::GetLabFrameData"); @@ -555,7 +554,7 @@ MultiParticleContainer * calls virtual function ContinuousInjection. */ void -MultiParticleContainer::ContinuousInjection(const RealBox& injection_box) const +MultiParticleContainer::ContinuousInjection (const RealBox& injection_box) const { for (int i=0; i<nspecies+nlasers; i++){ auto& pc = allcontainers[i]; @@ -571,7 +570,7 @@ MultiParticleContainer::ContinuousInjection(const RealBox& injection_box) const * a position to update (PhysicalParticleContainer does not do anything). */ void -MultiParticleContainer::UpdateContinuousInjectionPosition(Real dt) const +MultiParticleContainer::UpdateContinuousInjectionPosition (Real dt) const { for (int i=0; i<nspecies+nlasers; i++){ auto& pc = allcontainers[i]; @@ -642,77 +641,45 @@ 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 +#pragma omp parallel if (Gpu::notInLaunchRegion()) #endif - // Loop over all grids (if not tiling) or grids and tiles (if tiling) 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_dst = dst_tile.numParticles(); + auto num_added = filterCopyTransformParticles(dst_tile, src_tile, np_dst, + Filter, Copy, Transform); + + setNewParticleIDs(dst_tile, np_dst, num_added); } - } // lev - } // pc_source + } + } } void @@ -735,7 +702,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 +716,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..2d9832c30 --- /dev/null +++ b/Source/Particles/ParticleCreation/DefaultInitialization.H @@ -0,0 +1,83 @@ +/* Copyright 2019-2020 Andrew Myers, Axel Huebl, + * Maxence Thevenet + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ +#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 + * 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, 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::RandomTau : { + return amrex::Random(); + } + 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; + 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..6dcb2a1e5 --- /dev/null +++ b/Source/Particles/ParticleCreation/FilterCopyTransform.H @@ -0,0 +1,141 @@ +/* Copyright 2019-2020 Andrew Myers, Axel Huebl, + * Maxence Thevenet + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ +#ifndef FILTER_COPY_TRANSFORM_H_ +#define FILTER_COPY_TRANSFORM_H_ + +#include <AMReX_GpuContainers.H> +#include <AMReX_TypeTraits.H> + +/** + * \brief Apply a filter, copy, and transform operation to the particles + * in src, in that order, writing the result to dst, starting at dst_index. + * The dst tile will be extended so all the particles will fit, if needed. + * + * Note that the transform function operates on both the src and the dst, + * so both can be modified. + * + * \tparam DstTile the dst particle tile type + * \tparam SrcTile the src particle tile type + * \tparam Index the index type, e.g. unsigned int + * \tparam TransFunc the transform function type + * \tparam CopyFunc the copy function type + * + * \param dst the destination tile + * \param src the source tile + * \param mask pointer to the mask - 1 means copy, 0 means don't copy + * \param dst_index the location at which to starting writing the result to dst + * \param copy callable that defines what will be done for the "copy" step. + * \param transform callable that defines the transformation to apply on dst and src. + * + * The form of the callable should model: + * template <typename DstData, typename SrcData> + * AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + * void operator() (DstData& dst, SrcData& src, int i_src, int i_dst); + * + * where dst and src refer to the destination and source tiles and + * i_src and i_dst and the particle indices in each tile. + * + * \return num_added the number of particles that were written to dst. + */ +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; + + const auto np = src.numParticles(); + if (np == 0) return 0; + + 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); + + const Index num_added = last_mask + last_offset; + dst.resize(std::max(dst_index + num_added, dst.numParticles())); + + const 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 num_added; +} + +/** + * \brief Apply a filter, copy, and transform operation to the particles + * in src, in that order, writing the result to dst, starting at dst_index. + * The dst tile will be extended so all the particles will fit, if needed. + * + * Note that the transform function operates on both the src and the dst, + * so both can be modified. + * + * \tparam DstTile the dst particle tile type + * \tparam SrcTile the src particle tile type + * \tparam Index the index type, e.g. unsigned int + * \tparam Filter the filter function type + * \tparam TransFunc the transform function type + * \tparam CopyFunc the copy function type + * + * \param dst the destination tile + * \param src the source tile + * \param dst_index the location at which to starting writing the result to dst + * \param filter a callable returning true if that particle is to be copied and transformed + * \param copy callable that defines what will be done for the "copy" step. + * \param transform callable that defines the transformation to apply on dst and src. + * + * The form of the callable should model: + * template <typename DstData, typename SrcData> + * AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + * void operator() (DstData& dst, SrcData& src, int i_src, int i_dst); + * + * where dst and src refer to the destination and source tiles and + * i_src and i_dst and the particle indices in each tile. + * + * \return num_added the number of particles that were written to dst. + */ +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; + + const auto np = src.numParticles(); + if (np == 0) return 0; + + 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..0668321c9 --- /dev/null +++ b/Source/Particles/ParticleCreation/SmartCopy.H @@ -0,0 +1,226 @@ +/* Copyright 2019-2020 Andrew Myers, Axel Huebl, + * Maxence Thevenet + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ +#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; + +/** + * \brief Sets the ids of newly created particles to the next values. + * + * \tparam PTile the particle tile type + * + * \param ptile the particle tile + * \param old_size the index of the first new particle + * \param num_added the number of particles to set the ids for. + */ +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::ParallelFor(num_added, [=] AMREX_GPU_DEVICE (int ip) noexcept + { + auto& p = pp[ip]; + p.id() = pid+ip; + p.cpu() = cpuid; + }); +} + +/** + * \brief This is a functor for performing a "smart copy" that works + * in both host and device code. + * + * A "smart" copy does the following. First, the destination particle + * components are initialized to the default values for that component + * type. Second, if a given component name is found in both the src + * and the dst, then the src value is copied. + * + * Particle structs - positions and id numbers - are always copied. + * + * You don't create this directly - use the SmartCopyFactory object below. + */ +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_runtime_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_runtime_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, dst_comp; + amrex::ParticleReal* AMREX_RESTRICT dst_data; + const amrex::ParticleReal* AMREX_RESTRICT src_data; + + if (m_src_comps_r[j] < SrcData::NAR) + { + // This is a compile-time attribute of the src + src_comp = m_src_comps_r[j]; + src_data = src.m_rdata[src_comp]; + } + else + { + // This is a runtime attribute of the src + src_comp = m_src_comps_r[j] - SrcData::NAR; + src_data = src.m_runtime_rdata[src_comp]; + } + + if (m_dst_comps_r[j] < DstData::NAR) + { + // This is a compile-time attribute of the dst + dst_comp = m_dst_comps_r[j]; + dst_data = dst.m_rdata[dst_comp]; + } + else + { + // This is a runtime attribute of the dst + dst_comp = m_dst_comps_r[j] - DstData::NAR; + dst_data = dst.m_runtime_rdata[dst_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, dst_comp; + int* AMREX_RESTRICT dst_data; + int* AMREX_RESTRICT src_data; + + if (m_src_comps_i[j] < SrcData::NAI) + { + src_comp = m_src_comps_i[j]; + src_data = src.m_idata[src_comp]; + } + else + { + src_comp = m_src_comps_i[j] - SrcData::NAI; + src_data = src.m_runtime_idata[src_comp]; + } + + if (m_dst_comps_i[j] < DstData::NAI) + { + dst_comp = m_dst_comps_i[j]; + dst_data = dst.m_idata[dst_comp]; + } + else + { + dst_comp = m_dst_comps_i[j] - DstData::NAI; + dst_data = dst.m_runtime_idata[dst_comp]; + } + + dst_data[i_dst] = src_data[i_src]; + } + } +}; + +/** + * \brief A factory for creating SmartCopy functors. + * + * Given two particle containers, this can create a functor + * that will perform the smart copy operation between those + * particle container's tiles. + */ +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..7d323ce39 --- /dev/null +++ b/Source/Particles/ParticleCreation/SmartCopy.cpp @@ -0,0 +1,52 @@ +/* Copyright 2019-2020 Andrew Myers, Axel Huebl, + * Maxence Thevenet + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ +#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..0080706f7 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,14 @@ protected: amrex::Vector<std::map<PairIndex, std::array<DataContainer, TmpIdx::nattribs> > > tmp_particle_data; + /** + * 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 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, |