aboutsummaryrefslogtreecommitdiff
path: root/Source/Particles
diff options
context:
space:
mode:
Diffstat (limited to 'Source/Particles')
-rw-r--r--Source/Particles/Deposition/CurrentDeposition.H121
-rw-r--r--Source/Particles/Make.package1
-rw-r--r--Source/Particles/MultiParticleContainer.H25
-rw-r--r--Source/Particles/MultiParticleContainer.cpp181
-rw-r--r--Source/Particles/ParticleCreation/CopyParticle.H90
-rw-r--r--Source/Particles/ParticleCreation/ElementaryProcess.H257
-rw-r--r--Source/Particles/ParticleCreation/Make.package6
-rw-r--r--Source/Particles/ParticleCreation/TransformParticle.H44
-rw-r--r--Source/Particles/PhotonParticleContainer.cpp2
-rw-r--r--Source/Particles/PhysicalParticleContainer.cpp71
-rw-r--r--Source/Particles/RigidInjectedParticleContainer.cpp12
-rw-r--r--Source/Particles/WarpXParticleContainer.H5
12 files changed, 550 insertions, 265 deletions
diff --git a/Source/Particles/Deposition/CurrentDeposition.H b/Source/Particles/Deposition/CurrentDeposition.H
index 6da0f1155..2737eb008 100644
--- a/Source/Particles/Deposition/CurrentDeposition.H
+++ b/Source/Particles/Deposition/CurrentDeposition.H
@@ -4,6 +4,9 @@
#include "ShapeFactors.H"
#include <WarpX_Complex.H>
+#include <AMReX_Array4.H>
+#include <AMReX_REAL.H>
+
/* \brief Current Deposition for thread thread_num
* /param xp, yp, zp : Pointer to arrays of particle positions.
* \param wp : Pointer to array of particle weights.
@@ -208,69 +211,71 @@ void doEsirkepovDepositionShapeN (const amrex::ParticleReal * const xp,
const amrex::Real q,
const long n_rz_azimuthal_modes)
{
+ using namespace amrex;
+
// Whether ion_lev is a null pointer (do_ionization=0) or a real pointer
// (do_ionization=1)
- const bool do_ionization = ion_lev;
- const amrex::Real dxi = 1.0/dx[0];
- const amrex::Real dtsdx0 = dt*dxi;
- const amrex::Real xmin = xyzmin[0];
+ bool const do_ionization = ion_lev;
+ Real const dxi = 1.0_rt / dx[0];
+ Real const dtsdx0 = dt*dxi;
+ Real const xmin = xyzmin[0];
#if (defined WARPX_DIM_3D)
- const amrex::Real dyi = 1.0/dx[1];
- const amrex::Real dtsdy0 = dt*dyi;
- const amrex::Real ymin = xyzmin[1];
+ Real const dyi = 1.0_rt / dx[1];
+ Real const dtsdy0 = dt*dyi;
+ Real const ymin = xyzmin[1];
#endif
- const amrex::Real dzi = 1.0/dx[2];
- const amrex::Real dtsdz0 = dt*dzi;
- const amrex::Real zmin = xyzmin[2];
+ Real const dzi = 1.0_rt / dx[2];
+ Real const dtsdz0 = dt*dzi;
+ Real const zmin = xyzmin[2];
#if (defined WARPX_DIM_3D)
- const amrex::Real invdtdx = 1.0/(dt*dx[1]*dx[2]);
- const amrex::Real invdtdy = 1.0/(dt*dx[0]*dx[2]);
- const amrex::Real invdtdz = 1.0/(dt*dx[0]*dx[1]);
+ Real const invdtdx = 1.0_rt / (dt*dx[1]*dx[2]);
+ Real const invdtdy = 1.0_rt / (dt*dx[0]*dx[2]);
+ Real const invdtdz = 1.0_rt / (dt*dx[0]*dx[1]);
#elif (defined WARPX_DIM_XZ) || (defined WARPX_DIM_RZ)
- const amrex::Real invdtdx = 1.0/(dt*dx[2]);
- const amrex::Real invdtdz = 1.0/(dt*dx[0]);
- const amrex::Real invvol = 1.0/(dx[0]*dx[2]);
+ Real const invdtdx = 1.0_rt / (dt*dx[2]);
+ Real const invdtdz = 1.0_rt / (dt*dx[0]);
+ Real const invvol = 1.0_rt / (dx[0]*dx[2]);
#endif
#if (defined WARPX_DIM_RZ)
- const Complex I = Complex{0., 1.};
+ Complex const I = Complex{0., 1.};
#endif
- const amrex::Real clightsq = 1.0/PhysConst::c/PhysConst::c;
+ Real const clightsq = 1.0_rt / ( PhysConst::c * PhysConst::c );
// Loop over particles and deposit into Jx_arr, Jy_arr and Jz_arr
amrex::ParallelFor(
np_to_depose,
- [=] AMREX_GPU_DEVICE (long ip) {
+ [=] AMREX_GPU_DEVICE (long const ip) {
// --- Get particle quantities
- const amrex::Real gaminv = 1.0/std::sqrt(1.0 + uxp[ip]*uxp[ip]*clightsq
+ Real const gaminv = 1.0/std::sqrt(1.0 + uxp[ip]*uxp[ip]*clightsq
+ uyp[ip]*uyp[ip]*clightsq
+ uzp[ip]*uzp[ip]*clightsq);
// wqx, wqy wqz are particle current in each direction
- amrex::Real wq = q*wp[ip];
+ Real wq = q*wp[ip];
if (do_ionization){
wq *= ion_lev[ip];
}
- const amrex::Real wqx = wq*invdtdx;
+ Real const wqx = wq*invdtdx;
#if (defined WARPX_DIM_3D)
- const amrex::Real wqy = wq*invdtdy;
+ Real const wqy = wq*invdtdy;
#endif
- const amrex::Real wqz = wq*invdtdz;
+ Real const wqz = wq*invdtdz;
// computes current and old position in grid units
#if (defined WARPX_DIM_RZ)
- const amrex::Real xp_mid = xp[ip] - 0.5*dt*uxp[ip]*gaminv;
- const amrex::Real yp_mid = yp[ip] - 0.5*dt*uyp[ip]*gaminv;
- const amrex::Real xp_old = xp[ip] - dt*uxp[ip]*gaminv;
- const amrex::Real yp_old = yp[ip] - dt*uyp[ip]*gaminv;
- const amrex::Real rp_new = std::sqrt(xp[ip]*xp[ip] + yp[ip]*yp[ip]);
- const amrex::Real rp_mid = std::sqrt(xp_mid*xp_mid + yp_mid*yp_mid);
- const amrex::Real rp_old = std::sqrt(xp_old*xp_old + yp_old*yp_old);
- amrex::Real costheta_new, sintheta_new;
- if (rp_new > 0.) {
+ Real const xp_mid = xp[ip] - 0.5_rt * dt*uxp[ip]*gaminv;
+ Real const yp_mid = yp[ip] - 0.5_rt * dt*uyp[ip]*gaminv;
+ Real const xp_old = xp[ip] - dt*uxp[ip]*gaminv;
+ Real const yp_old = yp[ip] - dt*uyp[ip]*gaminv;
+ Real const rp_new = std::sqrt(xp[ip]*xp[ip] + yp[ip]*yp[ip]);
+ Real const rp_mid = std::sqrt(xp_mid*xp_mid + yp_mid*yp_mid);
+ Real const rp_old = std::sqrt(xp_old*xp_old + yp_old*yp_old);
+ Real costheta_new, sintheta_new;
+ if (rp_new > 0._rt) {
costheta_new = xp[ip]/rp_new;
sintheta_new = yp[ip]/rp_new;
} else {
@@ -278,7 +283,7 @@ void doEsirkepovDepositionShapeN (const amrex::ParticleReal * const xp,
sintheta_new = 0.;
}
amrex::Real costheta_mid, sintheta_mid;
- if (rp_mid > 0.) {
+ if (rp_mid > 0._rt) {
costheta_mid = xp_mid/rp_mid;
sintheta_mid = yp_mid/rp_mid;
} else {
@@ -286,7 +291,7 @@ void doEsirkepovDepositionShapeN (const amrex::ParticleReal * const xp,
sintheta_mid = 0.;
}
amrex::Real costheta_old, sintheta_old;
- if (rp_old > 0.) {
+ if (rp_old > 0._rt) {
costheta_old = xp_old/rp_old;
sintheta_old = yp_old/rp_old;
} else {
@@ -296,37 +301,37 @@ void doEsirkepovDepositionShapeN (const amrex::ParticleReal * const xp,
const Complex xy_new0 = Complex{costheta_new, sintheta_new};
const Complex xy_mid0 = Complex{costheta_mid, sintheta_mid};
const Complex xy_old0 = Complex{costheta_old, sintheta_old};
- const amrex::Real x_new = (rp_new - xmin)*dxi;
- const amrex::Real x_old = (rp_old - xmin)*dxi;
+ Real const x_new = (rp_new - xmin)*dxi;
+ Real const x_old = (rp_old - xmin)*dxi;
#else
- const amrex::Real x_new = (xp[ip] - xmin)*dxi;
- const amrex::Real x_old = x_new - dtsdx0*uxp[ip]*gaminv;
+ Real const x_new = (xp[ip] - xmin)*dxi;
+ Real const x_old = x_new - dtsdx0*uxp[ip]*gaminv;
#endif
#if (defined WARPX_DIM_3D)
- const amrex::Real y_new = (yp[ip] - ymin)*dyi;
- const amrex::Real y_old = y_new - dtsdy0*uyp[ip]*gaminv;
+ Real const y_new = (yp[ip] - ymin)*dyi;
+ Real const y_old = y_new - dtsdy0*uyp[ip]*gaminv;
#endif
- const amrex::Real z_new = (zp[ip] - zmin)*dzi;
- const amrex::Real z_old = z_new - dtsdz0*uzp[ip]*gaminv;
+ Real const z_new = (zp[ip] - zmin)*dzi;
+ Real const z_old = z_new - dtsdz0*uzp[ip]*gaminv;
#if (defined WARPX_DIM_RZ)
- const amrex::Real vy = (-uxp[ip]*sintheta_mid + uyp[ip]*costheta_mid)*gaminv;
+ Real const vy = (-uxp[ip]*sintheta_mid + uyp[ip]*costheta_mid)*gaminv;
#elif (defined WARPX_DIM_XZ)
- const amrex::Real vy = uyp[ip]*gaminv;
+ Real const vy = uyp[ip]*gaminv;
#endif
// Shape factor arrays
// Note that there are extra values above and below
// to possibly hold the factor for the old particle
// which can be at a different grid location.
- amrex::Real sx_new[depos_order + 3] = {0.};
- amrex::Real sx_old[depos_order + 3] = {0.};
+ Real sx_new[depos_order + 3] = {0.};
+ Real sx_old[depos_order + 3] = {0.};
#if (defined WARPX_DIM_3D)
- amrex::Real sy_new[depos_order + 3] = {0.};
- amrex::Real sy_old[depos_order + 3] = {0.};
+ Real sy_new[depos_order + 3] = {0.};
+ Real sy_old[depos_order + 3] = {0.};
#endif
- amrex::Real sz_new[depos_order + 3] = {0.};
- amrex::Real sz_old[depos_order + 3] = {0.};
+ Real sz_new[depos_order + 3] = {0.};
+ Real sz_old[depos_order + 3] = {0.};
// --- Compute shape factors
// Compute shape factors for position as they are now and at old positions
@@ -397,7 +402,7 @@ void doEsirkepovDepositionShapeN (const amrex::ParticleReal * const xp,
Complex xy_mid = xy_mid0; // Throughout the following loop, xy_mid takes the value e^{i m theta}
for (int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) {
// The factor 2 comes from the normalization of the modes
- const Complex djr_cmplx = amrex::Real(2.)*sdxi*xy_mid;
+ const Complex djr_cmplx = 2._rt *sdxi*xy_mid;
amrex::Gpu::Atomic::Add( &Jx_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode-1), djr_cmplx.real());
amrex::Gpu::Atomic::Add( &Jx_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode), djr_cmplx.imag());
xy_mid = xy_mid*xy_mid0;
@@ -407,8 +412,8 @@ void doEsirkepovDepositionShapeN (const amrex::ParticleReal * const xp,
}
for (int k=dkl; k<=depos_order+2-dku; k++) {
for (int i=dil; i<=depos_order+2-diu; i++) {
- const amrex::Real sdyj = wq*vy*invvol*((sz_new[k] + 0.5*(sz_old[k] - sz_new[k]))*sx_new[i] +
- (0.5*sz_new[k] + 1./3.*(sz_old[k] - sz_new[k]))*(sx_old[i] - sx_new[i]));
+ Real const sdyj = wq*vy*invvol*((sz_new[k] + 0.5_rt * (sz_old[k] - sz_new[k]))*sx_new[i] +
+ (0.5_rt * sz_new[k] + 1._rt / 3._rt *(sz_old[k] - sz_new[k]))*(sx_old[i] - sx_new[i]));
amrex::Gpu::Atomic::Add( &Jy_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 0), sdyj);
#if (defined WARPX_DIM_RZ)
Complex xy_new = xy_new0;
@@ -418,7 +423,7 @@ void doEsirkepovDepositionShapeN (const amrex::ParticleReal * const xp,
for (int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) {
// The factor 2 comes from the normalization of the modes
// The minus sign comes from the different convention with respect to Davidson et al.
- const Complex djt_cmplx = -amrex::Real(2.)*I*(i_new-1 + i + xmin*dxi)*wq*invdtdx/(amrex::Real)imode*
+ const Complex djt_cmplx = -2._rt * I*(i_new-1 + i + xmin*dxi)*wq*invdtdx/(amrex::Real)imode*
(sx_new[i]*sz_new[k]*(xy_new - xy_mid) + sx_old[i]*sz_old[k]*(xy_mid - xy_old));
amrex::Gpu::Atomic::Add( &Jy_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode-1), djt_cmplx.real());
amrex::Gpu::Atomic::Add( &Jy_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode), djt_cmplx.imag());
@@ -430,15 +435,15 @@ void doEsirkepovDepositionShapeN (const amrex::ParticleReal * const xp,
}
}
for (int i=dil; i<=depos_order+2-diu; i++) {
- amrex::Real sdzk = 0.;
+ Real sdzk = 0.;
for (int k=dkl; k<=depos_order+1-dku; k++) {
- sdzk += wqz*(sz_old[k] - sz_new[k])*(sx_new[i] + 0.5*(sx_old[i] - sx_new[i]));
+ sdzk += wqz*(sz_old[k] - sz_new[k])*(sx_new[i] + 0.5_rt * (sx_old[i] - sx_new[i]));
amrex::Gpu::Atomic::Add( &Jz_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 0), sdzk);
#if (defined WARPX_DIM_RZ)
Complex xy_mid = xy_mid0; // Throughout the following loop, xy_mid takes the value e^{i m theta}
for (int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) {
// The factor 2 comes from the normalization of the modes
- const Complex djz_cmplx = amrex::Real(2.)*sdzk*xy_mid;
+ const Complex djz_cmplx = 2._rt * sdzk * xy_mid;
amrex::Gpu::Atomic::Add( &Jz_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode-1), djz_cmplx.real());
amrex::Gpu::Atomic::Add( &Jz_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode), djz_cmplx.imag());
xy_mid = xy_mid*xy_mid0;
diff --git a/Source/Particles/Make.package b/Source/Particles/Make.package
index c9d520873..a6c124ddd 100644
--- a/Source/Particles/Make.package
+++ b/Source/Particles/Make.package
@@ -19,6 +19,7 @@ include $(WARPX_HOME)/Source/Particles/Pusher/Make.package
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_LOCATIONS += $(WARPX_HOME)/Source/Particles
VPATH_LOCATIONS += $(WARPX_HOME)/Source/Particles
diff --git a/Source/Particles/MultiParticleContainer.H b/Source/Particles/MultiParticleContainer.H
index 949173052..3341a8fe8 100644
--- a/Source/Particles/MultiParticleContainer.H
+++ b/Source/Particles/MultiParticleContainer.H
@@ -1,14 +1,15 @@
-
#ifndef WARPX_ParticleContainer_H_
#define WARPX_ParticleContainer_H_
-#include <AMReX_Particles.H>
+#include "ElementaryProcess.H"
+
#include <WarpXParticleContainer.H>
#include <PhysicalParticleContainer.H>
#include <RigidInjectedParticleContainer.H>
#include <PhotonParticleContainer.H>
#include <LaserParticleContainer.H>
+#include <AMReX_Particles.H>
#ifdef WARPX_QED
#include <BreitWheelerEngineWrapper.H>
#include <QuantumSyncEngineWrapper.H>
@@ -162,9 +163,9 @@ public:
int nSpecies() const {return nspecies;}
- int nSpeciesBoostedFrameDiags() const {return nspecies_boosted_frame_diags;}
- int mapSpeciesBoostedFrameDiags(int i) const {return map_species_boosted_frame_diags[i];}
- int doBoostedFrameDiags() const {return do_boosted_frame_diags;}
+ int nSpeciesBackTransformedDiagnostics() const {return nspecies_back_transformed_diagnostics;}
+ int mapSpeciesBackTransformedDiagnostics(int i) const {return map_species_back_transformed_diagnostics[i];}
+ int doBackTransformedDiagnostics() const {return do_back_transformed_diagnostics;}
int nSpeciesDepositOnMainGrid () const {
bool const onMainGrid = true;
@@ -196,6 +197,8 @@ public:
PhysicalParticleContainer& GetPCtmp () { return *pc_tmp; }
+ IonizationProcess ionization_process;
+
protected:
// Particle container types
@@ -236,12 +239,12 @@ private:
void mapSpeciesProduct ();
int getSpeciesID (std::string product_str);
- // Number of species dumped in BoostedFrameDiagnostics
- int nspecies_boosted_frame_diags = 0;
- // 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;
+ // Number of species dumped in BackTransformedDiagnostics
+ int nspecies_back_transformed_diagnostics = 0;
+ // map_species_back_transformed_diagnostics[i] is the species ID in
+ // MultiParticleContainer for 0<i<nspecies_back_transformed_diagnostics
+ std::vector<int> map_species_back_transformed_diagnostics;
+ int do_back_transformed_diagnostics = 0;
// runtime parameters
int nlasers = 0;
diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp
index f4c00404b..de2b4583d 100644
--- a/Source/Particles/MultiParticleContainer.cpp
+++ b/Source/Particles/MultiParticleContainer.cpp
@@ -38,16 +38,17 @@ MultiParticleContainer::MultiParticleContainer (AmrCore* amr_core)
// Compute the number of species for which lab-frame data is dumped
// nspecies_lab_frame_diags, and map their ID to MultiParticleContainer
// particle IDs in map_species_lab_diags.
- map_species_boosted_frame_diags.resize(nspecies);
- nspecies_boosted_frame_diags = 0;
+ map_species_back_transformed_diagnostics.resize(nspecies);
+ nspecies_back_transformed_diagnostics = 0;
for (int i=0; i<nspecies; i++){
auto& pc = allcontainers[i];
- if (pc->do_boosted_frame_diags){
- map_species_boosted_frame_diags[nspecies_boosted_frame_diags] = i;
- do_boosted_frame_diags = 1;
- nspecies_boosted_frame_diags += 1;
+ if (pc->do_back_transformed_diagnostics){
+ map_species_back_transformed_diagnostics[nspecies_back_transformed_diagnostics] = i;
+ do_back_transformed_diagnostics = 1;
+ nspecies_back_transformed_diagnostics += 1;
}
}
+ ionization_process = IonizationProcess();
}
void
@@ -315,7 +316,11 @@ void
MultiParticleContainer::Redistribute ()
{
for (auto& pc : allcontainers) {
- pc->Redistribute();
+ if ( (pc->NumRuntimeRealComps()>0) || (pc->NumRuntimeIntComps()>0) ) {
+ pc->RedistributeCPU();
+ } else {
+ pc->Redistribute();
+ }
}
}
@@ -323,7 +328,11 @@ void
MultiParticleContainer::RedistributeLocal (const int num_ghost)
{
for (auto& pc : allcontainers) {
- pc->Redistribute(0, 0, 0, num_ghost);
+ if ( (pc->NumRuntimeRealComps()>0) || (pc->NumRuntimeIntComps()>0) ) {
+ pc->RedistributeCPU(0, 0, 0, num_ghost);
+ } else {
+ pc->Redistribute(0, 0, 0, num_ghost);
+ }
}
}
@@ -387,8 +396,8 @@ 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];
+ for (int i = 0; i < nspecies_back_transformed_diagnostics; ++i){
+ int isp = map_species_back_transformed_diagnostics[i];
WarpXParticleContainer* pc = allcontainers[isp].get();
WarpXParticleContainer::DiagnosticParticles diagnostic_particles;
pc->GetParticleSlice(direction, z_old, z_new, t_boost, t_lab, dt, diagnostic_particles);
@@ -526,146 +535,6 @@ MultiParticleContainer::getSpeciesID (std::string product_str)
return i_product;
}
-namespace
-{
- // For particle i in mfi, if is_ionized[i]=1, copy particle
- // particle i from container pc_source into pc_product
- void createIonizedParticles (
- int lev, const MFIter& mfi,
- std::unique_ptr< WarpXParticleContainer>& pc_source,
- std::unique_ptr< WarpXParticleContainer>& pc_product,
- amrex::Gpu::ManagedDeviceVector<int>& is_ionized)
- {
- BL_PROFILE("createIonizedParticles");
-
- const int * const AMREX_RESTRICT p_is_ionized = is_ionized.dataPtr();
-
- 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();
- GpuArray<ParticleReal*,PIdx::nattribs> attribs_source;
- for (int ia = 0; ia < PIdx::nattribs; ++ia) {
- attribs_source[ia] = soa_source.GetRealData(ia).data();
- }
- // --- source runtime attribs
- GpuArray<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();
-
- // Indices of product particle for each ionized 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_ionized
- // Strictly speaking, i_product should be an exclusive_scan of
- // is_ionized. However, for indices where is_ionized 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_ionized
- // 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_ionized.begin(), is_ionized.end(), i_product.begin());
- int np_ionized = i_product[np_source-1];
- if (np_ionized == 0) return;
- int* AMREX_RESTRICT p_i_product = i_product.dataPtr();
-
- // Get product particle data
- auto& ptile_product = pc_product->GetParticles(lev)[std::make_pair(grid_id,tile_id)];
- // old and new (i.e., including ionized 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_ionized;
- // Allocate extra space in product species for ionized particles.
- ptile_product.resize(np_product_new);
- // --- product AoS particle data
- // First element is the first newly-created product particle
- WarpXParticleContainer::ParticleType* particles_product = ptile_product.GetArrayOfStructs()().data() + np_product_old;
- // --- product SoA particle data
- auto& soa_product = ptile_product.GetStructOfArrays();
- GpuArray<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
- GpuArray<ParticleReal*,6> runtime_attribs_product;
- bool do_boosted_product = WarpX::do_boosted_frame_diagnostic
- && pc_product->DoBoostedFrameDiags();
- if (do_boosted_product) {
- 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;
- }
-
- int pid_product;
-#pragma omp critical (doFieldIonization_nextid)
- {
- // 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_ionized);
- }
- const int cpuid = ParallelDescriptor::MyProc();
-
- // Loop over all source particles. If is_ionized, copy particle data
- // to corresponding product particle.
- amrex::For(
- np_source, [=] AMREX_GPU_DEVICE (int is) noexcept
- {
- if(p_is_ionized[is]){
- // offset of 1 due to inclusive scan
- int ip = p_i_product[is]-1;
- // is: index of ionized particle in source species
- // ip: index of corresponding new particle in product species
- WarpXParticleContainer::ParticleType& p_product = particles_product[ip];
- WarpXParticleContainer::ParticleType& p_source = particles_source[is];
- // Copy particle from source to product: AoS
- p_product.id() = pid_product + ip;
- p_product.cpu() = 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) {
- attribs_product[ia][ip] = 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 (do_boosted_product) {
- runtime_attribs_product[0][ip] = p_source.pos(0);
- runtime_attribs_product[1][ip] = p_source.pos(1);
- runtime_attribs_product[2][ip] = p_source.pos(2);
- runtime_attribs_product[3][ip] = runtime_uold_source[0][ip];
- runtime_attribs_product[4][ip] = runtime_uold_source[1][ip];
- runtime_attribs_product[5][ip] = runtime_uold_source[2][ip];
- }
- }
- }
- );
- }
-}
-
void
MultiParticleContainer::doFieldIonization ()
{
@@ -727,7 +596,17 @@ MultiParticleContainer::doFieldIonization ()
amrex::Gpu::ManagedDeviceVector<int> is_ionized;
pc_source->buildIonizationMask(mfi, lev, is_ionized);
// Create particles in pc_product
- createIonizedParticles(lev, mfi, pc_source, pc_product, is_ionized);
+ 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();
}
} // lev
} // pc_source
diff --git a/Source/Particles/ParticleCreation/CopyParticle.H b/Source/Particles/ParticleCreation/CopyParticle.H
new file mode 100644
index 000000000..bd2d530fb
--- /dev/null
+++ b/Source/Particles/ParticleCreation/CopyParticle.H
@@ -0,0 +1,90 @@
+#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_source 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/ElementaryProcess.H b/Source/Particles/ParticleCreation/ElementaryProcess.H
new file mode 100644
index 000000000..fa791c244
--- /dev/null
+++ b/Source/Particles/ParticleCreation/ElementaryProcess.H
@@ -0,0 +1,257 @@
+#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::Vector<copyParticle> v_copy_functor;
+ 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;
+#pragma omp critical (doParticleCreation_nextid)
+ {
+ // 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::Vector<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 p_copy_functor = v_copy_functor.data();
+ 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/Make.package b/Source/Particles/ParticleCreation/Make.package
new file mode 100644
index 000000000..6e32f4a77
--- /dev/null
+++ b/Source/Particles/ParticleCreation/Make.package
@@ -0,0 +1,6 @@
+CEXE_headers += ElementaryProcess.H
+CEXE_headers += CopyParticle.H
+CEXE_headers += TransformParticle.H
+
+INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Particles/ParticleCreation/
+VPATH_LOCATIONS += $(WARPX_HOME)/Source/Particles/ParticleCreation/
diff --git a/Source/Particles/ParticleCreation/TransformParticle.H b/Source/Particles/ParticleCreation/TransformParticle.H
new file mode 100644
index 000000000..14e534d0e
--- /dev/null
+++ b/Source/Particles/ParticleCreation/TransformParticle.H
@@ -0,0 +1,44 @@
+#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 particle particle struct
+ */
+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/PhotonParticleContainer.cpp b/Source/Particles/PhotonParticleContainer.cpp
index cf2ffbec4..612da01ca 100644
--- a/Source/Particles/PhotonParticleContainer.cpp
+++ b/Source/Particles/PhotonParticleContainer.cpp
@@ -73,7 +73,7 @@ PhotonParticleContainer::PushPX(WarpXParIter& pti,
const ParticleReal* const AMREX_RESTRICT By = attribs[PIdx::By].dataPtr();
const ParticleReal* const AMREX_RESTRICT Bz = attribs[PIdx::Bz].dataPtr();
- if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags)
+ if (WarpX::do_back_transformed_diagnostics && do_back_transformed_diagnostics)
{
copy_attribs(pti, x, y, z);
}
diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp
index d12a4dbff..fed7266e1 100644
--- a/Source/Particles/PhysicalParticleContainer.cpp
+++ b/Source/Particles/PhysicalParticleContainer.cpp
@@ -43,7 +43,7 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp
pp.query("do_continuous_injection", do_continuous_injection);
// Whether to plot back-transformed (lab-frame) diagnostics
// for this species.
- pp.query("do_boosted_frame_diags", do_boosted_frame_diags);
+ pp.query("do_back_transformed_diagnostics", do_back_transformed_diagnostics);
pp.query("do_field_ionization", do_field_ionization);
@@ -62,15 +62,16 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp
"Radiation reaction can be enabled only if Boris pusher is used");
//_____________________________
-
#ifdef AMREX_USE_GPU
- AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
- do_field_ionization == 0,
- "Field ionization does not work on GPU so far, because the current "
- "version of Redistribute in AMReX does not work with runtime parameters");
+ Print()<<"\n-----------------------------------------------------\n";
+ Print()<<"WARNING: field ionization on GPU uses RedistributeCPU\n";
+ Print()<<"-----------------------------------------------------\n\n";
+ //AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
+ //do_field_ionization == 0,
+ //"Field ionization does not work on GPU so far, because the current "
+ //"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);
@@ -86,7 +87,7 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp
//variable to set plot_flags size
int plot_flag_size = PIdx::nattribs;
- if(WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags)
+ if(WarpX::do_back_transformed_diagnostics && do_back_transformed_diagnostics)
plot_flag_size += 6;
#ifdef WARPX_QED
@@ -441,13 +442,13 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox)
for (int dir=0; dir<AMREX_SPACEDIM; dir++) {
if ( tile_realbox.lo(dir) <= part_realbox.hi(dir) ) {
Real ncells_adjust = std::floor( (tile_realbox.lo(dir) - part_realbox.lo(dir))/dx[dir] );
- overlap_realbox.setLo( dir, part_realbox.lo(dir) + std::max(ncells_adjust, Real(0.)) * dx[dir]);
+ overlap_realbox.setLo( dir, part_realbox.lo(dir) + std::max(ncells_adjust, 0._rt) * dx[dir]);
} else {
no_overlap = true; break;
}
if ( tile_realbox.hi(dir) >= part_realbox.lo(dir) ) {
Real ncells_adjust = std::floor( (part_realbox.hi(dir) - tile_realbox.hi(dir))/dx[dir] );
- overlap_realbox.setHi( dir, part_realbox.hi(dir) - std::max(ncells_adjust, Real(0.)) * dx[dir]);
+ overlap_realbox.setHi( dir, part_realbox.hi(dir) - std::max(ncells_adjust, 0._rt) * dx[dir]);
} else {
no_overlap = true; break;
}
@@ -1012,12 +1013,12 @@ PhysicalParticleContainer::FieldGather (int lev,
const FArrayBox& byfab = By[pti];
const FArrayBox& bzfab = Bz[pti];
- 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]);
+ Exp.assign(np,WarpX::E_external_particle[0]);
+ Eyp.assign(np,WarpX::E_external_particle[1]);
+ Ezp.assign(np,WarpX::E_external_particle[2]);
+ Bxp.assign(np,WarpX::B_external_particle[0]);
+ Byp.assign(np,WarpX::B_external_particle[1]);
+ Bzp.assign(np,WarpX::B_external_particle[2]);
//
// copy data from particle container to temp arrays
@@ -1078,7 +1079,7 @@ PhysicalParticleContainer::Evolve (int lev,
bool has_buffer = cEx || cjx;
- if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags)
+ if (WarpX::do_back_transformed_diagnostics && do_back_transformed_diagnostics)
{
for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti)
{
@@ -1148,13 +1149,13 @@ PhysicalParticleContainer::Evolve (int lev,
exfab, eyfab, ezfab, bxfab, byfab, bzfab);
}
- Exp.assign(np,WarpX::E_external[0]);
- Eyp.assign(np,WarpX::E_external[1]);
- Ezp.assign(np,WarpX::E_external[2]);
+ Exp.assign(np,WarpX::E_external_particle[0]);
+ Eyp.assign(np,WarpX::E_external_particle[1]);
+ Ezp.assign(np,WarpX::E_external_particle[2]);
- Bxp.assign(np,WarpX::B_external[0]);
- Byp.assign(np,WarpX::B_external[1]);
- Bzp.assign(np,WarpX::B_external[2]);
+ Bxp.assign(np,WarpX::B_external_particle[0]);
+ Byp.assign(np,WarpX::B_external_particle[1]);
+ Bzp.assign(np,WarpX::B_external_particle[2]);
// Determine which particles deposit/gather in the buffer, and
// which particles deposit/gather in the fine patch
@@ -1427,9 +1428,9 @@ PhysicalParticleContainer::SplitParticles(int lev)
// before splitting results in a uniform distribution after splitting
const amrex::Vector<int> ppc_nd = plasma_injector->num_particles_per_cell_each_dim;
const std::array<Real,3>& dx = WarpX::CellSize(lev);
- amrex::Vector<amrex::Real> split_offset = {dx[0]/2./ppc_nd[0],
- dx[1]/2./ppc_nd[1],
- dx[2]/2./ppc_nd[2]};
+ amrex::Vector<Real> split_offset = {dx[0]/2._rt/ppc_nd[0],
+ dx[1]/2._rt/ppc_nd[1],
+ dx[2]/2._rt/ppc_nd[2]};
// particle Array Of Structs data
auto& particles = pti.GetArrayOfStructs();
@@ -1585,7 +1586,7 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti,
const ParticleReal* const AMREX_RESTRICT By = attribs[PIdx::By].dataPtr();
const ParticleReal* const AMREX_RESTRICT Bz = attribs[PIdx::Bz].dataPtr();
- if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags && (a_dt_type!=DtType::SecondHalf))
+ if (WarpX::do_back_transformed_diagnostics && do_back_transformed_diagnostics && (a_dt_type!=DtType::SecondHalf))
{
copy_attribs(pti, x, y, z);
}
@@ -1701,13 +1702,13 @@ PhysicalParticleContainer::PushP (int lev, Real dt,
const FArrayBox& byfab = By[pti];
const FArrayBox& bzfab = Bz[pti];
- Exp.assign(np,WarpX::E_external[0]);
- Eyp.assign(np,WarpX::E_external[1]);
- Ezp.assign(np,WarpX::E_external[2]);
+ Exp.assign(np,WarpX::E_external_particle[0]);
+ Eyp.assign(np,WarpX::E_external_particle[1]);
+ Ezp.assign(np,WarpX::E_external_particle[2]);
- Bxp.assign(np,WarpX::B_external[0]);
- Byp.assign(np,WarpX::B_external[1]);
- Bzp.assign(np,WarpX::B_external[2]);
+ Bxp.assign(np,WarpX::B_external_particle[0]);
+ Byp.assign(np,WarpX::B_external_particle[1]);
+ Bzp.assign(np,WarpX::B_external_particle[2]);
//
// copy data from particle container to temp arrays
@@ -1842,7 +1843,7 @@ void PhysicalParticleContainer::GetParticleSlice(const int direction, const Real
// Note the the slice should always move in the negative boost direction.
AMREX_ALWAYS_ASSERT(z_new < z_old);
- AMREX_ALWAYS_ASSERT(do_boosted_frame_diags == 1);
+ AMREX_ALWAYS_ASSERT(do_back_transformed_diagnostics == 1);
const int nlevs = std::max(0, finestLevel()+1);
@@ -2218,8 +2219,6 @@ PhysicalParticleContainer::buildIonizationMask (const amrex::MFIter& mfi, const
Real p = 1. - std::exp( - w_dtau );
if (random_draw < p){
- // increment particle's ionization level
- ion_lev[i] += 1;
// update mask
p_ionization_mask[i] = 1;
}
diff --git a/Source/Particles/RigidInjectedParticleContainer.cpp b/Source/Particles/RigidInjectedParticleContainer.cpp
index 535ffec6f..507206041 100644
--- a/Source/Particles/RigidInjectedParticleContainer.cpp
+++ b/Source/Particles/RigidInjectedParticleContainer.cpp
@@ -392,12 +392,12 @@ RigidInjectedParticleContainer::PushP (int lev, Real dt,
const FArrayBox& byfab = By[pti];
const FArrayBox& bzfab = Bz[pti];
- 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]);
+ Exp.assign(np,WarpX::E_external_particle[0]);
+ Eyp.assign(np,WarpX::E_external_particle[1]);
+ Ezp.assign(np,WarpX::E_external_particle[2]);
+ Bxp.assign(np,WarpX::B_external_particle[0]);
+ Byp.assign(np,WarpX::B_external_particle[1]);
+ Bzp.assign(np,WarpX::B_external_particle[2]);
//
// copy data from particle container to temp arrays
diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H
index 879f4782e..b41dcede3 100644
--- a/Source/Particles/WarpXParticleContainer.H
+++ b/Source/Particles/WarpXParticleContainer.H
@@ -273,13 +273,14 @@ public:
AddIntComp(comm);
}
- int DoBoostedFrameDiags () const { return do_boosted_frame_diags; }
+ 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;}
protected:
@@ -316,7 +317,7 @@ protected:
amrex::Gpu::ManagedVector<amrex::Real> adk_exp_prefactor;
std::string physical_element;
- int do_boosted_frame_diags = 1;
+ int do_back_transformed_diagnostics = 1;
#ifdef WARPX_QED
bool do_qed = false;
generated by cgit v1.2.3 (git 2.46.0) at 2025-08-14 10:25:47 +0000