diff options
Diffstat (limited to 'Source/Particles')
| -rwxr-xr-x | Source/Particles/Deposition/ChargeDeposition.H | 13 | ||||
| -rw-r--r-- | Source/Particles/Deposition/CurrentDeposition.H | 26 | ||||
| -rw-r--r-- | Source/Particles/MultiParticleContainer.H | 6 | ||||
| -rw-r--r-- | Source/Particles/MultiParticleContainer.cpp | 256 | ||||
| -rw-r--r-- | Source/Particles/PhysicalParticleContainer.H | 5 | ||||
| -rw-r--r-- | Source/Particles/PhysicalParticleContainer.cpp | 245 | ||||
| -rw-r--r-- | Source/Particles/WarpXParticleContainer.H | 28 | ||||
| -rw-r--r-- | Source/Particles/WarpXParticleContainer.cpp | 125 |
8 files changed, 632 insertions, 72 deletions
diff --git a/Source/Particles/Deposition/ChargeDeposition.H b/Source/Particles/Deposition/ChargeDeposition.H index 26d42370e..4253f5e76 100755 --- a/Source/Particles/Deposition/ChargeDeposition.H +++ b/Source/Particles/Deposition/ChargeDeposition.H @@ -6,6 +6,10 @@ /* \brief Charge Deposition for thread thread_num * /param xp, yp, zp : Pointer to arrays of particle positions. * \param wp : Pointer to array of particle weights. + * \param ion_lev : Pointer to array of particle ionization level. This is + required to have the charge of each macroparticle + since q is a scalar. For non-ionizable species, + ion_lev is a null pointer. * \param rho_arr : Array4 of charge density, either full array or tile. * \param np_to_depose : Number of particles for which current is deposited. * \param dx : 3D cell size @@ -18,6 +22,7 @@ void doChargeDepositionShapeN(const amrex::Real * const xp, const amrex::Real * const yp, const amrex::Real * const zp, const amrex::Real * const wp, + const int * const ion_lev, const amrex::Array4<amrex::Real>& rho_arr, const long np_to_depose, const std::array<amrex::Real,3>& dx, @@ -25,6 +30,9 @@ void doChargeDepositionShapeN(const amrex::Real * const xp, const amrex::Dim3 lo, const amrex::Real q) { + // 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 dzi = 1.0/dx[2]; #if (AMREX_SPACEDIM == 2) @@ -43,7 +51,10 @@ void doChargeDepositionShapeN(const amrex::Real * const xp, np_to_depose, [=] AMREX_GPU_DEVICE (long ip) { // --- Get particle quantities - const amrex::Real wq = q*wp[ip]*invvol; + amrex::Real wq = q*wp[ip]*invvol; + if (do_ionization){ + wq *= ion_lev[ip]; + } // --- Compute shape factors // x direction diff --git a/Source/Particles/Deposition/CurrentDeposition.H b/Source/Particles/Deposition/CurrentDeposition.H index f1e30e1be..e8e7fab0a 100644 --- a/Source/Particles/Deposition/CurrentDeposition.H +++ b/Source/Particles/Deposition/CurrentDeposition.H @@ -7,6 +7,10 @@ * /param xp, yp, zp : Pointer to arrays of particle positions. * \param wp : Pointer to array of particle weights. * \param uxp uyp uzp : Pointer to arrays of particle momentum. + * \param ion_lev : Pointer to array of particle ionization level. This is + required to have the charge of each macroparticle + since q is a scalar. For non-ionizable species, + ion_lev is a null pointer. * \param jx_arr : Array4 of current density, either full array or tile. * \param jy_arr : Array4 of current density, either full array or tile. * \param jz_arr : Array4 of current density, either full array or tile. @@ -26,6 +30,7 @@ void doDepositionShapeN(const amrex::Real * const xp, const amrex::Real * const uxp, const amrex::Real * const uyp, const amrex::Real * const uzp, + const int * const ion_lev, const amrex::Array4<amrex::Real>& jx_arr, const amrex::Array4<amrex::Real>& jy_arr, const amrex::Array4<amrex::Real>& jz_arr, @@ -36,6 +41,9 @@ void doDepositionShapeN(const amrex::Real * const xp, const amrex::Real stagger_shift, const amrex::Real q) { + // 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 dzi = 1.0/dx[2]; const amrex::Real dts2dx = 0.5*dt*dxi; @@ -61,7 +69,10 @@ void doDepositionShapeN(const amrex::Real * const xp, const amrex::Real gaminv = 1.0/std::sqrt(1.0 + uxp[ip]*uxp[ip]*clightsq + uyp[ip]*uyp[ip]*clightsq + uzp[ip]*uzp[ip]*clightsq); - const amrex::Real wq = q*wp[ip]; + amrex::Real wq = q*wp[ip]; + if (do_ionization){ + wq *= ion_lev[ip]; + } const amrex::Real vx = uxp[ip]*gaminv; const amrex::Real vy = uyp[ip]*gaminv; const amrex::Real vz = uzp[ip]*gaminv; @@ -161,6 +172,10 @@ void doDepositionShapeN(const amrex::Real * const xp, * /param xp, yp, zp : Pointer to arrays of particle positions. * \param wp : Pointer to array of particle weights. * \param uxp uyp uzp : Pointer to arrays of particle momentum. + * \param ion_lev : Pointer to array of particle ionization level. This is + required to have the charge of each macroparticle + since q is a scalar. For non-ionizable species, + ion_lev is a null pointer. * \param Jx_arr : Array4 of current density, either full array or tile. * \param Jy_arr : Array4 of current density, either full array or tile. * \param Jz_arr : Array4 of current density, either full array or tile. @@ -179,6 +194,7 @@ void doEsirkepovDepositionShapeN (const amrex::Real * const xp, const amrex::Real * const uxp, const amrex::Real * const uyp, const amrex::Real * const uzp, + const int * ion_lev, const amrex::Array4<amrex::Real>& Jx_arr, const amrex::Array4<amrex::Real>& Jy_arr, const amrex::Array4<amrex::Real>& Jz_arr, @@ -189,6 +205,9 @@ void doEsirkepovDepositionShapeN (const amrex::Real * const xp, const amrex::Dim3 lo, const amrex::Real q) { + // 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]; @@ -224,7 +243,10 @@ void doEsirkepovDepositionShapeN (const amrex::Real * const xp, + uzp[ip]*uzp[ip]*clightsq); // wqx, wqy wqz are particle current in each direction - const amrex::Real wq = q*wp[ip]; + amrex::Real wq = q*wp[ip]; + if (do_ionization){ + wq *= ion_lev[ip]; + } const amrex::Real wqx = wq*invdtdx; #if (defined WARPX_DIM_3D) const amrex::Real wqy = wq*invdtdy; diff --git a/Source/Particles/MultiParticleContainer.H b/Source/Particles/MultiParticleContainer.H index 7c9ede411..aaed96423 100644 --- a/Source/Particles/MultiParticleContainer.H +++ b/Source/Particles/MultiParticleContainer.H @@ -7,7 +7,6 @@ #include <string> #include <AMReX_Particles.H> - #include <WarpXParticleContainer.H> #include <PhysicalParticleContainer.H> #include <RigidInjectedParticleContainer.H> @@ -128,6 +127,8 @@ public: /// std::unique_ptr<amrex::MultiFab> GetChargeDensity(int lev, bool local = false); + void doFieldIonization (); + void Checkpoint (const std::string& dir) const; void WritePlotFile (const std::string& dir) const; @@ -207,6 +208,9 @@ private: void ReadParameters (); + 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 diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 982e04e39..9e2cd097f 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -135,6 +135,9 @@ MultiParticleContainer::InitData () pc->InitData(); } pc_tmp->InitData(); + // For each species, get the ID of its product species. + // This is used for ionization and pair creation processes. + mapSpeciesProduct(); } @@ -450,7 +453,7 @@ MultiParticleContainer::UpdateContinuousInjectionPosition(Real dt) const } int -MultiParticleContainer::doContinuousInjection() const +MultiParticleContainer::doContinuousInjection () const { int warpx_do_continuous_injection = 0; for (int i=0; i<nspecies+nlasers; i++){ @@ -461,3 +464,254 @@ MultiParticleContainer::doContinuousInjection() const } return warpx_do_continuous_injection; } + +/* \brief Get ID of product species of each species. + * The users specifies the name of the product species, + * this routine get its ID. + */ +void +MultiParticleContainer::mapSpeciesProduct () +{ + for (int i=0; i<nspecies; i++){ + auto& pc = allcontainers[i]; + // If species pc has ionization on, find species with name + // pc->ionization_product_name and store its ID into + // pc->ionization_product. + if (pc->do_field_ionization){ + int i_product = getSpeciesID(pc->ionization_product_name); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + i != i_product, + "ERROR: ionization product cannot be the same species"); + pc->ionization_product = i_product; + } + } +} + +/* \brief Given a species name, return its ID. + */ +int +MultiParticleContainer::getSpeciesID (std::string product_str) +{ + int i_product; + bool found = 0; + // Loop over species + for (int i=0; i<nspecies; i++){ + // If species name matches, store its ID + // into i_product + if (species_names[i] == product_str){ + found = 1; + i_product = i; + } + } + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + found != 0, + "ERROR: could not find product species ID for ionization. Wrong name?"); + 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<Real*,PIdx::nattribs> attribs_source; + for (int ia = 0; ia < PIdx::nattribs; ++ia) { + attribs_source[ia] = soa_source.GetRealData(ia).data(); + } + // --- source runtime attribs + GpuArray<Real*,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<Real*,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<Real*,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 () +{ + BL_PROFILE("MPC::doFieldIonization"); + // Loop over all species. + // Ionized particles in pc_source create particles in pc_product + for (auto& pc_source : allcontainers){ + + // Skip if not ionizable + if (!pc_source->do_field_ionization){ continue; } + + // Get product species + auto& pc_product = allcontainers[pc_source->ionization_product]; + + for (int lev = 0; lev <= pc_source->finestLevel(); ++lev){ + + // When using runtime components, AMReX requires to touch all tiles + // in serial and create particles tiles with runtime components if + // they do not exist (or if they were defined by default, i.e., + // without runtime component). +#ifdef _OPENMP + // Touch all tiles of source species in serial if runtime attribs + for (MFIter mfi = pc_source->MakeMFIter(lev); mfi.isValid(); ++mfi) { + const int grid_id = mfi.index(); + const int tile_id = mfi.LocalTileIndex(); + pc_source->GetParticles(lev)[std::make_pair(grid_id,tile_id)]; + if ( (pc_source->NumRuntimeRealComps()>0) || (pc_source->NumRuntimeIntComps()>0) ) { + pc_source->DefineAndReturnParticleTile(lev, grid_id, tile_id); + } + } +#endif + // Touch all tiles of product species in serial + for (MFIter mfi = pc_source->MakeMFIter(lev); mfi.isValid(); ++mfi) { + const int grid_id = mfi.index(); + const int tile_id = mfi.LocalTileIndex(); + pc_product->GetParticles(lev)[std::make_pair(grid_id,tile_id)]; + pc_product->DefineAndReturnParticleTile(lev, grid_id, tile_id); + } + + // Enable tiling + MFItInfo info; + if (pc_source->do_tiling && Gpu::notInLaunchRegion()) { + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + pc_product->do_tiling, + "For ionization, either all or none of the " + "particle species must use tiling."); + info.EnableTiling(pc_source->tile_size); + } + +#ifdef _OPENMP + info.SetDynamic(true); +#pragma omp parallel +#endif + // Loop over all grids (if not tiling) or grids and tiles (if tiling) + for (MFIter mfi = pc_source->MakeMFIter(lev, info); mfi.isValid(); ++mfi) + { + // Ionization mask: one element per source particles. + // 0 if not ionized, 1 if ionized. + amrex::Gpu::ManagedDeviceVector<int> is_ionized; + pc_source->buildIonizationMask(mfi, lev, is_ionized); + // Create particles in pc_product + createIonizedParticles(lev, mfi, pc_source, pc_product, is_ionized); + } + } // lev + } // pc_source +} diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index b80619733..c7494fbdf 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -22,6 +22,8 @@ public: virtual void InitData () override; + void InitIonizationModule (); + #ifdef WARPX_DO_ELECTROSTATIC virtual void FieldGatherES(const amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>, 3> >& E, const amrex::Vector<std::unique_ptr<amrex::FabArray<amrex::BaseFab<int> > > >& masks) override; @@ -104,6 +106,9 @@ public: virtual void PostRestart () final {} void SplitParticles(int lev); + + virtual void buildIonizationMask (const amrex::MFIter& mfi, const int lev, + amrex::Gpu::ManagedDeviceVector<int>& ionization_mask) override; // Inject particles in Box 'part_box' virtual void AddParticles (int lev); diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index d10390204..1d1d6c1b7 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -6,6 +6,7 @@ #include <WarpX.H> #include <WarpXConst.H> #include <WarpXWrappers.h> +#include <IonizationEnergiesTable.H> #include <FieldGather.H> #include <WarpXAlgorithmSelection.H> @@ -37,6 +38,14 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp // for this species. pp.query("do_boosted_frame_diags", do_boosted_frame_diags); + pp.query("do_field_ionization", do_field_ionization); +#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"); +#endif + pp.query("plot_species", plot_species); int do_user_plot_vars; do_user_plot_vars = pp.queryarr("plot_vars", plot_vars); @@ -77,6 +86,9 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core) 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 } @@ -196,17 +208,20 @@ PhysicalParticleContainer::CheckAndAddParticle(Real x, Real y, Real z, attribs[PIdx::uz] = u[2]; attribs[PIdx::w ] = weight; - if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) + if ( (NumRuntimeRealComps()>0) || (NumRuntimeIntComps()>0) ) { // need to create old values auto& particle_tile = DefineAndReturnParticleTile(0, 0, 0); - particle_tile.push_back_real(particle_comps["xold"], x); - particle_tile.push_back_real(particle_comps["yold"], y); - particle_tile.push_back_real(particle_comps["zold"], z); - - particle_tile.push_back_real(particle_comps["uxold"], u[0]); - particle_tile.push_back_real(particle_comps["uyold"], u[1]); - particle_tile.push_back_real(particle_comps["uzold"], u[2]); + if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) + { + particle_tile.push_back_real(particle_comps["xold"], x); + particle_tile.push_back_real(particle_comps["yold"], y); + particle_tile.push_back_real(particle_comps["zold"], z); + + particle_tile.push_back_real(particle_comps["uxold"], u[0]); + particle_tile.push_back_real(particle_comps["uyold"], u[1]); + particle_tile.push_back_real(particle_comps["uzold"], u[2]); + } } // add particle AddOneParticle(0, 0, 0, x, y, z, attribs); @@ -289,7 +304,7 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) const int grid_id = mfi.index(); const int tile_id = mfi.LocalTileIndex(); GetParticles(lev)[std::make_pair(grid_id, tile_id)]; - if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { + if ( (NumRuntimeRealComps()>0) || (NumRuntimeIntComps()>0) ) { DefineAndReturnParticleTile(lev, grid_id, tile_id); } } @@ -416,10 +431,12 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) auto& particle_tile = GetParticles(lev)[std::make_pair(grid_id,tile_id)]; bool do_boosted = false; - if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { - do_boosted = true; + if ( (NumRuntimeRealComps()>0) || (NumRuntimeIntComps()>0) ) { DefineAndReturnParticleTile(lev, grid_id, tile_id); } + do_boosted = WarpX::do_boosted_frame_diagnostic + && do_boosted_frame_diags; + auto old_size = particle_tile.GetArrayOfStructs().size(); auto new_size = old_size + max_new_particles; particle_tile.resize(new_size); @@ -439,7 +456,11 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) pb[4] = soa.GetRealData(particle_comps["uyold"]).data() + old_size; pb[5] = soa.GetRealData(particle_comps["uzold"]).data() + old_size; } - + int* pi; + if (do_field_ionization) { + pi = soa.GetIntData(particle_icomps["ionization_level"]).data() + old_size; + } + const GpuArray<Real,AMREX_SPACEDIM> overlap_corner {AMREX_D_DECL(overlap_realbox.lo(0), overlap_realbox.lo(1), @@ -448,6 +469,9 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) std::size_t shared_mem_bytes = plasma_injector->sharedMemoryNeeded(); int lrrfac = rrfac; + bool loc_do_field_ionization = do_field_ionization; + int loc_ionization_initial_level = ionization_initial_level; + // Loop over all new particles and inject them (creates too many // particles, in particular does not consider xmin, xmax etc.). // The invalid ones are given negative ID and are deleted during the @@ -568,6 +592,10 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) u.z = gamma_boost * ( u.z -beta_boost*gamma_lab ); } + if (loc_do_field_ionization) { + pi[ip] = loc_ionization_initial_level; + } + u.x *= PhysConst::c; u.y *= PhysConst::c; u.z *= PhysConst::c; @@ -1157,9 +1185,18 @@ PhysicalParticleContainer::Evolve (int lev, BL_PROFILE_VAR_STOP(blp_copy); if (rho) { - DepositCharge(pti, wp, rho, 0, 0, np_current, thread_num, lev, lev); + // Deposit charge before particle push, in component 0 of MultiFab rho. + int* AMREX_RESTRICT ion_lev; + if (do_field_ionization){ + ion_lev = pti.GetiAttribs(particle_icomps["ionization_level"]).dataPtr(); + } else { + ion_lev = nullptr; + } + DepositCharge(pti, wp, ion_lev, rho, 0, 0, + np_current, thread_num, lev, lev); if (has_buffer){ - DepositCharge(pti, wp, crho, 0, np_current, np-np_current, thread_num, lev, lev-1); + DepositCharge(pti, wp, ion_lev, crho, 0, np_current, + np-np_current, thread_num, lev, lev-1); } } @@ -1277,13 +1314,20 @@ PhysicalParticleContainer::Evolve (int lev, lev, lev-1, dt); } } else { + int* AMREX_RESTRICT ion_lev; + if (do_field_ionization){ + ion_lev = pti.GetiAttribs(particle_icomps["ionization_level"]).dataPtr(); + } else { + ion_lev = nullptr; + } + // Deposit inside domains - DepositCurrent(pti, wp, uxp, uyp, uzp, &jx, &jy, &jz, + DepositCurrent(pti, wp, uxp, uyp, uzp, ion_lev, &jx, &jy, &jz, 0, np_current, thread_num, lev, lev, dt); if (has_buffer){ // Deposit in buffers - DepositCurrent(pti, wp, uxp, uyp, uzp, cjx, cjy, cjz, + DepositCurrent(pti, wp, uxp, uyp, uzp, ion_lev, cjx, cjy, cjz, np_current, np-np_current, thread_num, lev, lev-1, dt); } @@ -1298,9 +1342,18 @@ PhysicalParticleContainer::Evolve (int lev, } if (rho) { - DepositCharge(pti, wp, rho, 1, 0, np_current, thread_num, lev, lev); + // Deposit charge after particle push, in component 1 of MultiFab rho. + int* AMREX_RESTRICT ion_lev; + if (do_field_ionization){ + ion_lev = pti.GetiAttribs(particle_icomps["ionization_level"]).dataPtr(); + } else { + ion_lev = nullptr; + } + DepositCharge(pti, wp, ion_lev, rho, 1, 0, + np_current, thread_num, lev, lev); if (has_buffer){ - DepositCharge(pti, wp, crho, 1, np_current, np-np_current, thread_num, lev, lev-1); + DepositCharge(pti, wp, ion_lev, crho, 1, np_current, + np-np_current, thread_num, lev, lev-1); } } @@ -1505,25 +1558,38 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti, copy_attribs(pti, x, y, z); } + int* AMREX_RESTRICT ion_lev = nullptr; + if (do_field_ionization){ + ion_lev = pti.GetiAttribs(particle_icomps["ionization_level"]).dataPtr(); + } + // Loop over the particles and update their momentum const Real q = this->charge; const Real m = this-> mass; if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Boris){ - amrex::ParallelFor( pti.numParticles(), + amrex::ParallelFor( + pti.numParticles(), [=] AMREX_GPU_DEVICE (long i) { + Real qp = q; + if (ion_lev){ qp *= ion_lev[i]; } UpdateMomentumBoris( ux[i], uy[i], uz[i], gi[i], - Ex[i], Ey[i], Ez[i], Bx[i], By[i], Bz[i], q, m, dt); + 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 if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Vay) { - amrex::ParallelFor( pti.numParticles(), + amrex::ParallelFor( + pti.numParticles(), [=] AMREX_GPU_DEVICE (long i) { + Real qp = q; + if (ion_lev){ qp *= ion_lev[i]; } UpdateMomentumVay( ux[i], uy[i], uz[i], gi[i], - Ex[i], Ey[i], Ez[i], Bx[i], By[i], Bz[i], q, m, dt); + 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 ); + ux[i], uy[i], uz[i], dt ); } ); } else { @@ -1931,3 +1997,136 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti, } } } + + +void PhysicalParticleContainer::InitIonizationModule () +{ + if (!do_field_ionization) return; + ParmParse pp(species_name); + if (charge != PhysConst::q_e){ + amrex::Warning( + "charge != q_e for ionizable species: overriding user value and setting charge = q_e."); + charge = PhysConst::q_e; + } + pp.query("ionization_initial_level", ionization_initial_level); + pp.get("ionization_product_species", ionization_product_name); + pp.get("physical_element", physical_element); + // Add runtime integer component for ionization level + AddIntComp("ionization_level"); + // Get atomic number and ionization energies from file + int ion_element_id = ion_map_ids[physical_element]; + ion_atomic_number = ion_atomic_numbers[ion_element_id]; + ionization_energies.resize(ion_atomic_number); + int offset = ion_energy_offsets[ion_element_id]; + for(int i=0; i<ion_atomic_number; i++){ + ionization_energies[i] = table_ionization_energies[i+offset]; + } + // Compute ADK prefactors (See Chen, JCP 236 (2013), equation (2)) + // For now, we assume l=0 and m=0. + // The approximate expressions are used, + // without Gamma function + Real wa = std::pow(PhysConst::alpha,3) * PhysConst::c / PhysConst::r_e; + Real Ea = PhysConst::m_e * PhysConst::c*PhysConst::c /PhysConst::q_e * + std::pow(PhysConst::alpha,4)/PhysConst::r_e; + Real UH = table_ionization_energies[0]; + Real l_eff = std::sqrt(UH/ionization_energies[0]) - 1.; + + const Real dt = WarpX::GetInstance().getdt(0); + + adk_power.resize(ion_atomic_number); + adk_prefactor.resize(ion_atomic_number); + adk_exp_prefactor.resize(ion_atomic_number); + for (int i=0; i<ion_atomic_number; ++i){ + Real n_eff = (i+1) * std::sqrt(UH/ionization_energies[i]); + Real C2 = std::pow(2,2*n_eff)/(n_eff*tgamma(n_eff+l_eff+1)*tgamma(n_eff-l_eff)); + adk_power[i] = -(2*n_eff - 1); + Real Uion = ionization_energies[i]; + adk_prefactor[i] = dt * wa * C2 * ( Uion/(2*UH) ) + * std::pow(2*std::pow((Uion/UH),3./2)*Ea,2*n_eff - 1); + adk_exp_prefactor[i] = -2./3 * std::pow( Uion/UH,3./2) * Ea; + } +} + +/* \brief create mask of ionized particles (1 if ionized, 0 otherwise) + * + * \param mfi: tile or grid + * \param lev: MR level + * \param ionization_mask: Array with as many elements as particles in mfi. + * This function initialized the array, and set each element to 1 or 0 + * depending on whether the particle is ionized or not. + */ +void +PhysicalParticleContainer::buildIonizationMask (const amrex::MFIter& mfi, const int lev, + amrex::Gpu::ManagedDeviceVector<int>& ionization_mask) +{ + BL_PROFILE("PPC::buildIonizationMask"); + // Get pointers to ionization data from pc_source + const Real * const AMREX_RESTRICT p_ionization_energies = ionization_energies.dataPtr(); + const Real * const AMREX_RESTRICT p_adk_prefactor = adk_prefactor.dataPtr(); + const Real * const AMREX_RESTRICT p_adk_exp_prefactor = adk_exp_prefactor.dataPtr(); + const Real * const AMREX_RESTRICT p_adk_power = adk_power.dataPtr(); + + // Current tile info + const int grid_id = mfi.index(); + const int tile_id = mfi.LocalTileIndex(); + + // Get GPU-friendly arrays of particle data + auto& ptile = GetParticles(lev)[std::make_pair(grid_id,tile_id)]; + // Only need attribs (i.e., SoA data) + auto& soa = ptile.GetStructOfArrays(); + const int np = ptile.GetArrayOfStructs().size(); + + // If no particle, nothing to do. + if (np == 0) return; + // Otherwise, resize ionization_mask, and get poiters to attribs arrays. + ionization_mask.resize(np); + int * const AMREX_RESTRICT p_ionization_mask = ionization_mask.data(); + const Real * const AMREX_RESTRICT ux = soa.GetRealData(PIdx::ux).data(); + const Real * const AMREX_RESTRICT uy = soa.GetRealData(PIdx::uy).data(); + const Real * const AMREX_RESTRICT uz = soa.GetRealData(PIdx::uz).data(); + const Real * const AMREX_RESTRICT ex = soa.GetRealData(PIdx::Ex).data(); + const Real * const AMREX_RESTRICT ey = soa.GetRealData(PIdx::Ey).data(); + const Real * const AMREX_RESTRICT ez = soa.GetRealData(PIdx::Ez).data(); + const Real * const AMREX_RESTRICT bx = soa.GetRealData(PIdx::Bx).data(); + const Real * const AMREX_RESTRICT by = soa.GetRealData(PIdx::By).data(); + const Real * const AMREX_RESTRICT bz = soa.GetRealData(PIdx::Bz).data(); + int* ion_lev = soa.GetIntData(particle_icomps["ionization_level"]).data(); + + Real c = PhysConst::c; + Real c2_inv = 1./c/c; + int atomic_number = ion_atomic_number; + + // Loop over all particles in grid/tile. If ionized, set mask to 1 + // and increment ionization level. + ParallelFor( + np, + [=] AMREX_GPU_DEVICE (long i) { + // Get index of ionization_level + p_ionization_mask[i] = 0; + if ( ion_lev[i]<atomic_number ){ + Real random_draw = amrex::Random(); + // Compute electric field amplitude in the particle's frame of + // reference (particularly important when in boosted frame). + Real ga = std::sqrt(1. + (ux[i]*ux[i] + uy[i]*uy[i] + uz[i]*uz[i]) * c2_inv); + Real E = std::sqrt( + - ( ux[i]*ex[i] + uy[i]*ey[i] + uz[i]*ez[i] ) * ( ux[i]*ex[i] + uy[i]*ey[i] + uz[i]*ez[i] ) * c2_inv + + ( ga *ex[i] + uy[i]*bz[i] - uz[i]*by[i] ) * ( ga *ex[i] + uy[i]*bz[i] - uz[i]*by[i] ) + + ( ga *ey[i] + uz[i]*bx[i] - ux[i]*bz[i] ) * ( ga *ey[i] + uz[i]*bx[i] - ux[i]*bz[i] ) + + ( ga *ez[i] + ux[i]*by[i] - uy[i]*bx[i] ) * ( ga *ez[i] + ux[i]*by[i] - uy[i]*bx[i] ) + ); + // Compute probability of ionization p + Real w_dtau = 1./ ga * p_adk_prefactor[ion_lev[i]] * + std::pow(E,p_adk_power[ion_lev[i]]) * + std::exp( p_adk_exp_prefactor[ion_lev[i]]/E ); + Real p = 1. - std::exp( - w_dtau ); + + if (random_draw < p){ + // increment particle's ionization level + ion_lev[i] += 1; + // update mask + p_ionization_mask[i] = 1; + } + } + } + ); +} diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index ac5b47ada..c252eee6b 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -76,6 +76,10 @@ public: RealVector& GetAttribs (int comp) { return GetStructOfArrays().GetRealData(comp); } + + IntVector& GetiAttribs (int comp) { + return GetStructOfArrays().GetIntData(comp); + } }; class MultiParticleContainer; @@ -167,6 +171,7 @@ public: virtual void DepositCharge(WarpXParIter& pti, RealVector& wp, + const int * const ion_lev, amrex::MultiFab* rho, int icomp, const long offset, @@ -180,6 +185,7 @@ public: RealVector& uxp, RealVector& uyp, RealVector& uzp, + const int * const ion_lev, amrex::MultiFab* jx, amrex::MultiFab* jy, amrex::MultiFab* jz, @@ -249,6 +255,8 @@ public: static int NextID () { return ParticleType::NextID(); } + void setNextID(int next_id) { ParticleType::NextID(next_id); } + bool do_splitting = false; // split along axes (0) or diagonals (1) @@ -265,15 +273,22 @@ public: void AddIntComp (const std::string& name, bool comm=true) { - particle_comps[name] = NumIntComps(); + particle_icomps[name] = NumIntComps(); AddIntComp(comm); } int DoBoostedFrameDiags () const { return do_boosted_frame_diags; } + 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;} + protected: std::map<std::string, int> particle_comps; + std::map<std::string, int> particle_icomps; int species_id; @@ -290,6 +305,17 @@ protected: // support all features allowed by direct injection. int do_continuous_injection = 0; + int do_field_ionization = 0; + int ionization_product; + std::string ionization_product_name; + int ion_atomic_number; + int ionization_initial_level = 0; + amrex::Gpu::ManagedVector<amrex::Real> ionization_energies; + amrex::Gpu::ManagedVector<amrex::Real> adk_power; + amrex::Gpu::ManagedVector<amrex::Real> adk_prefactor; + amrex::Gpu::ManagedVector<amrex::Real> adk_exp_prefactor; + std::string physical_element; + int do_boosted_frame_diags = 1; amrex::Vector<amrex::FArrayBox> local_rho; diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index befa5cfed..af56c58fb 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -238,12 +238,14 @@ WarpXParticleContainer::AddNParticles (int lev, p.pos(1) = z[i]; #endif - if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) - { - auto& ptile = DefineAndReturnParticleTile(0, 0, 0); - ptile.push_back_real(particle_comps["xold"], x[i]); - ptile.push_back_real(particle_comps["yold"], y[i]); - ptile.push_back_real(particle_comps["zold"], z[i]); + if ( (NumRuntimeRealComps()>0) || (NumRuntimeIntComps()>0) ){ + auto& ptile = DefineAndReturnParticleTile(0, 0, 0); + if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) + { + ptile.push_back_real(particle_comps["xold"], x[i]); + ptile.push_back_real(particle_comps["yold"], y[i]); + ptile.push_back_real(particle_comps["zold"], z[i]); + } } particle_tile.push_back(p); @@ -256,12 +258,14 @@ WarpXParticleContainer::AddNParticles (int lev, particle_tile.push_back_real(PIdx::uy, vy + ibegin, vy + iend); particle_tile.push_back_real(PIdx::uz, vz + ibegin, vz + iend); - if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) - { - auto& ptile = DefineAndReturnParticleTile(0, 0, 0); - ptile.push_back_real(particle_comps["uxold"], vx + ibegin, vx + iend); - ptile.push_back_real(particle_comps["uyold"], vy + ibegin, vy + iend); - ptile.push_back_real(particle_comps["uzold"], vz + ibegin, vz + iend); + if ( (NumRuntimeRealComps()>0) || (NumRuntimeIntComps()>0) ){ + auto& ptile = DefineAndReturnParticleTile(0, 0, 0); + if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) + { + ptile.push_back_real(particle_comps["uxold"], vx + ibegin, vx + iend); + ptile.push_back_real(particle_comps["uyold"], vy + ibegin, vy + iend); + ptile.push_back_real(particle_comps["uzold"], vz + ibegin, vz + iend); + } } for (int comp = PIdx::uz+1; comp < PIdx::nattribs; ++comp) @@ -412,6 +416,10 @@ WarpXParticleContainer::DepositCurrentFortran(WarpXParIter& pti, * \param pti : Particle iterator * \param wp : Array of particle weights * \param uxp uyp uzp : Array of particle + * \param ion_lev : Pointer to array of particle ionization level. This is + required to have the charge of each macroparticle + since q is a scalar. For non-ionizable species, + ion_lev is a null pointer. * \param jx jy jz : Full array of current density * \param offset : Index of first particle for which current is deposited * \param np_to_depose: Number of particles for which current is deposited. @@ -425,6 +433,7 @@ void WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, RealVector& wp, RealVector& uxp, RealVector& uyp, RealVector& uzp, + const int * const ion_lev, MultiFab* jx, MultiFab* jy, MultiFab* jz, const long offset, const long np_to_depose, int thread_num, int lev, int depos_lev, @@ -507,37 +516,40 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, if (WarpX::current_deposition_algo == CurrentDepositionAlgo::Esirkepov) { if (WarpX::nox == 1){ - doEsirkepovDepositionShapeN<1>(xp, yp, zp, wp.dataPtr() + offset, uxp.dataPtr() + offset, - uyp.dataPtr() + offset, uzp.dataPtr() + offset, jx_arr, jy_arr, - jz_arr, np_to_depose, dt, dx, - xyzmin, lo, q); + doEsirkepovDepositionShapeN<1>( + xp, yp, zp, wp.dataPtr() + offset, uxp.dataPtr() + offset, + uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev, + jx_arr, jy_arr, jz_arr, np_to_depose, dt, dx, xyzmin, lo, q); } else if (WarpX::nox == 2){ - doEsirkepovDepositionShapeN<2>(xp, yp, zp, wp.dataPtr() + offset, uxp.dataPtr() + offset, - uyp.dataPtr() + offset, uzp.dataPtr() + offset, jx_arr, jy_arr, - jz_arr, np_to_depose, dt, dx, - xyzmin, lo, q); + doEsirkepovDepositionShapeN<2>( + xp, yp, zp, wp.dataPtr() + offset, uxp.dataPtr() + offset, + uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev, + jx_arr, jy_arr, jz_arr, np_to_depose, dt, dx, xyzmin, lo, q); } else if (WarpX::nox == 3){ - doEsirkepovDepositionShapeN<3>(xp, yp, zp, wp.dataPtr() + offset, uxp.dataPtr() + offset, - uyp.dataPtr() + offset, uzp.dataPtr() + offset, jx_arr, jy_arr, - jz_arr, np_to_depose, dt, dx, - xyzmin, lo, q); + doEsirkepovDepositionShapeN<3>( + xp, yp, zp, wp.dataPtr() + offset, uxp.dataPtr() + offset, + uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev, + jx_arr, jy_arr, jz_arr, np_to_depose, dt, dx, xyzmin, lo, q); } } else { if (WarpX::nox == 1){ - doDepositionShapeN<1>(xp, yp, zp, wp.dataPtr() + offset, uxp.dataPtr() + offset, - uyp.dataPtr() + offset, uzp.dataPtr() + offset, jx_arr, jy_arr, - jz_arr, np_to_depose, dt, dx, - xyzmin, lo, stagger_shift, q); + doDepositionShapeN<1>( + xp, yp, zp, wp.dataPtr() + offset, uxp.dataPtr() + offset, + uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev, + jx_arr, jy_arr, jz_arr, np_to_depose, dt, dx, xyzmin, lo, + stagger_shift, q); } else if (WarpX::nox == 2){ - doDepositionShapeN<2>(xp, yp, zp, wp.dataPtr() + offset, uxp.dataPtr() + offset, - uyp.dataPtr() + offset, uzp.dataPtr() + offset, jx_arr, jy_arr, - jz_arr, np_to_depose, dt, dx, - xyzmin, lo, stagger_shift, q); + doDepositionShapeN<2>( + xp, yp, zp, wp.dataPtr() + offset, uxp.dataPtr() + offset, + uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev, + jx_arr, jy_arr, jz_arr, np_to_depose, dt, dx, xyzmin, lo, + stagger_shift, q); } else if (WarpX::nox == 3){ - doDepositionShapeN<3>(xp, yp, zp, wp.dataPtr() + offset, uxp.dataPtr() + offset, - uyp.dataPtr() + offset, uzp.dataPtr() + offset, jx_arr, jy_arr, - jz_arr, np_to_depose, dt, dx, - xyzmin, lo, stagger_shift, q); + doDepositionShapeN<3>( + xp, yp, zp, wp.dataPtr() + offset, uxp.dataPtr() + offset, + uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev, + jx_arr, jy_arr, jz_arr, np_to_depose, dt, dx, xyzmin, lo, + stagger_shift, q); } } @@ -552,8 +564,27 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, #endif } +/* \brief Charge Deposition for thread thread_num + * \param pti : Particle iterator + * \param wp : Array of particle weights + * \param ion_lev : Pointer to array of particle ionization level. This is + required to have the charge of each macroparticle + since q is a scalar. For non-ionizable species, + ion_lev is a null pointer. + * \param rho : Full array of charge density + * \param icomp : Component of rho into which charge is deposited. + 0: old value (before particle push). + 1: new value (after particle push). + * \param offset : Index of first particle for which charge is deposited + * \param np_to_depose: Number of particles for which charge is deposited. + Particles [offset,offset+np_tp_depose] deposit charge + * \param thread_num : Thread number (if tiling) + * \param lev : Level of box that contains particles + * \param depos_lev : Level on which particles deposit (if buffers are used) + */ void WarpXParticleContainer::DepositCharge (WarpXParIter& pti, RealVector& wp, + const int * const ion_lev, MultiFab* rho, int icomp, const long offset, const long np_to_depose, int thread_num, int lev, int depos_lev) @@ -615,14 +646,14 @@ WarpXParticleContainer::DepositCharge (WarpXParIter& pti, RealVector& wp, BL_PROFILE_VAR_START(blp_ppc_chd); if (WarpX::nox == 1){ - doChargeDepositionShapeN<1>(xp, yp, zp, wp.dataPtr()+offset, rho_arr, - np_to_depose, dx, xyzmin, lo, q); + doChargeDepositionShapeN<1>(xp, yp, zp, wp.dataPtr()+offset, ion_lev, + rho_arr, np_to_depose, dx, xyzmin, lo, q); } else if (WarpX::nox == 2){ - doChargeDepositionShapeN<2>(xp, yp, zp, wp.dataPtr()+offset, rho_arr, - np_to_depose, dx, xyzmin, lo, q); + doChargeDepositionShapeN<2>(xp, yp, zp, wp.dataPtr()+offset, ion_lev, + rho_arr, np_to_depose, dx, xyzmin, lo, q); } else if (WarpX::nox == 3){ - doChargeDepositionShapeN<3>(xp, yp, zp, wp.dataPtr()+offset, rho_arr, - np_to_depose, dx, xyzmin, lo, q); + doChargeDepositionShapeN<3>(xp, yp, zp, wp.dataPtr()+offset, ion_lev, + rho_arr, np_to_depose, dx, xyzmin, lo, q); } BL_PROFILE_VAR_STOP(blp_ppc_chd); @@ -732,7 +763,15 @@ WarpXParticleContainer::GetChargeDensity (int lev, bool local) pti.GetPosition(m_xp[thread_num], m_yp[thread_num], m_zp[thread_num]); - DepositCharge(pti, wp, rho.get(), 0, 0, np, thread_num, lev, lev); + int* AMREX_RESTRICT ion_lev; + if (do_field_ionization){ + ion_lev = pti.GetiAttribs(particle_icomps["ionization_level"]).dataPtr(); + } else { + ion_lev = nullptr; + } + + DepositCharge(pti, wp, ion_lev, rho.get(), 0, 0, np, + thread_num, lev, lev); } #ifdef _OPENMP } |