diff options
Diffstat (limited to 'Source')
| -rw-r--r-- | Source/FortranInterface/WarpX_f.H | 1 | ||||
| -rw-r--r-- | Source/Make.WarpX | 19 | ||||
| -rw-r--r-- | Source/Particles/Make.package | 6 | ||||
| -rw-r--r-- | Source/Particles/MultiParticleContainer.H | 41 | ||||
| -rw-r--r-- | Source/Particles/MultiParticleContainer.cpp | 82 | ||||
| -rw-r--r-- | Source/Particles/PhotonParticleContainer.H | 83 | ||||
| -rw-r--r-- | Source/Particles/PhotonParticleContainer.cpp | 115 | ||||
| -rw-r--r-- | Source/Particles/WarpXParticleContainer.cpp | 14 | ||||
| -rw-r--r-- | Source/QED/Make.package | 6 | ||||
| -rw-r--r-- | Source/QED/amrex_rng_wrapper.cpp | 18 | ||||
| -rw-r--r-- | Source/QED/amrex_rng_wrapper.h | 20 | ||||
| -rw-r--r-- | Source/QED/breit_wheeler_engine_wrapper.h | 24 |
12 files changed, 395 insertions, 34 deletions
diff --git a/Source/FortranInterface/WarpX_f.H b/Source/FortranInterface/WarpX_f.H index 60a9817a2..26da42606 100644 --- a/Source/FortranInterface/WarpX_f.H +++ b/Source/FortranInterface/WarpX_f.H @@ -63,6 +63,7 @@ extern "C" { #endif + // Laser pusher void warpx_gaussian_laser( const long* np, amrex::Real* Xp, amrex::Real* Yp, amrex::Real* t, amrex::Real* wavelength, amrex::Real* e_max, amrex::Real* waist, diff --git a/Source/Make.WarpX b/Source/Make.WarpX index ea06f8e06..d369e38e2 100644 --- a/Source/Make.WarpX +++ b/Source/Make.WarpX @@ -81,6 +81,21 @@ ifeq ($(USE_SENSEI_INSITU),TRUE) include $(AMREX_HOME)/Src/Extern/SENSEI/Make.package endif +ifeq ($(QED),TRUE) + BOOST_ROOT ?= NOT_SET + ifneq ($(BOOST_ROOT),NOT_SET) + VPATH_LOCATIONS += $(BOOST_ROOT) + INCLUDE_LOCATIONS += $(BOOST_ROOT) + endif + INCLUDE_LOCATIONS += $(PICSAR_HOME)/src/multi_physics/QED/src + CXXFLAGS += -DWARPX_QED + CFLAGS += -DWARPX_QED + FFLAGS += -DWARPX_QED + F90FLAGS += -DWARPX_QED + include $(WARPX_HOME)/Source/QED/Make.package +endif + + include $(PICSAR_HOME)/src/Make.package WARPX_GIT_VERSION := $(shell cd $(WARPX_HOME); git describe --abbrev=12 --dirty --always --tags) @@ -160,7 +175,11 @@ ifeq ($(USE_HDF5),TRUE) LIBRARY_LOCATIONS += $(HDF5_HOME)/lib endif DEFINES += -DWARPX_USE_HDF5 +<<<<<<< HEAD + LIBRARIES += -lhdf5 -lz +======= libraries += -lhdf5 -lz +>>>>>>> dev endif # job_info support diff --git a/Source/Particles/Make.package b/Source/Particles/Make.package index db90de1dc..95f36cf65 100644 --- a/Source/Particles/Make.package +++ b/Source/Particles/Make.package @@ -1,14 +1,18 @@ F90EXE_sources += deposit_cic.F90 F90EXE_sources += interpolate_cic.F90 F90EXE_sources += push_particles_ES.F90 -CEXE_sources += PhysicalParticleContainer.cpp + CEXE_sources += MultiParticleContainer.cpp CEXE_sources += WarpXParticleContainer.cpp CEXE_sources += RigidInjectedParticleContainer.cpp +CEXE_sources += PhysicalParticleContainer.cpp +CEXE_sources += PhotonParticleContainer.cpp + CEXE_headers += MultiParticleContainer.H CEXE_headers += WarpXParticleContainer.H CEXE_headers += RigidInjectedParticleContainer.H CEXE_headers += PhysicalParticleContainer.H +CEXE_headers += PhotonParticleContainer.H CEXE_headers += ShapeFactors.H include $(WARPX_HOME)/Source/Particles/Pusher/Make.package diff --git a/Source/Particles/MultiParticleContainer.H b/Source/Particles/MultiParticleContainer.H index 1aff7edfb..67392c449 100644 --- a/Source/Particles/MultiParticleContainer.H +++ b/Source/Particles/MultiParticleContainer.H @@ -10,6 +10,7 @@ #include <WarpXParticleContainer.H> #include <PhysicalParticleContainer.H> #include <RigidInjectedParticleContainer.H> +#include <PhotonParticleContainer.H> #include <LaserParticleContainer.H> // @@ -20,6 +21,11 @@ // concrete class dervied from WarpXParticlContainer. // +// +#ifdef WARPX_QED + #include "breit_wheeler_engine_wrapper.h" +#endif + class MultiParticleContainer { @@ -50,24 +56,24 @@ public: const amrex::Vector<std::unique_ptr<amrex::FabArray<amrex::BaseFab<int> > > >& masks); /// - /// This evolves all the particles by one PIC time step, including charge deposition, the + /// This evolves all the particles by one PIC time step, including charge deposition, the /// field solve, and pushing the particles, for all the species in the MultiParticleContainer. /// This is the electrostatic version. /// - void EvolveES (const amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>, 3> >& E, - amrex::Vector<std::unique_ptr<amrex::MultiFab> >& rho, + void EvolveES (const amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>, 3> >& E, + amrex::Vector<std::unique_ptr<amrex::MultiFab> >& rho, amrex::Real t, amrex::Real dt); /// - /// This pushes the particle positions by one half time step for all the species in the + /// This pushes the particle positions by one half time step for all the species in the /// MultiParticleContainer. It is used to desynchronize the particles after initializaton /// or when restarting from a checkpoint. This is the electrostatic version. - /// + /// void PushXES (amrex::Real dt); /// /// This deposits the particle charge onto rho, accumulating the value for all the species - /// in the MultiParticleContainer. rho is assumed to contain node-centered multifabs. + /// in the MultiParticleContainer. rho is assumed to contain node-centered multifabs. /// This version is hard-coded for CIC deposition. /// void DepositCharge(amrex::Vector<std::unique_ptr<amrex::MultiFab> >& rho, bool local = false); @@ -78,7 +84,7 @@ public: /// amrex::Real sumParticleCharge(bool local = false); #endif // WARPX_DO_ELECTROSTATIC - + /// /// Performs the field gather operation using the input fields E and B, for all the species /// in the MultiParticleContainer. This is the electromagnetic version of the field gather. @@ -96,8 +102,8 @@ public: void Evolve (int lev, const amrex::MultiFab& Ex, const amrex::MultiFab& Ey, const amrex::MultiFab& Ez, const amrex::MultiFab& Bx, const amrex::MultiFab& By, const amrex::MultiFab& Bz, - amrex::MultiFab& jx, amrex::MultiFab& jy, amrex::MultiFab& jz, - amrex::MultiFab* cjx, amrex::MultiFab* cjy, amrex::MultiFab* cjz, + amrex::MultiFab& jx, amrex::MultiFab& jy, amrex::MultiFab& jz, + amrex::MultiFab* cjx, amrex::MultiFab* cjy, amrex::MultiFab* cjz, amrex::MultiFab* rho, amrex::MultiFab* crho, const amrex::MultiFab* cEx, const amrex::MultiFab* cEy, const amrex::MultiFab* cEz, const amrex::MultiFab* cBx, const amrex::MultiFab* cBy, const amrex::MultiFab* cBz, @@ -107,7 +113,7 @@ public: /// This pushes the particle positions by one half time step for all the species in the /// MultiParticleContainer. It is used to desynchronize the particles after initializaton /// or when restarting from a checkpoint. This is the electromagnetic version. - /// + /// void PushX (amrex::Real dt); /// @@ -115,7 +121,7 @@ public: /// MultiParticleContainer. It is used to desynchronize the particles after initializaton /// or when restarting from a checkpoint. It is also used to synchronize particles at the /// the end of the run. This is the electromagnetic version. - /// + /// void PushP (int lev, amrex::Real dt, const amrex::MultiFab& Ex, const amrex::MultiFab& Ey, const amrex::MultiFab& Ez, const amrex::MultiFab& Bx, const amrex::MultiFab& By, const amrex::MultiFab& Bz); @@ -173,7 +179,7 @@ public: const amrex::Real t_boost, const amrex::Real t_lab, const amrex::Real dt, amrex::Vector<WarpXParticleContainer::DiagnosticParticleData>& parts) const; - // Inject particles during the simulation (for particles entering the + // Inject particles during the simulation (for particles entering the // simulation domain after some iterations, due to flowing plasma and/or // moving window). void ContinuousInjection(const amrex::RealBox& injection_box) const; @@ -184,11 +190,16 @@ public: std::vector<std::string> GetSpeciesNames() const { return species_names; } PhysicalParticleContainer& GetPCtmp () { return *pc_tmp; } - + + //bw_engine is a public member of the MultiParticleContainer object +#ifdef WARPX_QED + warpx_breit_wheeler_engine bw_engine; +#endif + protected: // Particle container types - enum struct PCTypes {Physical, RigidInjected}; + enum struct PCTypes {Physical, RigidInjected, Photon}; std::vector<std::string> species_names; @@ -212,7 +223,7 @@ private: // Number of species dumped in BoostedFrameDiagnostics int nspecies_boosted_frame_diags = 0; - // map_species_boosted_frame_diags[i] is the species ID in + // map_species_boosted_frame_diags[i] is the species ID in // MultiParticleContainer for 0<i<nspecies_boosted_frame_diags std::vector<int> map_species_boosted_frame_diags; int do_boosted_frame_diags = 0; diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 143f4b7f9..da09adbf2 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -8,8 +8,16 @@ using namespace amrex; +#ifdef WARPX_QED +MultiParticleContainer::MultiParticleContainer (AmrCore* amr_core): + bw_engine{std::move(init_warpx_breit_wheeler_engine())} +#else MultiParticleContainer::MultiParticleContainer (AmrCore* amr_core) +#endif { + + + ReadParameters(); allcontainers.resize(nspecies + nlasers); @@ -20,9 +28,12 @@ MultiParticleContainer::MultiParticleContainer (AmrCore* amr_core) else if (species_types[i] == PCTypes::RigidInjected) { allcontainers[i].reset(new RigidInjectedParticleContainer(amr_core, i, species_names[i])); } + else if (species_types[i] == PCTypes::Photon) { + allcontainers[i].reset(new PhotonParticleContainer(amr_core, i, species_names[i])); + } allcontainers[i]->deposit_on_main_grid = deposit_on_main_grid[i]; } - + for (int i = nspecies; i < nspecies+nlasers; ++i) { allcontainers[i].reset(new LaserParticleContainer(amr_core,i, lasers_names[i-nspecies])); } @@ -56,9 +67,11 @@ MultiParticleContainer::ReadParameters () BL_ASSERT(nspecies >= 0); if (nspecies > 0) { + // Get species names pp.getarr("species_names", species_names); BL_ASSERT(species_names.size() == nspecies); + // Get species to deposit on main grid deposit_on_main_grid.resize(nspecies, 0); std::vector<std::string> tmp; pp.queryarr("deposit_on_main_grid", tmp); @@ -71,9 +84,9 @@ MultiParticleContainer::ReadParameters () species_types.resize(nspecies, PCTypes::Physical); + // Get rigid-injected species std::vector<std::string> rigid_injected_species; pp.queryarr("rigid_injected_species", rigid_injected_species); - if (!rigid_injected_species.empty()) { for (auto const& name : rigid_injected_species) { auto it = std::find(species_names.begin(), species_names.end(), name); @@ -82,6 +95,18 @@ MultiParticleContainer::ReadParameters () species_types[i] = PCTypes::RigidInjected; } } + // Get photon species + std::vector<std::string> photon_species; + pp.queryarr("photon_species", photon_species); + if (!photon_species.empty()) { + for (auto const& name : photon_species) { + auto it = std::find(species_names.begin(), species_names.end(), name); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(it != species_names.end(), "ERROR: species in particles.rigid_injected_species must be part of particles.species_names"); + int i = std::distance(species_names.begin(), it); + species_types[i] = PCTypes::Photon; + } + } + } pp.query("use_fdtd_nci_corr", WarpX::use_fdtd_nci_corr); @@ -118,6 +143,31 @@ MultiParticleContainer::InitData () // For each species, get the ID of its product species. // This is used for ionization and pair creation processes. mapSpeciesProduct(); + +#ifdef WARPX_QED + + //Helper function + auto does_file_exist = [](const char *fileName) + {return (std::ifstream{fileName}).good(); }; + + //Initialize the lookup tables + //Generates tables if they do not exist + if(!does_file_exist("bw_engine_dndt.bin")){ + bw_engine.compute_dN_dt_lookup_table(&std::cout); + bw_engine.write_dN_dt_table("bw_engine_dndt.bin"); + } + else{ + bw_engine.read_dN_dt_table("bw_engine_dndt.bin"); + } + + if(!does_file_exist("bw_engine_pair.bin")){ + bw_engine.compute_cumulative_pair_table(&std::cout); + bw_engine.write_cumulative_pair_table("bw_engine_pair.bin"); + } + else{ + bw_engine.read_cumulative_pair_table("bw_engine_pair.bin"); + } +#endif } @@ -348,7 +398,7 @@ MultiParticleContainer { BL_PROFILE("MultiParticleContainer::GetLabFrameData"); - + // Loop over particle species for (int i = 0; i < nspecies_boosted_frame_diags; ++i){ int isp = map_species_boosted_frame_diags[i]; @@ -357,41 +407,41 @@ MultiParticleContainer pc->GetParticleSlice(direction, z_old, z_new, t_boost, t_lab, dt, diagnostic_particles); // Here, diagnostic_particles[lev][index] is a WarpXParticleContainer::DiagnosticParticleData // where "lev" is the AMR level and "index" is a [grid index][tile index] pair. - + // Loop over AMR levels for (int lev = 0; lev <= pc->finestLevel(); ++lev){ - // Loop over [grid index][tile index] pairs - // and Fills parts[species number i] with particle data from all grids and - // tiles in diagnostic_particles. parts contains particles from all + // Loop over [grid index][tile index] pairs + // and Fills parts[species number i] with particle data from all grids and + // tiles in diagnostic_particles. parts contains particles from all // AMR levels indistinctly. for (auto it = diagnostic_particles[lev].begin(); it != diagnostic_particles[lev].end(); ++it){ // it->first is the [grid index][tile index] key - // it->second is the corresponding + // it->second is the corresponding // WarpXParticleContainer::DiagnosticParticleData value parts[i].GetRealData(DiagIdx::w).insert( parts[i].GetRealData(DiagIdx::w ).end(), it->second.GetRealData(DiagIdx::w ).begin(), it->second.GetRealData(DiagIdx::w ).end()); - + parts[i].GetRealData(DiagIdx::x).insert( parts[i].GetRealData(DiagIdx::x ).end(), it->second.GetRealData(DiagIdx::x ).begin(), it->second.GetRealData(DiagIdx::x ).end()); - + parts[i].GetRealData(DiagIdx::y).insert( parts[i].GetRealData(DiagIdx::y ).end(), it->second.GetRealData(DiagIdx::y ).begin(), it->second.GetRealData(DiagIdx::y ).end()); - + parts[i].GetRealData(DiagIdx::z).insert( parts[i].GetRealData(DiagIdx::z ).end(), it->second.GetRealData(DiagIdx::z ).begin(), it->second.GetRealData(DiagIdx::z ).end()); - + parts[i].GetRealData(DiagIdx::ux).insert( parts[i].GetRealData(DiagIdx::ux).end(), it->second.GetRealData(DiagIdx::ux).begin(), it->second.GetRealData(DiagIdx::ux).end()); - + parts[i].GetRealData(DiagIdx::uy).insert( parts[i].GetRealData(DiagIdx::uy).end(), it->second.GetRealData(DiagIdx::uy).begin(), it->second.GetRealData(DiagIdx::uy).end()); - + parts[i].GetRealData(DiagIdx::uz).insert( parts[i].GetRealData(DiagIdx::uz).end(), it->second.GetRealData(DiagIdx::uz).begin(), it->second.GetRealData(DiagIdx::uz).end()); @@ -402,7 +452,7 @@ MultiParticleContainer /* \brief Continuous injection for particles initially outside of the domain. * \param injection_box: Domain where new particles should be injected. - * Loop over all WarpXParticleContainer in MultiParticleContainer and + * Loop over all WarpXParticleContainer in MultiParticleContainer and * calls virtual function ContinuousInjection. */ void @@ -418,7 +468,7 @@ MultiParticleContainer::ContinuousInjection(const RealBox& injection_box) const /* \brief Update position of continuous injection parameters. * \param dt: simulation time step (level 0) - * All classes inherited from WarpXParticleContainer do not have + * All classes inherited from WarpXParticleContainer do not have * a position to update (PhysicalParticleContainer does not do anything). */ void diff --git a/Source/Particles/PhotonParticleContainer.H b/Source/Particles/PhotonParticleContainer.H new file mode 100644 index 000000000..966fc1e2d --- /dev/null +++ b/Source/Particles/PhotonParticleContainer.H @@ -0,0 +1,83 @@ +#ifndef WARPX_PhotonParticleContainer_H_ +#define WARPX_PhotonParticleContainer_H_ + +#include <PhysicalParticleContainer.H> +#include <AMReX_Vector.H> + +#ifdef WARPX_QED + #include "breit_wheeler_engine_wrapper.h" +#endif + +class PhotonParticleContainer + : public PhysicalParticleContainer +{ +public: + PhotonParticleContainer (amrex::AmrCore* amr_core, + int ispecies, + const std::string& name); + virtual ~PhotonParticleContainer () {} + + virtual void InitData() override; + + virtual void Evolve (int lev, + const amrex::MultiFab& Ex, + const amrex::MultiFab& Ey, + const amrex::MultiFab& Ez, + const amrex::MultiFab& Bx, + const amrex::MultiFab& By, + const amrex::MultiFab& Bz, + amrex::MultiFab& jx, + amrex::MultiFab& jy, + amrex::MultiFab& jz, + amrex::MultiFab* cjx, + amrex::MultiFab* cjy, + amrex::MultiFab* cjz, + amrex::MultiFab* rho, + amrex::MultiFab* crho, + const amrex::MultiFab* cEx, + const amrex::MultiFab* cEy, + const amrex::MultiFab* cEz, + const amrex::MultiFab* cBx, + const amrex::MultiFab* cBy, + const amrex::MultiFab* cBz, + amrex::Real t, + amrex::Real dt) override; + + virtual void PushPX(WarpXParIter& pti, + amrex::Cuda::ManagedDeviceVector<amrex::Real>& xp, + amrex::Cuda::ManagedDeviceVector<amrex::Real>& yp, + amrex::Cuda::ManagedDeviceVector<amrex::Real>& zp, + amrex::Cuda::ManagedDeviceVector<amrex::Real>& giv, + amrex::Real dt) override; + + virtual void DepositCurrent(WarpXParIter& pti, + RealVector& wp, + RealVector& uxp, + RealVector& uyp, + RealVector& uzp, + amrex::MultiFab& jx, + amrex::MultiFab& jy, + amrex::MultiFab& jz, + amrex::MultiFab* cjx, + amrex::MultiFab* cjy, + amrex::MultiFab* cjz, + const long np_current, + const long np, + int thread_num, + int lev, + amrex::Real dt ) override {}; + +//This function initialises the optical depth of the photons +#ifdef WARPX_QED + virtual void InitOpticalDepth( + WarpXParIter& pti, + warpx_breit_wheeler_engine& engine); +#endif + +private: + + amrex::Real size_in_inches = 1.; + +}; + +#endif // #ifndef WARPX_PhotonParticleContainer_H_ diff --git a/Source/Particles/PhotonParticleContainer.cpp b/Source/Particles/PhotonParticleContainer.cpp new file mode 100644 index 000000000..2bcf87d6a --- /dev/null +++ b/Source/Particles/PhotonParticleContainer.cpp @@ -0,0 +1,115 @@ +#include <limits> +#include <sstream> +#include <algorithm> + +#ifdef _OPENMP +#include <omp.h> +#endif + +#include <PhotonParticleContainer.H> +#include <WarpX_f.H> +#include <WarpX.H> +#include <WarpXConst.H> + +using namespace amrex; + +PhotonParticleContainer::PhotonParticleContainer (AmrCore* amr_core, int ispecies, + const std::string& name) + : PhysicalParticleContainer(amr_core, ispecies, name) +{ + + // This will read <species>.[...] from the inputs file + // where <species> is the name of your species + ParmParse pp(species_name); + + // read <species>.size_in_inches in the input file, and + // store it into member data. + pp.query("size_in_inches", size_in_inches); + +#ifdef WARPX_QED + AddRealComp("tau"); + plot_flags.resize(PIdx::nattribs + 1, 1); +#endif + +} + +void PhotonParticleContainer::InitData() +{ + AddParticles(0); // Note - add on level 0 + + if (maxLevel() > 0) { + Redistribute(); // We then redistribute + } +} + +void +PhotonParticleContainer::PushPX(WarpXParIter& pti, + Cuda::ManagedDeviceVector<Real>& xp, + Cuda::ManagedDeviceVector<Real>& yp, + Cuda::ManagedDeviceVector<Real>& zp, + Cuda::ManagedDeviceVector<Real>& giv, + Real dt) +{ + // This wraps the call to warpx_particle_pusher so that inheritors can modify the call. + auto& attribs = pti.GetAttribs(); + auto& uxp = attribs[PIdx::ux]; + auto& uyp = attribs[PIdx::uy]; + auto& uzp = attribs[PIdx::uz]; + auto& Exp = attribs[PIdx::Ex]; + auto& Eyp = attribs[PIdx::Ey]; + auto& Ezp = attribs[PIdx::Ez]; + auto& Bxp = attribs[PIdx::Bx]; + auto& Byp = attribs[PIdx::By]; + auto& Bzp = attribs[PIdx::Bz]; + const long np = pti.numParticles(); + + // Using new pusher for positions + const amrex_real zero_mass = 0.0; + warpx_particle_pusher_positions(&np, + xp.dataPtr(), + yp.dataPtr(), + zp.dataPtr(), + uxp.dataPtr(), uyp.dataPtr(), uzp.dataPtr(), + giv.dataPtr(), + &zero_mass, &dt, + &WarpX::particle_pusher_algo); +} + +void +PhotonParticleContainer::Evolve (int lev, + const MultiFab& Ex, const MultiFab& Ey, const MultiFab& Ez, + const MultiFab& Bx, const MultiFab& By, const MultiFab& Bz, + MultiFab& jx, MultiFab& jy, MultiFab& jz, + MultiFab* cjx, MultiFab* cjy, MultiFab* cjz, + MultiFab* rho, MultiFab* crho, + const MultiFab* cEx, const MultiFab* cEy, const MultiFab* cEz, + const MultiFab* cBx, const MultiFab* cBy, const MultiFab* cBz, + Real t, Real dt) +{ + // This does gather, push and depose. + // Push and depose have been re-written for photon, + // so they do not do anything. + // Currently, I guess photons do gather fields from the mesh. + PhysicalParticleContainer::Evolve (lev, + Ex, Ey, Ez, + Bx, By, Bz, + jx, jy, jz, + cjx, cjy, cjz, + rho, crho, + cEx, cEy, cEz, + cBx, cBy, cBz, + t, dt); + +} + + +#ifdef WARPX_QED +void PhotonParticleContainer::InitOpticalDepth( + WarpXParIter& pti, + warpx_breit_wheeler_engine& engine) +{ + auto& taus = pti.GetAttribs(particle_comps["tau"]); + for(auto& tau: taus) + tau = engine.get_optical_depth(); +} +#endif diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 80d1caa0f..6d5b51264 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -136,7 +136,7 @@ WarpXParticleContainer::AddOneParticle (int lev, int grid, int tile, Real x, Real y, Real z, std::array<Real,PIdx::nattribs>& attribs) { - auto& particle_tile = GetParticles(lev)[std::make_pair(grid,tile)]; + auto& particle_tile = DefineAndReturnParticleTile(lev, grid, tile); AddOneParticle(particle_tile, x, y, z, attribs); } @@ -163,6 +163,11 @@ WarpXParticleContainer::AddOneParticle (ParticleTileType& particle_tile, particle_tile.push_back(p); particle_tile.push_back_real(attribs); + + for (int i = PIdx::nattribs; i < NumRealComps(); ++i) + { + particle_tile.push_back_real(i, 0.0); + } } void @@ -195,7 +200,7 @@ WarpXParticleContainer::AddNParticles (int lev, // Add to grid 0 and tile 0 // Redistribute() will move them to proper places. std::pair<int,int> key {0,0}; - auto& particle_tile = GetParticles(lev)[key]; + auto& particle_tile = DefineAndReturnParticleTile(0, 0, 0); std::size_t np = iend-ibegin; @@ -258,6 +263,11 @@ WarpXParticleContainer::AddNParticles (int lev, particle_tile.push_back_real(comp, np, 0.0); #endif } + + for (int i = PIdx::nattribs; i < NumRealComps(); ++i) + { + particle_tile.push_back_real(i, 0.0); + } } Redistribute(); diff --git a/Source/QED/Make.package b/Source/QED/Make.package new file mode 100644 index 000000000..6e38efb9e --- /dev/null +++ b/Source/QED/Make.package @@ -0,0 +1,6 @@ +CEXE_headers += amrex_rng_wrapper.h +CEXE_headers += breit_wheeler_engine_wrapper.h +CEXE_sources += amrex_rng_wrapper.cpp + +INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/QED +VPATH_LOCATIONS += $(WARPX_HOME)/Source/QED diff --git a/Source/QED/amrex_rng_wrapper.cpp b/Source/QED/amrex_rng_wrapper.cpp new file mode 100644 index 000000000..955144c4c --- /dev/null +++ b/Source/QED/amrex_rng_wrapper.cpp @@ -0,0 +1,18 @@ +//This file provides a wrapper aroud the RNG +//provided by the amrex library + +#include "amrex_rng_wrapper.h" + +//RNG wrapper BW engine +//Get rnd number uniformly distributed in [a,b) +amrex::Real amrex_rng_wrapper::unf(amrex::Real a, amrex::Real b) +{ + return (b-a)*amrex::Random() + a; +} + +//Get rnd number with exponential distribution +amrex::Real amrex_rng_wrapper::exp(amrex::Real l) +{ + amrex::Real zero_plus_to_one = 1.0 - unf(0.0, 1.0); + return -log(zero_plus_to_one)/l; +} diff --git a/Source/QED/amrex_rng_wrapper.h b/Source/QED/amrex_rng_wrapper.h new file mode 100644 index 000000000..f9818c589 --- /dev/null +++ b/Source/QED/amrex_rng_wrapper.h @@ -0,0 +1,20 @@ +#ifndef WARPX_amrex_rng_wrapper_h_ +#define WARPX_rng_wrapper_h_ + +//This file provides a wrapper aroud the RNG +//provided by the amrex library + +#include <AMReX_AmrCore.H> + +//RNG wrapper BW engine +class amrex_rng_wrapper +{ +public: + //Get rnd number uniformly distributed in [a,b) + amrex::Real unf(amrex::Real a, amrex::Real b); + + //Get rnd number with exponential distribution + amrex::Real exp(amrex::Real l); +}; + +#endif //amrex_rng_wrapper_hpp_ diff --git a/Source/QED/breit_wheeler_engine_wrapper.h b/Source/QED/breit_wheeler_engine_wrapper.h new file mode 100644 index 000000000..3bb3ea740 --- /dev/null +++ b/Source/QED/breit_wheeler_engine_wrapper.h @@ -0,0 +1,24 @@ +#ifndef WARPX_breit_wheeler_engine_wrapper_h_ +#define WARPX_breit_wheeler_engine_wrapper_h_ + +//This file provides a wrapper aroud the breit_wheeler engine +//provided by the standard template library + +//BW ENGINE +//#define PXRMP_GPU __host__ __device__ +#define PXRMP_WITH_SI_UNITS +#include "breit_wheeler_engine.hpp" + +#include "amrex_rng_wrapper.h" + +using warpx_breit_wheeler_engine = + picsar::multi_physics::breit_wheeler_engine<amrex::Real, amrex_rng_wrapper>; + +//Helper function to initialize the engine +inline warpx_breit_wheeler_engine init_warpx_breit_wheeler_engine(){ + return warpx_breit_wheeler_engine{std::move(amrex_rng_wrapper{})}; +} + +//___________________________________________ + +#endif //WARPX_breit_wheeler_engine_wrapper_H_ |