diff options
Diffstat (limited to 'Source')
24 files changed, 573 insertions, 40 deletions
diff --git a/Source/Diagnostics/ParticleIO.cpp b/Source/Diagnostics/ParticleIO.cpp index 2a9c16aa8..d55660b39 100644 --- a/Source/Diagnostics/ParticleIO.cpp +++ b/Source/Diagnostics/ParticleIO.cpp @@ -112,6 +112,12 @@ MultiParticleContainer::WritePlotFile (const std::string& dir) const int_flags.resize(1, 1); } +#ifdef WARPX_QED + if(pc->do_qed){ + real_names.push_back("tau"); + } +#endif + // Convert momentum to SI pc->ConvertUnits(ConvertDirection::WarpX_to_SI); // real_names contains a list of all particle attributes. diff --git a/Source/Make.WarpX b/Source/Make.WarpX index c7a4ef752..54baecbf6 100644 --- a/Source/Make.WarpX +++ b/Source/Make.WarpX @@ -4,6 +4,7 @@ OPENBC_HOME ?= ../../../openbc_poisson USE_MPI = TRUE USE_PARTICLES = TRUE +USE_RPATH = TRUE ifeq ($(USE_GPU),TRUE) USE_OMP = FALSE @@ -93,6 +94,7 @@ ifeq ($(QED),TRUE) FFLAGS += -DWARPX_QED F90FLAGS += -DWARPX_QED include $(WARPX_HOME)/Source/QED/Make.package + USERSuffix := $(USERSuffix).QED endif diff --git a/Source/Particles/MultiParticleContainer.H b/Source/Particles/MultiParticleContainer.H index 4ce83685d..58546a106 100644 --- a/Source/Particles/MultiParticleContainer.H +++ b/Source/Particles/MultiParticleContainer.H @@ -9,6 +9,11 @@ #include <PhotonParticleContainer.H> #include <LaserParticleContainer.H> +#ifdef WARPX_QED + #include <BreitWheelerEngineWrapper.H> + #include <QuantumSyncEngineWrapper.H> +#endif + #include <memory> #include <map> #include <string> @@ -208,6 +213,17 @@ protected: std::vector<PCTypes> species_types; +#ifdef WARPX_QED + // The QED engines + BreitWheelerEngine bw_engine; + QuantumSynchrotronEngine qs_engine; + //_______________________________ + + //Initialize QED engines and provides smart pointers + //to species who need QED processes + void InitQED (); +#endif + private: // physical particles (+ laser) diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 715c97b99..c860d21f5 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -1,11 +1,12 @@ -#include <limits> -#include <algorithm> -#include <string> - #include <MultiParticleContainer.H> #include <WarpX_f.H> #include <WarpX.H> +#include <limits> +#include <algorithm> +#include <string> +#include <memory> + using namespace amrex; MultiParticleContainer::MultiParticleContainer (AmrCore* amr_core) @@ -149,6 +150,11 @@ 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 + InitQED(); +#endif + } @@ -726,3 +732,19 @@ MultiParticleContainer::doFieldIonization () } // lev } // pc_source } + +#ifdef WARPX_QED +void MultiParticleContainer::InitQED () +{ + for (auto& pc : allcontainers) { + if(pc->has_quantum_sync()){ + pc->set_quantum_sync_engine_ptr + (std::make_shared<QuantumSynchrotronEngine>(qs_engine)); + } + if(pc->has_breit_wheeler()){ + pc->set_breit_wheeler_engine_ptr + (std::make_shared<BreitWheelerEngine>(bw_engine)); + } + } +} +#endif diff --git a/Source/Particles/PhotonParticleContainer.H b/Source/Particles/PhotonParticleContainer.H index 9b53a04ba..407cf26f3 100644 --- a/Source/Particles/PhotonParticleContainer.H +++ b/Source/Particles/PhotonParticleContainer.H @@ -46,6 +46,15 @@ public: amrex::Cuda::ManagedDeviceVector<amrex::ParticleReal>& zp, amrex::Real dt, DtType a_dt_type=DtType::Full) override; + // Don't push momenta for photons + virtual 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) override {}; // DepositCurrent should do nothing for photons @@ -64,7 +73,6 @@ public: int lev, int depos_lev, amrex::Real dt) override {}; - }; #endif // #ifndef WARPX_PhotonParticleContainer_H_ diff --git a/Source/Particles/PhotonParticleContainer.cpp b/Source/Particles/PhotonParticleContainer.cpp index 4a75ec9f3..3c70a957f 100644 --- a/Source/Particles/PhotonParticleContainer.cpp +++ b/Source/Particles/PhotonParticleContainer.cpp @@ -21,7 +21,25 @@ using namespace amrex; PhotonParticleContainer::PhotonParticleContainer (AmrCore* amr_core, int ispecies, const std::string& name) : PhysicalParticleContainer(amr_core, ispecies, name) -{} +{ + + ParmParse pp(species_name); + +#ifdef WARPX_QED + //IF do_qed is enabled, find out if Breit Wheeler process is enabled + if(do_qed) + pp.query("do_qed_breit_wheeler", do_qed_breit_wheeler); + + //Check for processes which do not make sense for photons + bool test_quantum_sync; + pp.query("do_qed_quantum_sync", test_quantum_sync); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + test_quantum_sync == 0, + "ERROR: do_qed_quantum_sync can be 1 for species NOT listed in particles.photon_species only!"); + //_________________________________________________________ +#endif + +} void PhotonParticleContainer::InitData() { diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index b4081e959..c953aa2d7 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -7,6 +7,11 @@ #include <AMReX_IArrayBox.H> +#ifdef WARPX_QED + #include <QuantumSyncEngineWrapper.H> + #include <BreitWheelerEngineWrapper.H> +#endif + #include <map> class PhysicalParticleContainer @@ -180,6 +185,16 @@ public: amrex::FArrayBox const * & ezfab, amrex::FArrayBox const * & bxfab, amrex::FArrayBox const * & byfab, amrex::FArrayBox const * & bzfab); +#ifdef WARPX_QED + bool has_quantum_sync() override; + bool has_breit_wheeler() override; + + void set_breit_wheeler_engine_ptr + (std::shared_ptr<BreitWheelerEngine>) override; + void set_quantum_sync_engine_ptr + (std::shared_ptr<QuantumSynchrotronEngine>) override; +#endif + protected: std::string species_name; @@ -194,6 +209,20 @@ protected: // Inject particles during the whole simulation void ContinuousInjection (const amrex::RealBox& injection_box) override; +#ifdef WARPX_QED + // A flag to enable quantum_synchrotron process for leptons + bool do_qed_quantum_sync = false; + + // A flag to enable breit_wheeler process [photons only!!] + bool do_qed_breit_wheeler = false; + + // A smart pointer to an instance of a Quantum Synchrotron engine + std::shared_ptr<QuantumSynchrotronEngine> shr_ptr_qs_engine; + + // A smart pointer to an instance of a Breit Wheeler engine [photons only!] + std::shared_ptr<BreitWheelerEngine> shr_ptr_bw_engine; +#endif + }; #endif diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 75e438454..66331db26 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1,6 +1,8 @@ #include <limits> #include <sstream> +#include <PhysicalParticleContainer.H> + #include <MultiParticleContainer.H> #include <WarpX_f.H> #include <WarpX.H> @@ -15,6 +17,7 @@ #include <UpdatePosition.H> #include <UpdateMomentumBoris.H> #include <UpdateMomentumVay.H> +#include <UpdateMomentumHigueraCary.H> using namespace amrex; @@ -49,6 +52,33 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp "version of Redistribute in AMReX does not work with runtime parameters"); #endif + +#ifdef WARPX_QED + //Add real component if QED is enabled + pp.query("do_qed", do_qed); + if(do_qed) + AddRealComp("tau"); + + //IF do_qed is enabled, find out if Quantum Synchrotron process is enabled + if(do_qed) + pp.query("do_qed_quantum_sync", do_qed_quantum_sync); + + //TODO: SHOULD CHECK IF SPECIES IS EITHER ELECTRONS OR POSITRONS!! +#endif + + //variable to set plot_flags size + int plot_flag_size = PIdx::nattribs; + if(WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) + plot_flag_size += 6; + +#ifdef WARPX_QED + if(do_qed){ + // plot_flag will have an entry for the optical depth + plot_flag_size++; + } +#endif + //_______________________________ + pp.query("plot_species", plot_species); int do_user_plot_vars; do_user_plot_vars = pp.queryarr("plot_vars", plot_vars); @@ -56,18 +86,11 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp // By default, all particle variables are dumped to plotfiles, // including {x,y,z,ux,uy,uz}old variables when running in a // boosted frame - if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags){ - plot_flags.resize(PIdx::nattribs + 6, 1); - } else { - plot_flags.resize(PIdx::nattribs, 1); - } + plot_flags.resize(plot_flag_size, 1); } else { // Set plot_flag to 0 for all attribs - if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags){ - plot_flags.resize(PIdx::nattribs + 6, 0); - } else { - plot_flags.resize(PIdx::nattribs, 0); - } + plot_flags.resize(plot_flag_size, 0); + // If not none, set plot_flags values to 1 for elements in plot_vars. if (plot_vars[0] != "none"){ for (const auto& var : plot_vars){ @@ -79,6 +102,13 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp } } } + + #ifdef WARPX_QED + if(do_qed){ + //Optical depths is always plotted if QED is on + plot_flags[plot_flag_size-1] = 1; + } + #endif } PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core) @@ -92,7 +122,9 @@ void PhysicalParticleContainer::InitData() // Init ionization module here instead of in the PhysicalParticleContainer // constructor because dt is required if (do_field_ionization) {InitIonizationModule();} + AddParticles(0); // Note - add on level 0 + Redistribute(); // We then redistribute } @@ -477,6 +509,30 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) pi = soa.GetIntData(particle_icomps["ionization_level"]).data() + old_size; } +#ifdef WARPX_QED + //Pointer to the optical depth component + amrex::Real* p_tau; + + //If a QED effect is enabled, tau has to be initialized + bool loc_has_quantum_sync = has_quantum_sync(); + bool loc_has_breit_wheeler = has_breit_wheeler(); + if(loc_has_quantum_sync || loc_has_breit_wheeler){ + p_tau = soa.GetRealData(particle_comps["tau"]).data() + old_size; + } + + //If needed, get the appropriate functors from the engines + QuantumSynchrotronGetOpticalDepth quantum_sync_get_opt; + BreitWheelerGetOpticalDepth breit_wheeler_get_opt; + if(loc_has_quantum_sync){ + quantum_sync_get_opt = + shr_ptr_qs_engine->build_optical_depth_functor(); + } + if(loc_has_breit_wheeler){ + breit_wheeler_get_opt = + shr_ptr_bw_engine->build_optical_depth_functor(); + } +#endif + const GpuArray<Real,AMREX_SPACEDIM> overlap_corner {AMREX_D_DECL(overlap_realbox.lo(0), overlap_realbox.lo(1), @@ -624,6 +680,16 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) pi[ip] = loc_ionization_initial_level; } +#ifdef WARPX_QED + if(loc_has_quantum_sync){ + p_tau[ip] = quantum_sync_get_opt(); + } + + if(loc_has_breit_wheeler){ + p_tau[ip] = breit_wheeler_get_opt(); + } +#endif + u.x *= PhysConst::c; u.y *= PhysConst::c; u.z *= PhysConst::c; @@ -928,12 +994,12 @@ PhysicalParticleContainer::FieldGather (int lev, const FArrayBox& byfab = By[pti]; const FArrayBox& bzfab = Bz[pti]; - Exp.assign(np,0.0); - Eyp.assign(np,0.0); - Ezp.assign(np,0.0); - Bxp.assign(np,0.0); - Byp.assign(np,0.0); - Bzp.assign(np,0.0); + Exp.assign(np,WarpX::E_external[0]); + Eyp.assign(np,WarpX::E_external[1]); + Ezp.assign(np,WarpX::E_external[2]); + Bxp.assign(np,WarpX::B_external[0]); + Byp.assign(np,WarpX::B_external[1]); + Bzp.assign(np,WarpX::B_external[2]); // // copy data from particle container to temp arrays @@ -1064,9 +1130,10 @@ PhysicalParticleContainer::Evolve (int lev, exfab, eyfab, ezfab, bxfab, byfab, bzfab); } - Exp.assign(np,0.0); - Eyp.assign(np,0.0); - Ezp.assign(np,0.0); + Exp.assign(np,WarpX::E_external[0]); + Eyp.assign(np,WarpX::E_external[1]); + Ezp.assign(np,WarpX::E_external[2]); + Bxp.assign(np,WarpX::B_external[0]); Byp.assign(np,WarpX::B_external[1]); Bzp.assign(np,WarpX::B_external[2]); @@ -1539,6 +1606,19 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti, ux[i], uy[i], uz[i], dt ); } ); + } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::HigueraCary) { + amrex::ParallelFor( + pti.numParticles(), + [=] AMREX_GPU_DEVICE (long i) { + Real qp = q; + if (ion_lev){ qp *= ion_lev[i]; } + UpdateMomentumHigueraCary( ux[i], uy[i], uz[i], + Ex[i], Ey[i], Ez[i], Bx[i], + By[i], Bz[i], qp, m, dt); + UpdatePosition( x[i], y[i], z[i], + ux[i], uy[i], uz[i], dt ); + } + ); } else { amrex::Abort("Unknown particle pusher"); }; @@ -1587,9 +1667,10 @@ PhysicalParticleContainer::PushP (int lev, Real dt, const FArrayBox& byfab = By[pti]; const FArrayBox& bzfab = Bz[pti]; - Exp.assign(np,0.0); - Eyp.assign(np,0.0); - Ezp.assign(np,0.0); + Exp.assign(np,WarpX::E_external[0]); + Eyp.assign(np,WarpX::E_external[1]); + Ezp.assign(np,WarpX::E_external[2]); + Bxp.assign(np,WarpX::B_external[0]); Byp.assign(np,WarpX::B_external[1]); Bzp.assign(np,WarpX::B_external[2]); @@ -1634,6 +1715,13 @@ PhysicalParticleContainer::PushP (int lev, Real dt, Expp[i], Eypp[i], Ezpp[i], Bxpp[i], Bypp[i], Bzpp[i], q, m, dt); } ); + } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::HigueraCary) { + amrex::ParallelFor( pti.numParticles(), + [=] AMREX_GPU_DEVICE (long i) { + UpdateMomentumHigueraCary( ux[i], uy[i], uz[i], + Expp[i], Eypp[i], Ezpp[i], Bxpp[i], Bypp[i], Bzpp[i], q, m, dt); + } + ); } else { amrex::Abort("Unknown particle pusher"); }; @@ -2074,3 +2162,30 @@ PhysicalParticleContainer::buildIonizationMask (const amrex::MFIter& mfi, const } ); } + +#ifdef WARPX_QED + +bool PhysicalParticleContainer::has_quantum_sync() +{ + return do_qed_quantum_sync; +} + +bool PhysicalParticleContainer::has_breit_wheeler() +{ + return do_qed_breit_wheeler; +} + +void +PhysicalParticleContainer:: +set_breit_wheeler_engine_ptr(std::shared_ptr<BreitWheelerEngine> ptr) +{ + shr_ptr_bw_engine = ptr; +} + +void +PhysicalParticleContainer:: +set_quantum_sync_engine_ptr(std::shared_ptr<QuantumSynchrotronEngine> ptr) +{ + shr_ptr_qs_engine = ptr; +} +#endif diff --git a/Source/Particles/Pusher/Make.package b/Source/Particles/Pusher/Make.package index 45886702e..b033bfb57 100644 --- a/Source/Particles/Pusher/Make.package +++ b/Source/Particles/Pusher/Make.package @@ -2,6 +2,7 @@ CEXE_headers += GetAndSetPosition.H CEXE_headers += UpdatePosition.H CEXE_headers += UpdateMomentumBoris.H CEXE_headers += UpdateMomentumVay.H +CEXE_headers += UpdateMomentumHigueraCary.H CEXE_headers += UpdatePositionPhoton.H INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Particles/Pusher VPATH_LOCATIONS += $(WARPX_HOME)/Source/Particles/Pusher diff --git a/Source/Particles/Pusher/UpdateMomentumHigueraCary.H b/Source/Particles/Pusher/UpdateMomentumHigueraCary.H new file mode 100644 index 000000000..51d7fd620 --- /dev/null +++ b/Source/Particles/Pusher/UpdateMomentumHigueraCary.H @@ -0,0 +1,59 @@ +#ifndef WARPX_PARTICLES_PUSHER_UPDATEMOMENTUM_HIGUERACARY_H_ +#define WARPX_PARTICLES_PUSHER_UPDATEMOMENTUM_HIGUERACARY_H_ + +#include "WarpXConst.H" +#include <AMReX_FArrayBox.H> +#include <AMReX_REAL.H> + +/** + * \brief Push the particle's positions over one timestep, + * given the value of its momenta `ux`, `uy`, `uz` + */ + +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void UpdateMomentumHigueraCary( + amrex::ParticleReal& ux, amrex::ParticleReal& uy, amrex::ParticleReal& uz, + const amrex::ParticleReal Ex, const amrex::ParticleReal Ey, const amrex::ParticleReal Ez, + const amrex::ParticleReal Bx, const amrex::ParticleReal By, const amrex::ParticleReal Bz, + const amrex::Real q, const amrex::Real m, const amrex::Real dt ) +{ + // Constants + const amrex::Real qmt = 0.5*q*dt/m; + constexpr amrex::Real invclight = 1./PhysConst::c; + constexpr amrex::Real invclightsq = 1./(PhysConst::c*PhysConst::c); + // Compute u_minus + const amrex::Real umx = ux + qmt*Ex; + const amrex::Real umy = uy + qmt*Ey; + const amrex::Real umz = uz + qmt*Ez; + // Compute gamma squared of u_minus + amrex::Real gamma = 1. + (umx*umx + umy*umy + umz*umz)*invclightsq; + // Compute beta and betam squared + const amrex::Real betax = qmt*Bx; + const amrex::Real betay = qmt*By; + const amrex::Real betaz = qmt*Bz; + const amrex::Real betam = betax*betax + betay*betay + betaz*betaz; + // Compute sigma + const amrex::Real sigma = gamma - betam; + // Get u* + const amrex::Real ust = (umx*betax + umy*betay + umz*betaz)*invclight; + // Get new gamma inversed + gamma = 1./std::sqrt(0.5*(sigma + std::sqrt(sigma*sigma + 4.*(betam + ust*ust)) )); + // Compute t + const amrex::Real tx = gamma*betax; + const amrex::Real ty = gamma*betay; + const amrex::Real tz = gamma*betaz; + // Compute s + const amrex::Real s = 1./(1.+(tx*tx + ty*ty + tz*tz)); + // Compute um dot t + const amrex::Real umt = umx*tx + umy*ty + umz*tz; + // Compute u_plus + const amrex::Real upx = s*( umx + umt*tx + umy*tz - umz*ty ); + const amrex::Real upy = s*( umy + umt*ty + umz*tx - umx*tz ); + const amrex::Real upz = s*( umz + umt*tz + umx*ty - umy*tx ); + // Get new u + ux = upx + qmt*Ex + upy*tz - upz*ty; + uy = upy + qmt*Ey + upz*tx - upx*tz; + uz = upz + qmt*Ez + upx*ty - upy*tx; +} + +#endif // WARPX_PARTICLES_PUSHER_UPDATEMOMENTUM_HIGUERACARY_H_ diff --git a/Source/Particles/RigidInjectedParticleContainer.cpp b/Source/Particles/RigidInjectedParticleContainer.cpp index 891ade76d..2ef833151 100644 --- a/Source/Particles/RigidInjectedParticleContainer.cpp +++ b/Source/Particles/RigidInjectedParticleContainer.cpp @@ -13,6 +13,7 @@ #include <WarpXAlgorithmSelection.H> #include <UpdateMomentumBoris.H> #include <UpdateMomentumVay.H> +#include <UpdateMomentumHigueraCary.H> using namespace amrex; @@ -390,9 +391,9 @@ RigidInjectedParticleContainer::PushP (int lev, Real dt, const FArrayBox& byfab = By[pti]; const FArrayBox& bzfab = Bz[pti]; - Exp.assign(np,0.0); - Eyp.assign(np,0.0); - Ezp.assign(np,0.0); + Exp.assign(np,WarpX::E_external[0]); + Eyp.assign(np,WarpX::E_external[1]); + Ezp.assign(np,WarpX::E_external[2]); Bxp.assign(np,WarpX::B_external[0]); Byp.assign(np,WarpX::B_external[1]); Bzp.assign(np,WarpX::B_external[2]); @@ -443,6 +444,13 @@ RigidInjectedParticleContainer::PushP (int lev, Real dt, Expp[i], Eypp[i], Ezpp[i], Bxpp[i], Bypp[i], Bzpp[i], q, m, dt); } ); + } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::HigueraCary) { + amrex::ParallelFor( pti.numParticles(), + [=] AMREX_GPU_DEVICE (long i) { + UpdateMomentumHigueraCary( uxpp[i], uypp[i], uzpp[i], + Expp[i], Eypp[i], Ezpp[i], Bxpp[i], Bypp[i], Bzpp[i], q, m, dt); + } + ); } else { amrex::Abort("Unknown particle pusher"); }; diff --git a/Source/Particles/Sorting/Partition.cpp b/Source/Particles/Sorting/Partition.cpp index d1455249a..e88af017f 100644 --- a/Source/Particles/Sorting/Partition.cpp +++ b/Source/Particles/Sorting/Partition.cpp @@ -56,7 +56,7 @@ PhysicalParticleContainer::PartitionParticlesInBuffers( gather_masks : current_masks; // - For each particle, find whether it is in the larger buffer, // by looking up the mask. Store the answer in `inexflag`. - amrex::ParallelFor( np, fillBufferFlag(pti, bmasks, inexflag, Geom(lev), 0) ); + amrex::ParallelFor( np, fillBufferFlag(pti, bmasks, inexflag, Geom(lev)) ); // - Find the indices that reorder particles so that the last particles // are in the larger buffer fillWithConsecutiveIntegers( pid ); @@ -93,7 +93,7 @@ PhysicalParticleContainer::PartitionParticlesInBuffers( // - For each particle in the large buffer, find whether it is in // the smaller buffer, by looking up the mask. Store the answer in `inexflag`. amrex::ParallelFor( np - n_fine, - fillBufferFlag(pti, bmasks, inexflag, Geom(lev), n_fine) ); + fillBufferFlagRemainingParticles(pti, bmasks, inexflag, Geom(lev), pid, n_fine) ); auto const sep2 = stablePartition( sep, pid.end(), inexflag ); if (bmasks == gather_masks) { diff --git a/Source/Particles/Sorting/SortingUtils.H b/Source/Particles/Sorting/SortingUtils.H index 852b8a73f..28a1b73b7 100644 --- a/Source/Particles/Sorting/SortingUtils.H +++ b/Source/Particles/Sorting/SortingUtils.H @@ -78,27 +78,84 @@ int iteratorDistance(ForwardIterator const first, /* \brief Functor that fills the elements of the particle array `inexflag` * with the value of the spatial array `bmasks`, at the corresponding particle position. * - * This is done only for the elements from `start_index` to the end of `inexflag` - * * \param[in] pti Contains information on the particle positions * \param[in] bmasks Spatial array, that contains a flag indicating * whether each cell is part of the gathering/deposition buffers * \param[out] inexflag Vector to be filled with the value of `bmasks` * \param[in] geom Geometry object, necessary to locate particles within the array `bmasks` - * \param[in] start_index Index that which elements start to be modified + * */ class fillBufferFlag { public: fillBufferFlag( WarpXParIter const& pti, amrex::iMultiFab const* bmasks, amrex::Gpu::DeviceVector<int>& inexflag, - amrex::Geometry const& geom, long const start_index=0 ) : + amrex::Geometry const& geom ) { + + // Extract simple structure that can be used directly on the GPU + m_particles = &(pti.GetArrayOfStructs()[0]); + m_buffer_mask = (*bmasks)[pti].array(); + m_inexflag_ptr = inexflag.dataPtr(); + m_domain = geom.Domain(); + for (int idim=0; idim<AMREX_SPACEDIM; idim++) { + m_prob_lo[idim] = geom.ProbLo(idim); + m_inv_cell_size[idim] = geom.InvCellSize(idim); + } + }; + + + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + void operator()( const long i ) const { + // Select a particle + auto const& p = m_particles[i]; + // Find the index of the cell where this particle is located + amrex::IntVect const iv = amrex::getParticleCell( p, + m_prob_lo, m_inv_cell_size, m_domain ); + // Find the value of the buffer flag in this cell and + // store it at the corresponding particle position in the array `inexflag` + m_inexflag_ptr[i] = m_buffer_mask(iv); + }; + + private: + amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> m_prob_lo; + amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> m_inv_cell_size; + amrex::Box m_domain; + int* m_inexflag_ptr; + WarpXParticleContainer::ParticleType const* m_particles; + amrex::Array4<int const> m_buffer_mask; +}; + +/* \brief Functor that fills the elements of the particle array `inexflag` + * with the value of the spatial array `bmasks`, at the corresponding particle position. + * + * Contrary to `fillBufferFlag`, here this is done only for the particles that + * the last elements of `particle_indices` point to (from the element at + * index `start_index` in `particle_indices`, to the last element of `particle_indices`) + * + * \param[in] pti Contains information on the particle positions + * \param[in] bmasks Spatial array, that contains a flag indicating + * whether each cell is part of the gathering/deposition buffers + * \param[out] inexflag Vector to be filled with the value of `bmasks` + * \param[in] geom Geometry object, necessary to locate particles within the array `bmasks` + * \param[in] start_index Index that which elements start to be modified + */ +class fillBufferFlagRemainingParticles +{ + public: + fillBufferFlagRemainingParticles( + WarpXParIter const& pti, + amrex::iMultiFab const* bmasks, + amrex::Gpu::DeviceVector<int>& inexflag, + amrex::Geometry const& geom, + amrex::Gpu::DeviceVector<long> const& particle_indices, + long const start_index ) : m_start_index(start_index) { // Extract simple structure that can be used directly on the GPU m_particles = &(pti.GetArrayOfStructs()[0]); m_buffer_mask = (*bmasks)[pti].array(); m_inexflag_ptr = inexflag.dataPtr(); + m_indices_ptr = particle_indices.dataPtr(); m_domain = geom.Domain(); for (int idim=0; idim<AMREX_SPACEDIM; idim++) { m_prob_lo[idim] = geom.ProbLo(idim); @@ -110,13 +167,13 @@ class fillBufferFlag AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void operator()( const long i ) const { // Select a particle - auto const& p = m_particles[i+m_start_index]; + auto const& p = m_particles[m_indices_ptr[i+m_start_index]]; // Find the index of the cell where this particle is located amrex::IntVect const iv = amrex::getParticleCell( p, m_prob_lo, m_inv_cell_size, m_domain ); // Find the value of the buffer flag in this cell and // store it at the corresponding particle position in the array `inexflag` - m_inexflag_ptr[i+m_start_index] = m_buffer_mask(iv); + m_inexflag_ptr[m_indices_ptr[i+m_start_index]] = m_buffer_mask(iv); }; private: @@ -127,6 +184,7 @@ class fillBufferFlag WarpXParticleContainer::ParticleType const* m_particles; amrex::Array4<int const> m_buffer_mask; long const m_start_index; + long const* m_indices_ptr; }; /* \brief Functor that copies the elements of `src` into `dst`, diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index ceb88d024..879f4782e 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -6,6 +6,11 @@ #include <AMReX_Particles.H> #include <AMReX_AmrCore.H> +#ifdef WARPX_QED + #include <QuantumSyncEngineWrapper.H> + #include <BreitWheelerEngineWrapper.H> +#endif + #include <memory> enum struct ConvertDirection{WarpX_to_SI, SI_to_WarpX}; @@ -313,6 +318,18 @@ protected: int do_boosted_frame_diags = 1; +#ifdef WARPX_QED + bool do_qed = false; + + virtual bool has_quantum_sync(){return false;}; + virtual bool has_breit_wheeler(){return false;}; + + virtual void + set_breit_wheeler_engine_ptr(std::shared_ptr<BreitWheelerEngine>){}; + virtual void + set_quantum_sync_engine_ptr(std::shared_ptr<QuantumSynchrotronEngine>){}; +#endif + amrex::Vector<amrex::FArrayBox> local_rho; amrex::Vector<amrex::FArrayBox> local_jx; amrex::Vector<amrex::FArrayBox> local_jy; diff --git a/Source/QED/BreitWheelerEngineWrapper.H b/Source/QED/BreitWheelerEngineWrapper.H new file mode 100644 index 000000000..7cdc66b44 --- /dev/null +++ b/Source/QED/BreitWheelerEngineWrapper.H @@ -0,0 +1,56 @@ +#ifndef WARPX_breit_wheeler_engine_wrapper_h_ +#define WARPX_breit_wheeler_engine_wrapper_h_ + +#include "QedWrapperCommons.H" + +//BW ENGINE from PICSAR +#define PXRMP_CORE_ONLY +#include "breit_wheeler_engine.hpp" + +using WarpXBreitWheelerWrapper = + picsar::multi_physics::breit_wheeler_engine<amrex::Real, DummyStruct>; + +// Functors ================================== + +// These functors provide the core elementary functions of the library +// Can be included in GPU kernels + +/** + * \brief Functor to initialize the optical depth of photons for the + * Breit-Wheeler process + */ +class BreitWheelerGetOpticalDepth +{ +public: + BreitWheelerGetOpticalDepth () + {}; + + AMREX_GPU_DEVICE + AMREX_FORCE_INLINE + amrex::Real operator() () const + { + return WarpXBreitWheelerWrapper:: + internal_get_optical_depth(amrex::Random()); + } +}; +//____________________________________________ + +// Factory class ============================= + +/** + * \brief Wrapper for the Breit Wheeler engine of the PICSAR library + */ +class BreitWheelerEngine +{ +public: + BreitWheelerEngine (); + + /** + * \brief Builds the functor to initialize the optical depth + */ + BreitWheelerGetOpticalDepth build_optical_depth_functor (); +}; + +//============================================ + +#endif //WARPX_breit_wheeler_engine_wrapper_H_ diff --git a/Source/QED/BreitWheelerEngineWrapper.cpp b/Source/QED/BreitWheelerEngineWrapper.cpp new file mode 100644 index 000000000..33e6d48d7 --- /dev/null +++ b/Source/QED/BreitWheelerEngineWrapper.cpp @@ -0,0 +1,16 @@ +#include "BreitWheelerEngineWrapper.H" +//This file provides a wrapper aroud the breit_wheeler engine +//provided by the PICSAR library + +using namespace picsar::multi_physics; + +// Factory class ============================= + +BreitWheelerEngine::BreitWheelerEngine (){} + +BreitWheelerGetOpticalDepth BreitWheelerEngine::build_optical_depth_functor () +{ + return BreitWheelerGetOpticalDepth(); +} + +//============================================ diff --git a/Source/QED/Make.package b/Source/QED/Make.package new file mode 100644 index 000000000..d4bad3f18 --- /dev/null +++ b/Source/QED/Make.package @@ -0,0 +1,8 @@ +CEXE_headers += QedWrapperCommons.H +CEXE_headers += BreitWheelerEngineWrapper.H +CEXE_headers += QuantumSyncsEngineWrapper.H +CEXE_sources += BreitWheelerEngineWrapper.cpp +CEXE_sources += QuantumSyncEngineWrapper.cpp + +INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/QED +VPATH_LOCATIONS += $(WARPX_HOME)/Source/QED diff --git a/Source/QED/QedWrapperCommons.H b/Source/QED/QedWrapperCommons.H new file mode 100644 index 000000000..821034c06 --- /dev/null +++ b/Source/QED/QedWrapperCommons.H @@ -0,0 +1,16 @@ +#ifndef WARPX_amrex_qed_wrapper_commons_h_ +#define WARPX_amrex_qed_wrapper_commons_h_ + +//Common definitions for the QED library wrappers + +#include <AMReX_AmrCore.H> + +//Sets the decorator for GPU +#define PXRMP_GPU AMREX_GPU_DEVICE +//Sets SI units in the library +#define PXRMP_WITH_SI_UNITS + +//An empty data type +struct DummyStruct{}; + +#endif //WARPX_amrex_qed_wrapper_commons_h_ diff --git a/Source/QED/QuantumSyncEngineWrapper.H b/Source/QED/QuantumSyncEngineWrapper.H new file mode 100644 index 000000000..a285dc8b6 --- /dev/null +++ b/Source/QED/QuantumSyncEngineWrapper.H @@ -0,0 +1,56 @@ +#ifndef WARPX_quantum_sync_engine_wrapper_h_ +#define WARPX_quantum_sync_engine_wrapper_h_ + +#include "QedWrapperCommons.H" + +//QS ENGINE from PICSAR +#define PXRMP_CORE_ONLY +#include "quantum_sync_engine.hpp" + +using WarpXQuantumSynchrotronWrapper = + picsar::multi_physics::quantum_synchrotron_engine<amrex::Real, DummyStruct>; + +// Functors ================================== + +// These functors provide the core elementary functions of the library +// Can be included in GPU kernels + +/** +* Functor to initialize the optical depth of leptons for the +* Quantum Synchrotron process +*/ +class QuantumSynchrotronGetOpticalDepth +{ +public: + QuantumSynchrotronGetOpticalDepth () + {}; + + AMREX_GPU_DEVICE + AMREX_FORCE_INLINE + amrex::Real operator() () const + { + return WarpXQuantumSynchrotronWrapper:: + internal_get_optical_depth(amrex::Random()); + } +}; +//____________________________________________ + +// Factory class ============================= + +/** + * \brief Wrapper for the Quantum Synchrotron engine of the PICSAR library + */ +class QuantumSynchrotronEngine +{ +public: + QuantumSynchrotronEngine (); + + /** + * \brief Builds the functor to initialize the optical depth + */ + QuantumSynchrotronGetOpticalDepth build_optical_depth_functor (); +}; + +//============================================ + +#endif //WARPX_quantum_sync_engine_wrapper_h_ diff --git a/Source/QED/QuantumSyncEngineWrapper.cpp b/Source/QED/QuantumSyncEngineWrapper.cpp new file mode 100644 index 000000000..b55149187 --- /dev/null +++ b/Source/QED/QuantumSyncEngineWrapper.cpp @@ -0,0 +1,17 @@ +#include "QuantumSyncEngineWrapper.H" +//This file provides a wrapper aroud the quantum_sync engine +//provided by the PICSAR library + +using namespace picsar::multi_physics; + +// Factory class ============================= + +QuantumSynchrotronEngine::QuantumSynchrotronEngine (){} + +QuantumSynchrotronGetOpticalDepth +QuantumSynchrotronEngine::build_optical_depth_functor () +{ + return QuantumSynchrotronGetOpticalDepth(); +} + +//============================================ diff --git a/Source/Utils/WarpXAlgorithmSelection.H b/Source/Utils/WarpXAlgorithmSelection.H index 54c721abf..269171a5f 100644 --- a/Source/Utils/WarpXAlgorithmSelection.H +++ b/Source/Utils/WarpXAlgorithmSelection.H @@ -14,7 +14,8 @@ struct MaxwellSolverAlgo { struct ParticlePusherAlgo { enum { Boris = 0, - Vay = 1 + Vay = 1, + HigueraCary = 2 }; }; diff --git a/Source/Utils/WarpXAlgorithmSelection.cpp b/Source/Utils/WarpXAlgorithmSelection.cpp index 8a6ff6dbf..be9cdd118 100644 --- a/Source/Utils/WarpXAlgorithmSelection.cpp +++ b/Source/Utils/WarpXAlgorithmSelection.cpp @@ -18,6 +18,7 @@ const std::map<std::string, int> maxwell_solver_algo_to_int = { const std::map<std::string, int> particle_pusher_algo_to_int = { {"boris", ParticlePusherAlgo::Boris }, {"vay", ParticlePusherAlgo::Vay }, + {"higuera", ParticlePusherAlgo::HigueraCary }, {"default", ParticlePusherAlgo::Boris } }; diff --git a/Source/WarpX.H b/Source/WarpX.H index 5e400a029..29b89686e 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -75,6 +75,7 @@ public: // External fields static amrex::Vector<amrex::Real> B_external; + static amrex::Vector<amrex::Real> E_external; // Algorithms static long current_deposition_algo; diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 5c3c94c5f..d94541f17 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -25,6 +25,7 @@ using namespace amrex; Vector<Real> WarpX::B_external(3, 0.0); +Vector<Real> WarpX::E_external(3, 0.0); int WarpX::do_moving_window = 0; int WarpX::moving_window_dir = -1; @@ -294,6 +295,7 @@ WarpX::ReadParameters () zmax_plasma_to_compute_max_step); pp.queryarr("B_external", B_external); + pp.queryarr("E_external", E_external); pp.query("do_moving_window", do_moving_window); if (do_moving_window) |