From 58d728b3160c69ce8e9f117b255e7b057cc419fe Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Sat, 3 Aug 2019 17:45:23 -0700 Subject: add ionization initialization --- Source/Particles/MultiParticleContainer.cpp | 51 +++++++++++++++++++++++++++++ 1 file changed, 51 insertions(+) (limited to 'Source/Particles/MultiParticleContainer.cpp') diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 982e04e39..16c9b9762 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -1,10 +1,12 @@ #include #include #include +#include #include #include #include +#include using namespace amrex; @@ -29,6 +31,10 @@ MultiParticleContainer::MultiParticleContainer (AmrCore* amr_core) pc_tmp.reset(new PhysicalParticleContainer(amr_core)); + // For each species, get the ID of its product species. + // This is used for ionization and pair creation processes. + mapSpeciesProduct(); + // 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. @@ -461,3 +467,48 @@ 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; iionization_product_name and store its ID into + // pc->ionization_product. + if (pc->do_field_ionization){ + int i_product = getSpeciesID(pc->ionization_product_name); + 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 Date: Sun, 4 Aug 2019 15:47:12 -0700 Subject: add target species handling --- Source/Particles/MultiParticleContainer.H | 2 +- Source/Particles/MultiParticleContainer.cpp | 155 ++++++++++++++++++++++++- Source/Particles/PhysicalParticleContainer.H | 2 + Source/Particles/PhysicalParticleContainer.cpp | 21 +++- Source/Particles/WarpXParticleContainer.H | 5 +- 5 files changed, 179 insertions(+), 6 deletions(-) (limited to 'Source/Particles/MultiParticleContainer.cpp') diff --git a/Source/Particles/MultiParticleContainer.H b/Source/Particles/MultiParticleContainer.H index 4ea9c2763..8f7fd56e3 100644 --- a/Source/Particles/MultiParticleContainer.H +++ b/Source/Particles/MultiParticleContainer.H @@ -7,7 +7,7 @@ #include #include - +#include #include #include #include diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 16c9b9762..736b906df 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -3,6 +3,9 @@ #include #include +#include +#include +#include #include #include #include @@ -508,7 +511,155 @@ MultiParticleContainer::getSpeciesID(std::string product_str) return i_product; } - void MultiParticleContainer::doFieldIonization() -{} +{ + /* + amrex::Gpu::SharedMemory grid_ids; + amrex::Gpu::SharedMemory tile_ids; + amrex::Gpu::SharedMemory is_ionized; + */ + + for (auto& pc : allcontainers){ + if (!pc->do_field_ionization){ continue; } + const Real * const AMREX_RESTRICT p_ionization_energies = pc->ionization_energies.dataPtr(); + const Real * const AMREX_RESTRICT p_adk_prefactor = pc->adk_prefactor.dataPtr(); + const Real * const AMREX_RESTRICT p_adk_exp_prefactor = pc->adk_exp_prefactor.dataPtr(); + const Real * const AMREX_RESTRICT p_adk_power = pc->adk_power.dataPtr(); + for (int lev = 0; lev <= pc->finestLevel(); ++lev){ + +#ifdef _OPENMP + // First touch all tiles in the map in serial + for (MFIter mfi = pc->MakeMFIter(lev); mfi.isValid(); ++mfi) { + const int grid_id = mfi.index(); + const int tile_id = mfi.LocalTileIndex(); + pc->GetParticles(lev)[std::make_pair(grid_id,tile_id)]; + if ( (pc->NumRuntimeRealComps()>0) || (pc->NumRuntimeIntComps()>0) ) { + pc->DefineAndReturnParticleTile(lev, grid_id, tile_id); + } + } +#endif + MFItInfo info; + if (pc->do_tiling && Gpu::notInLaunchRegion()) { + info.EnableTiling(pc->tile_size); + } +#ifdef _OPENMP + info.SetDynamic(true); +#pragma omp parallel +#endif + for (MFIter mfi = pc->MakeMFIter(lev, info); mfi.isValid(); ++mfi) + { + const int grid_id = mfi.index(); + const int tile_id = mfi.LocalTileIndex(); + auto& particle_tile = pc->GetParticles(lev)[std::make_pair(grid_id,tile_id)]; + auto& soa = particle_tile.GetStructOfArrays(); + const int np = particle_tile.GetArrayOfStructs().size(); + if (np == 0) break; + amrex::Gpu::ManagedVector is_ionized; + // const int np_lev = pc->NumberOfParticlesAtLevel(lev); + is_ionized.resize(np); + int * const AMREX_RESTRICT p_is_ionized = is_ionized.dataPtr(); + 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(); + Real * const AMREX_RESTRICT ilev_real = soa.GetRealData(pc->particle_comps["ionization_level"]).data(); + Real c = PhysConst::c; + Real c2_inv = 1./c/c; + + ParallelFor( + np, + [=] AMREX_GPU_DEVICE (long i) { + Real random_draw = Random(); + 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] ) + ); + int ilev = (int) round(ilev_real[i]); + // int ilev = static_cast(round(ilev_real[i])); + Real p; + Real w_dtau; + Print()< is_ionized_cumsum_vector; + is_ionized_cumsum_vector.resize(np); + int np_ionized = p_is_ionized[0]; + for(int i=1; iionization_product]; + auto& prod_particle_tile = prod_pc->GetParticles(lev)[std::make_pair(grid_id,tile_id)]; + const int np_old = prod_particle_tile.GetArrayOfStructs().size(); + const int np_new = np_old + np_ionized; + prod_particle_tile.resize(np_new); + auto& prod_soa = prod_particle_tile.GetStructOfArrays(); + + GpuArray prod_pa; + for (int ia = 0; ia < PIdx::nattribs; ++ia) { + prod_pa[ia] = prod_soa.GetRealData(ia).data() + np_old; + } + WarpXParticleContainer::ParticleType* prod_pp = prod_particle_tile.GetArrayOfStructs()().data() + np_old; + WarpXParticleContainer::ParticleType* pp = particle_tile.GetArrayOfStructs()().data(); + + GpuArray pa; + for (int ia = 0; ia < PIdx::nattribs; ++ia) { + pa[ia] = soa.GetRealData(ia).data(); + } + + const int cpuid = ParallelDescriptor::MyProc(); + int* AMREX_RESTRICT is_ionized_cumsum = is_ionized_cumsum_vector.dataPtr(); + amrex::For( + np, [=] AMREX_GPU_DEVICE (int ip) noexcept + { + if(is_ionized[ip]){ + int i = is_ionized_cumsum[ip]; + WarpXParticleContainer::ParticleType& prod_p = prod_pp[i]; + WarpXParticleContainer::ParticleType& p = pp[ip]; + prod_p.id() = 9999; + prod_p.cpu() = cpuid; + prod_p.pos(0) = p.pos(0); + prod_p.pos(1) = p.pos(1); +#if (AMREX_SPACEDIM == 3) + prod_p.pos(20 = p.pos(2); +#endif + for (int ia = 0; ia < PIdx::nattribs; ++ia) { + prod_pa[ia][i] = pa[ia][ip]; + } + } + } + ); + } // MFIter + } + } +} diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index a6cf9d81f..6f9fdc485 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -108,6 +108,8 @@ public: virtual void PostRestart () final {} void SplitParticles(int lev); + + virtual int doFieldIonization (const int lev) 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 0183a44b9..8232b621f 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -215,7 +215,7 @@ PhysicalParticleContainer::CheckAndAddParticle(Real x, Real y, Real z, auto& particle_tile = DefineAndReturnParticleTile(0, 0, 0); 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["xold"], x); particle_tile.push_back_real(particle_comps["yold"], y); particle_tile.push_back_real(particle_comps["zold"], z); @@ -591,7 +591,7 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) } if (do_field_ionization) { - pi[ip] = 2.; // species_ionization_level; + pi[ip] = species_ionization_level; } u.x *= PhysConst::c; @@ -2105,4 +2105,21 @@ void PhysicalParticleContainer::InitIonizationModule() adk_power.resize(ion_atomic_number); adk_prefactor.resize(ion_atomic_number); adk_exp_prefactor.resize(ion_atomic_number); + for (int i=0; i Date: Sun, 4 Aug 2019 15:53:18 -0700 Subject: minor cleaning --- Source/Particles/MultiParticleContainer.cpp | 1 - 1 file changed, 1 deletion(-) (limited to 'Source/Particles/MultiParticleContainer.cpp') diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 736b906df..1af169d1c 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -590,7 +590,6 @@ MultiParticleContainer::doFieldIonization() // int ilev = static_cast(round(ilev_real[i])); Real p; Real w_dtau; - Print()< Date: Mon, 5 Aug 2019 09:52:54 -0700 Subject: typo --- Source/Particles/MultiParticleContainer.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'Source/Particles/MultiParticleContainer.cpp') diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 1af169d1c..329a43ca7 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -650,7 +650,7 @@ MultiParticleContainer::doFieldIonization() prod_p.pos(0) = p.pos(0); prod_p.pos(1) = p.pos(1); #if (AMREX_SPACEDIM == 3) - prod_p.pos(20 = p.pos(2); + prod_p.pos(2) = p.pos(2); #endif for (int ia = 0; ia < PIdx::nattribs; ++ia) { prod_pa[ia][i] = pa[ia][ip]; -- cgit v1.2.3 From dd12326a323019abda4cb20d1d5d49e332084fad Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Mon, 5 Aug 2019 17:48:59 -0700 Subject: get dt from WarpX instance, to initialize ionization prefactors --- Source/Particles/MultiParticleContainer.cpp | 45 ++++++++++++++------------ Source/Particles/PhysicalParticleContainer.cpp | 18 ++++++++--- 2 files changed, 38 insertions(+), 25 deletions(-) (limited to 'Source/Particles/MultiParticleContainer.cpp') diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 329a43ca7..a79670582 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -34,10 +34,6 @@ MultiParticleContainer::MultiParticleContainer (AmrCore* amr_core) pc_tmp.reset(new PhysicalParticleContainer(amr_core)); - // For each species, get the ID of its product species. - // This is used for ionization and pair creation processes. - mapSpeciesProduct(); - // 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. @@ -144,6 +140,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(); } @@ -511,21 +510,23 @@ MultiParticleContainer::getSpeciesID(std::string product_str) return i_product; } +//void +//MultiParticleContainer::InitIonizationModule() + void MultiParticleContainer::doFieldIonization() { - /* - amrex::Gpu::SharedMemory grid_ids; - amrex::Gpu::SharedMemory tile_ids; - amrex::Gpu::SharedMemory is_ionized; - */ - for (auto& pc : allcontainers){ + if (!pc->do_field_ionization){ continue; } + + auto& prod_pc = allcontainers[pc->ionization_product]; + const Real * const AMREX_RESTRICT p_ionization_energies = pc->ionization_energies.dataPtr(); const Real * const AMREX_RESTRICT p_adk_prefactor = pc->adk_prefactor.dataPtr(); const Real * const AMREX_RESTRICT p_adk_exp_prefactor = pc->adk_exp_prefactor.dataPtr(); const Real * const AMREX_RESTRICT p_adk_power = pc->adk_power.dataPtr(); + for (int lev = 0; lev <= pc->finestLevel(); ++lev){ #ifdef _OPENMP @@ -537,6 +538,10 @@ MultiParticleContainer::doFieldIonization() if ( (pc->NumRuntimeRealComps()>0) || (pc->NumRuntimeIntComps()>0) ) { pc->DefineAndReturnParticleTile(lev, grid_id, tile_id); } + prod_pc->GetParticles(lev)[std::make_pair(grid_id,tile_id)]; + if ( (prod_pc->NumRuntimeRealComps()>0) || (prod_pc->NumRuntimeIntComps()>0) ) { + prod_pc->DefineAndReturnParticleTile(lev, grid_id, tile_id); + } } #endif MFItInfo info; @@ -575,7 +580,7 @@ MultiParticleContainer::doFieldIonization() ParallelFor( np, [=] AMREX_GPU_DEVICE (long i) { - Real random_draw = Random(); + Real random_draw = amrex::Random(); Real ga = std::sqrt(1. + (ux[i]*ux[i] + uy[i]*uy[i] + @@ -589,15 +594,14 @@ MultiParticleContainer::doFieldIonization() int ilev = (int) round(ilev_real[i]); // int ilev = static_cast(round(ilev_real[i])); Real p; - Real w_dtau; - if (E<1.e-100*(p_ionization_energies[0])){ - p = 0.; - } else { - w_dtau = 1./ ga * p_adk_prefactor[ilev] * - std::pow(E,p_adk_power[ilev]) * - std::exp( p_adk_exp_prefactor[ilev]/E ); - p = 1. - std::exp( - w_dtau ); - } + //if (E<1.e-20*(p_ionization_energies[0])){ + // p = 0.; + //} else { + Real w_dtau = 1./ ga * p_adk_prefactor[ilev] * + std::pow(E,p_adk_power[ilev]) * + std::exp( p_adk_exp_prefactor[ilev]/E ); + p = 1. - std::exp( - w_dtau ); + //} p_is_ionized[i] = 0; if (random_draw < p){ ilev_real[i] += 1.; @@ -617,7 +621,6 @@ MultiParticleContainer::doFieldIonization() break; } - auto& prod_pc = allcontainers[pc->ionization_product]; auto& prod_particle_tile = prod_pc->GetParticles(lev)[std::make_pair(grid_id,tile_id)]; const int np_old = prod_particle_tile.GetArrayOfStructs().size(); const int np_new = np_old + np_ionized; diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 8232b621f..fa01f472d 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -41,8 +41,8 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp pp.query("do_field_ionization", do_field_ionization); // If do_field_ionization, read initialization data from input file and // read ionization energies from table. - if (do_field_ionization) - InitIonizationModule(); + //if (do_field_ionization) + // InitIonizationModule(); pp.query("plot_species", plot_species); int do_user_plot_vars; @@ -90,6 +90,7 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core) void PhysicalParticleContainer::InitData() { + if (do_field_ionization) {InitIonizationModule();} AddParticles(0); // Note - add on level 0 Redistribute(); // We then redistribute } @@ -2078,6 +2079,7 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti, void PhysicalParticleContainer::InitIonizationModule() { + if (!do_field_ionization) return; ParmParse pp(species_name); pp.query("ionization_level", species_ionization_level); pp.get("ionization_product", ionization_product_name); @@ -2092,16 +2094,22 @@ void PhysicalParticleContainer::InitIonizationModule() int offset = ion_energy_offsets[ion_element_id]; for(int i=0; i Date: Tue, 6 Aug 2019 13:07:57 -0700 Subject: cleaning, and add capability for boosted frame runtime particle components --- Source/Particles/MultiParticleContainer.cpp | 318 +++++++++++++++---------- Source/Particles/PhysicalParticleContainer.H | 3 + Source/Particles/PhysicalParticleContainer.cpp | 82 +++++++ Source/Particles/WarpXParticleContainer.H | 8 + 4 files changed, 283 insertions(+), 128 deletions(-) (limited to 'Source/Particles/MultiParticleContainer.cpp') diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index a79670582..6c6a03231 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -513,155 +513,217 @@ MultiParticleContainer::getSpeciesID(std::string product_str) //void //MultiParticleContainer::InitIonizationModule() +namespace +{ + static void createIonizedParticles( + int lev, const MFIter& mfi, + std::unique_ptr< WarpXParticleContainer>& pc_source, + std::unique_ptr< WarpXParticleContainer>& pc_product, + const int * const p_is_ionized) + { + + 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(); + // --- source AoS particle data + WarpXParticleContainer::ParticleType* particles_source = ptile_source.GetArrayOfStructs()().data(); + // --- source SoA particle data + auto& soa_source = ptile_source.GetStructOfArrays(); + GpuArray attribs_source; + for (int ia = 0; ia < PIdx::nattribs; ++ia) { + attribs_source[ia] = soa_source.GetRealData(ia).data(); + } + // --- source runtime attribs + bool do_boosted_product = WarpX::do_boosted_frame_diagnostic + && pc_product->DoBoostedFrameDiags(); + bool do_boosted_source = WarpX::do_boosted_frame_diagnostic + && pc_source->DoBoostedFrameDiags(); + GpuArray runtime_attribs_source; + // + if (do_boosted_product && do_boosted_source) { + // If boosted frame diagnostics for source species, store them + std::map comps_source = pc_source->getParticleComps(); + runtime_attribs_source[0] = soa_source.GetRealData(comps_source[ "xold"]).data(); + runtime_attribs_source[1] = soa_source.GetRealData(comps_source[ "yold"]).data(); + runtime_attribs_source[2] = soa_source.GetRealData(comps_source[ "zold"]).data(); + runtime_attribs_source[3] = soa_source.GetRealData(comps_source["uxold"]).data(); + runtime_attribs_source[4] = soa_source.GetRealData(comps_source["uyold"]).data(); + runtime_attribs_source[5] = soa_source.GetRealData(comps_source["uzold"]).data(); + } else if (do_boosted_product && !do_boosted_source){ + // Otherwise, store current particle momenta. + // Positions are copied from AoS data. + runtime_attribs_source[3] = soa_source.GetRealData(PIdx::ux).data(); + runtime_attribs_source[4] = soa_source.GetRealData(PIdx::uy).data(); + runtime_attribs_source[5] = soa_source.GetRealData(PIdx::uz).data(); + } + + // Indices of product particle for each ionized source particle. + // i_product[i] is the location in product tile of product particle + // from source particle i. + amrex::Gpu::ManagedVector i_product; + i_product.resize(np_source); + // 0GetParticles(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 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 runtime_attribs_product; + if (do_boosted_product) { + std::map 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]){ + int ip = p_i_product[is]; + // is: index of ionized particle in source species + // ip: index of corresponding created 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. For position, we need a different + // handling depending on do_boosted_source. For momentum, + // runtime_attribs_source[3-5] contains appropriate data. + if (do_boosted_product) { + if (do_boosted_source) { + runtime_attribs_product[0][ip] = runtime_attribs_source[0][ip]; + runtime_attribs_product[1][ip] = runtime_attribs_source[1][ip]; + runtime_attribs_product[2][ip] = runtime_attribs_source[2][ip]; + } else { + 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_attribs_source[3][ip]; + runtime_attribs_product[4][ip] = runtime_attribs_source[4][ip]; + runtime_attribs_product[5][ip] = runtime_attribs_source[5][ip]; + } + } + } + ); + } +} + void MultiParticleContainer::doFieldIonization() { - for (auto& pc : allcontainers){ - - if (!pc->do_field_ionization){ continue; } - auto& prod_pc = allcontainers[pc->ionization_product]; + // 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; } - const Real * const AMREX_RESTRICT p_ionization_energies = pc->ionization_energies.dataPtr(); - const Real * const AMREX_RESTRICT p_adk_prefactor = pc->adk_prefactor.dataPtr(); - const Real * const AMREX_RESTRICT p_adk_exp_prefactor = pc->adk_exp_prefactor.dataPtr(); - const Real * const AMREX_RESTRICT p_adk_power = pc->adk_power.dataPtr(); + // Get product species + auto& pc_product = allcontainers[pc_source->ionization_product]; - for (int lev = 0; lev <= pc->finestLevel(); ++lev){ + for (int lev = 0; lev <= pc_source->finestLevel(); ++lev){ #ifdef _OPENMP - // First touch all tiles in the map in serial - for (MFIter mfi = pc->MakeMFIter(lev); mfi.isValid(); ++mfi) { + // Touch all tiles of source species in serial if additional arguments + for (MFIter mfi = pc_source->MakeMFIter(lev); mfi.isValid(); ++mfi) { const int grid_id = mfi.index(); const int tile_id = mfi.LocalTileIndex(); - pc->GetParticles(lev)[std::make_pair(grid_id,tile_id)]; - if ( (pc->NumRuntimeRealComps()>0) || (pc->NumRuntimeIntComps()>0) ) { - pc->DefineAndReturnParticleTile(lev, grid_id, tile_id); - } - prod_pc->GetParticles(lev)[std::make_pair(grid_id,tile_id)]; - if ( (prod_pc->NumRuntimeRealComps()>0) || (prod_pc->NumRuntimeIntComps()>0) ) { - prod_pc->DefineAndReturnParticleTile(lev, grid_id, tile_id); + 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->do_tiling && Gpu::notInLaunchRegion()) { - info.EnableTiling(pc->tile_size); + 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 - for (MFIter mfi = pc->MakeMFIter(lev, info); mfi.isValid(); ++mfi) + // Loop over all grids (if not tiling) or grids and tiles (if tiling) + for (MFIter mfi = pc_source->MakeMFIter(lev, info); mfi.isValid(); ++mfi) { - const int grid_id = mfi.index(); - const int tile_id = mfi.LocalTileIndex(); - auto& particle_tile = pc->GetParticles(lev)[std::make_pair(grid_id,tile_id)]; - auto& soa = particle_tile.GetStructOfArrays(); - const int np = particle_tile.GetArrayOfStructs().size(); - if (np == 0) break; + // Ionization mask: one element per source particles. + // 0 if not ionized, 1 if ionized. amrex::Gpu::ManagedVector is_ionized; - // const int np_lev = pc->NumberOfParticlesAtLevel(lev); - is_ionized.resize(np); - int * const AMREX_RESTRICT p_is_ionized = is_ionized.dataPtr(); - 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(); - Real * const AMREX_RESTRICT ilev_real = soa.GetRealData(pc->particle_comps["ionization_level"]).data(); - Real c = PhysConst::c; - Real c2_inv = 1./c/c; - - ParallelFor( - np, - [=] AMREX_GPU_DEVICE (long i) { - Real random_draw = amrex::Random(); - 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] ) - ); - int ilev = (int) round(ilev_real[i]); - // int ilev = static_cast(round(ilev_real[i])); - Real p; - //if (E<1.e-20*(p_ionization_energies[0])){ - // p = 0.; - //} else { - Real w_dtau = 1./ ga * p_adk_prefactor[ilev] * - std::pow(E,p_adk_power[ilev]) * - std::exp( p_adk_exp_prefactor[ilev]/E ); - p = 1. - std::exp( - w_dtau ); - //} - p_is_ionized[i] = 0; - if (random_draw < p){ - ilev_real[i] += 1.; - p_is_ionized[i] = 1; - } - } - ); // ParallelFor - - amrex::Gpu::ManagedVector is_ionized_cumsum_vector; - is_ionized_cumsum_vector.resize(np); - int np_ionized = p_is_ionized[0]; - for(int i=1; iGetParticles(lev)[std::make_pair(grid_id,tile_id)]; - const int np_old = prod_particle_tile.GetArrayOfStructs().size(); - const int np_new = np_old + np_ionized; - prod_particle_tile.resize(np_new); - auto& prod_soa = prod_particle_tile.GetStructOfArrays(); - - GpuArray prod_pa; - for (int ia = 0; ia < PIdx::nattribs; ++ia) { - prod_pa[ia] = prod_soa.GetRealData(ia).data() + np_old; - } - WarpXParticleContainer::ParticleType* prod_pp = prod_particle_tile.GetArrayOfStructs()().data() + np_old; - WarpXParticleContainer::ParticleType* pp = particle_tile.GetArrayOfStructs()().data(); - - GpuArray pa; - for (int ia = 0; ia < PIdx::nattribs; ++ia) { - pa[ia] = soa.GetRealData(ia).data(); - } + pc_source->buildIonizationMask(mfi, lev, is_ionized); + const int * const AMREX_RESTRICT p_is_ionized = is_ionized.dataPtr(); - const int cpuid = ParallelDescriptor::MyProc(); - int* AMREX_RESTRICT is_ionized_cumsum = is_ionized_cumsum_vector.dataPtr(); - amrex::For( - np, [=] AMREX_GPU_DEVICE (int ip) noexcept - { - if(is_ionized[ip]){ - int i = is_ionized_cumsum[ip]; - WarpXParticleContainer::ParticleType& prod_p = prod_pp[i]; - WarpXParticleContainer::ParticleType& p = pp[ip]; - prod_p.id() = 9999; - prod_p.cpu() = cpuid; - prod_p.pos(0) = p.pos(0); - prod_p.pos(1) = p.pos(1); -#if (AMREX_SPACEDIM == 3) - prod_p.pos(2) = p.pos(2); -#endif - for (int ia = 0; ia < PIdx::nattribs; ++ia) { - prod_pa[ia][i] = pa[ia][ip]; - } - } - } - ); + createIonizedParticles(lev, mfi, pc_source, pc_product, p_is_ionized); } // MFIter - } - } + } // lev + } // pc_source } diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index 6f9fdc485..8dab1a2f7 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -111,6 +111,9 @@ public: virtual int doFieldIonization (const int lev) override; + virtual void buildIonizationMask (const amrex::MFIter& mfi, const int lev, + amrex::Gpu::ManagedVector& 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 fa01f472d..7870e683c 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -2133,3 +2133,85 @@ PhysicalParticleContainer::doFieldIonization(const int lev) Print()<<"in PhysicalParticleContainer::doFieldIonization\n"; return 0; } + +/* \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::ManagedVector& ionization_mask) +{ + // 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(); + Real* ilev_real = soa.GetRealData(particle_comps["ionization_level"]).data(); + + Real c = PhysConst::c; + Real c2_inv = 1./c/c; + + // Loop over all particles in grid/tile. If ionized, set mask to 1 + // and increment ionization level. + ParallelFor( + np, + [=] AMREX_GPU_DEVICE (long i) { + 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] ) + ); + // Get index of ionization_level + int ilev = (int) round(ilev_real[i]); + // Compute probability of ionization p + Real p; + Real w_dtau = 1./ ga * p_adk_prefactor[ilev] * + std::pow(E,p_adk_power[ilev]) * + std::exp( p_adk_exp_prefactor[ilev]/E ); + p = 1. - std::exp( - w_dtau ); + p_ionization_mask[i] = 0; + + if (random_draw < p){ + // increment particle's ionization level + ilev_real[i] += 1.; + // update mask + p_ionization_mask[i] = 1; + } + } + ); +} diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index 1c1b3fd22..017a857a7 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -250,6 +250,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) @@ -277,6 +279,12 @@ public: virtual int doFieldIonization ( const int lev) { return 0; }; + virtual void buildIonizationMask (const amrex::MFIter& mfi, const int lev, + amrex::Gpu::ManagedVector& ionization_mask) + {}; + + std::map getParticleComps() { return particle_comps;} + protected: std::map particle_comps; -- cgit v1.2.3 From 2245a69b89b0d3f632bf3ee4efcc1267308903d7 Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Tue, 6 Aug 2019 13:51:50 -0700 Subject: escape if no ionized particle --- Source/Particles/MultiParticleContainer.cpp | 30 +++++++++++++++++++++++++++-- 1 file changed, 28 insertions(+), 2 deletions(-) (limited to 'Source/Particles/MultiParticleContainer.cpp') diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 6c6a03231..84cd107f8 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -521,14 +521,16 @@ namespace std::unique_ptr< WarpXParticleContainer>& pc_product, const int * const p_is_ionized) { - + // Print()<<"in createIonizedParticles\n"; 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 + //Print()<<"la 1\n"; WarpXParticleContainer::ParticleType* particles_source = ptile_source.GetArrayOfStructs()().data(); // --- source SoA particle data auto& soa_source = ptile_source.GetStructOfArrays(); @@ -543,7 +545,9 @@ namespace && pc_source->DoBoostedFrameDiags(); GpuArray runtime_attribs_source; // + //Print()<<"la 2\n"; if (do_boosted_product && do_boosted_source) { + //Print()<<"do_boosted_product && do_boosted_source\n"; // If boosted frame diagnostics for source species, store them std::map comps_source = pc_source->getParticleComps(); runtime_attribs_source[0] = soa_source.GetRealData(comps_source[ "xold"]).data(); @@ -553,6 +557,7 @@ namespace runtime_attribs_source[4] = soa_source.GetRealData(comps_source["uyold"]).data(); runtime_attribs_source[5] = soa_source.GetRealData(comps_source["uzold"]).data(); } else if (do_boosted_product && !do_boosted_source){ + //Print()<<"do_boosted_product && !do_boosted_source\n"; // Otherwise, store current particle momenta. // Positions are copied from AoS data. runtime_attribs_source[3] = soa_source.GetRealData(PIdx::ux).data(); @@ -563,20 +568,27 @@ namespace // Indices of product particle for each ionized source particle. // i_product[i] is the location in product tile of product particle // from source particle i. + //Print()<<"la 3\n"; amrex::Gpu::ManagedVector i_product; i_product.resize(np_source); // 0GetParticles(lev)[std::make_pair(grid_id,tile_id)]; // old and new (i.e., including ionized particles) number of particles @@ -589,15 +601,18 @@ namespace // First element is the first newly-created product particle WarpXParticleContainer::ParticleType* particles_product = ptile_product.GetArrayOfStructs()().data() + np_product_old; // --- product SoA particle data + //Print()<<"la 8\n"; auto& soa_product = ptile_product.GetStructOfArrays(); GpuArray 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; } + //Print()<<"la 9\n"; // --- product runtime attribs GpuArray runtime_attribs_product; if (do_boosted_product) { + //Print()<<"do_boosted_product 1\n"; std::map 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; @@ -645,18 +660,23 @@ namespace // handling depending on do_boosted_source. For momentum, // runtime_attribs_source[3-5] contains appropriate data. if (do_boosted_product) { + //Print()<<"do_boosted_product\n"; if (do_boosted_source) { + //Print()<<"do_boosted_product do_boosted_source\n"; runtime_attribs_product[0][ip] = runtime_attribs_source[0][ip]; runtime_attribs_product[1][ip] = runtime_attribs_source[1][ip]; runtime_attribs_product[2][ip] = runtime_attribs_source[2][ip]; } else { + //Print()<<"do_boosted_product NO do_boosted_source\n"; 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); } + //Print()<<"momenta\n"; runtime_attribs_product[3][ip] = runtime_attribs_source[3][ip]; runtime_attribs_product[4][ip] = runtime_attribs_source[4][ip]; runtime_attribs_product[5][ip] = runtime_attribs_source[5][ip]; + //Print()<<"done\n"; } } } @@ -668,6 +688,8 @@ void MultiParticleContainer::doFieldIonization() { + Print()<<"in MultiParticleContainer::doFieldIonization\n"; + // Loop over all species. // Ionized particles in pc_source create particles in pc_product for (auto& pc_source : allcontainers){ @@ -718,11 +740,15 @@ MultiParticleContainer::doFieldIonization() { // Ionization mask: one element per source particles. // 0 if not ionized, 1 if ionized. + // Print()<<"here 1\n"; amrex::Gpu::ManagedVector is_ionized; + // Print()<<"here 2\n"; pc_source->buildIonizationMask(mfi, lev, is_ionized); + // Print()<<"here 3\n"; const int * const AMREX_RESTRICT p_is_ionized = is_ionized.dataPtr(); - + // Print()<<"here 4\n"; createIonizedParticles(lev, mfi, pc_source, pc_product, p_is_ionized); + // Print()<<"here 5\n"; } // MFIter } // lev } // pc_source -- cgit v1.2.3 From dab4a0081eac6929db26d16b79ef57c5882df0b4 Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Wed, 7 Aug 2019 15:06:49 -0700 Subject: use c++ deposition by default. includes old attribs --- Source/Particles/MultiParticleContainer.cpp | 28 +++----------- Source/Particles/PhysicalParticleContainer.cpp | 53 ++++++++++++-------------- Source/WarpX.cpp | 2 +- 3 files changed, 31 insertions(+), 52 deletions(-) (limited to 'Source/Particles/MultiParticleContainer.cpp') diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 84cd107f8..36b41e0aa 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -516,12 +516,11 @@ MultiParticleContainer::getSpeciesID(std::string product_str) namespace { static void createIonizedParticles( - int lev, const MFIter& mfi, + int lev, const MFIter& mfi, std::unique_ptr< WarpXParticleContainer>& pc_source, std::unique_ptr< WarpXParticleContainer>& pc_product, const int * const p_is_ionized) { - // Print()<<"in createIonizedParticles\n"; const int grid_id = mfi.index(); const int tile_id = mfi.LocalTileIndex(); @@ -530,7 +529,6 @@ namespace const int np_source = ptile_source.GetArrayOfStructs().size(); if (np_source == 0) return; // --- source AoS particle data - //Print()<<"la 1\n"; WarpXParticleContainer::ParticleType* particles_source = ptile_source.GetArrayOfStructs()().data(); // --- source SoA particle data auto& soa_source = ptile_source.GetStructOfArrays(); @@ -544,10 +542,10 @@ namespace bool do_boosted_source = WarpX::do_boosted_frame_diagnostic && pc_source->DoBoostedFrameDiags(); GpuArray runtime_attribs_source; - // - //Print()<<"la 2\n"; + // Prepare arrays for boosted frame diagnostics. + // If do_boosted_product, need different treatment + // depending on do_boosted_source if (do_boosted_product && do_boosted_source) { - //Print()<<"do_boosted_product && do_boosted_source\n"; // If boosted frame diagnostics for source species, store them std::map comps_source = pc_source->getParticleComps(); runtime_attribs_source[0] = soa_source.GetRealData(comps_source[ "xold"]).data(); @@ -557,7 +555,6 @@ namespace runtime_attribs_source[4] = soa_source.GetRealData(comps_source["uyold"]).data(); runtime_attribs_source[5] = soa_source.GetRealData(comps_source["uzold"]).data(); } else if (do_boosted_product && !do_boosted_source){ - //Print()<<"do_boosted_product && !do_boosted_source\n"; // Otherwise, store current particle momenta. // Positions are copied from AoS data. runtime_attribs_source[3] = soa_source.GetRealData(PIdx::ux).data(); @@ -568,27 +565,20 @@ namespace // Indices of product particle for each ionized source particle. // i_product[i] is the location in product tile of product particle // from source particle i. - //Print()<<"la 3\n"; amrex::Gpu::ManagedVector i_product; i_product.resize(np_source); // 0GetParticles(lev)[std::make_pair(grid_id,tile_id)]; // old and new (i.e., including ionized particles) number of particles @@ -601,18 +591,15 @@ namespace // First element is the first newly-created product particle WarpXParticleContainer::ParticleType* particles_product = ptile_product.GetArrayOfStructs()().data() + np_product_old; // --- product SoA particle data - //Print()<<"la 8\n"; auto& soa_product = ptile_product.GetStructOfArrays(); GpuArray 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; } - //Print()<<"la 9\n"; // --- product runtime attribs GpuArray runtime_attribs_product; if (do_boosted_product) { - //Print()<<"do_boosted_product 1\n"; std::map 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; @@ -640,7 +627,7 @@ namespace if(p_is_ionized[is]){ int ip = p_i_product[is]; // is: index of ionized particle in source species - // ip: index of corresponding created particle in product 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 @@ -660,23 +647,18 @@ namespace // handling depending on do_boosted_source. For momentum, // runtime_attribs_source[3-5] contains appropriate data. if (do_boosted_product) { - //Print()<<"do_boosted_product\n"; if (do_boosted_source) { - //Print()<<"do_boosted_product do_boosted_source\n"; runtime_attribs_product[0][ip] = runtime_attribs_source[0][ip]; runtime_attribs_product[1][ip] = runtime_attribs_source[1][ip]; runtime_attribs_product[2][ip] = runtime_attribs_source[2][ip]; } else { - //Print()<<"do_boosted_product NO do_boosted_source\n"; 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); } - //Print()<<"momenta\n"; runtime_attribs_product[3][ip] = runtime_attribs_source[3][ip]; runtime_attribs_product[4][ip] = runtime_attribs_source[4][ip]; runtime_attribs_product[5][ip] = runtime_attribs_source[5][ip]; - //Print()<<"done\n"; } } } diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 7870e683c..75e975131 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -2094,14 +2094,11 @@ void PhysicalParticleContainer::InitIonizationModule() int offset = ion_energy_offsets[ion_element_id]; for(int i=0; i WarpX::boost_direction = {0,0,0}; int WarpX::do_compute_max_step_from_zmax = 0; Real WarpX::zmax_plasma_to_compute_max_step = 0.; -long WarpX::use_picsar_deposition = 1; +long WarpX::use_picsar_deposition = 0; long WarpX::current_deposition_algo; long WarpX::charge_deposition_algo; long WarpX::field_gathering_algo; -- cgit v1.2.3 From f2d1b0f716f0b94612ce30d7c95c418a928c66b8 Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Wed, 7 Aug 2019 15:13:56 -0700 Subject: remove copy of old attribs, use current values instead for simplicity --- Source/Particles/MultiParticleContainer.cpp | 52 ++++++++--------------------- 1 file changed, 14 insertions(+), 38 deletions(-) (limited to 'Source/Particles/MultiParticleContainer.cpp') diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 36b41e0aa..136d23f6c 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -537,30 +537,11 @@ namespace attribs_source[ia] = soa_source.GetRealData(ia).data(); } // --- source runtime attribs - bool do_boosted_product = WarpX::do_boosted_frame_diagnostic - && pc_product->DoBoostedFrameDiags(); - bool do_boosted_source = WarpX::do_boosted_frame_diagnostic - && pc_source->DoBoostedFrameDiags(); - GpuArray runtime_attribs_source; + GpuArray runtime_uold_source; // Prepare arrays for boosted frame diagnostics. - // If do_boosted_product, need different treatment - // depending on do_boosted_source - if (do_boosted_product && do_boosted_source) { - // If boosted frame diagnostics for source species, store them - std::map comps_source = pc_source->getParticleComps(); - runtime_attribs_source[0] = soa_source.GetRealData(comps_source[ "xold"]).data(); - runtime_attribs_source[1] = soa_source.GetRealData(comps_source[ "yold"]).data(); - runtime_attribs_source[2] = soa_source.GetRealData(comps_source[ "zold"]).data(); - runtime_attribs_source[3] = soa_source.GetRealData(comps_source["uxold"]).data(); - runtime_attribs_source[4] = soa_source.GetRealData(comps_source["uyold"]).data(); - runtime_attribs_source[5] = soa_source.GetRealData(comps_source["uzold"]).data(); - } else if (do_boosted_product && !do_boosted_source){ - // Otherwise, store current particle momenta. - // Positions are copied from AoS data. - runtime_attribs_source[3] = soa_source.GetRealData(PIdx::ux).data(); - runtime_attribs_source[4] = soa_source.GetRealData(PIdx::uy).data(); - runtime_attribs_source[5] = soa_source.GetRealData(PIdx::uz).data(); - } + 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] is the location in product tile of product particle @@ -599,6 +580,8 @@ namespace } // --- product runtime attribs GpuArray runtime_attribs_product; + bool do_boosted_product = WarpX::do_boosted_frame_diagnostic + && pc_product->DoBoostedFrameDiags(); if (do_boosted_product) { std::map comps_product = pc_product->getParticleComps(); runtime_attribs_product[0] = soa_product.GetRealData(comps_product[ "xold"]).data() + np_product_old; @@ -643,22 +626,15 @@ namespace attribs_product[ia][ip] = attribs_source[ia][is]; } // Update xold etc. if boosted frame diagnostics required - // for product species. For position, we need a different - // handling depending on do_boosted_source. For momentum, - // runtime_attribs_source[3-5] contains appropriate data. + // for product species. Fill runtime attribs with a copy of + // current properties (xold = x etc.). if (do_boosted_product) { - if (do_boosted_source) { - runtime_attribs_product[0][ip] = runtime_attribs_source[0][ip]; - runtime_attribs_product[1][ip] = runtime_attribs_source[1][ip]; - runtime_attribs_product[2][ip] = runtime_attribs_source[2][ip]; - } else { - 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_attribs_source[3][ip]; - runtime_attribs_product[4][ip] = runtime_attribs_source[4][ip]; - runtime_attribs_product[5][ip] = runtime_attribs_source[5][ip]; + 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]; } } } -- cgit v1.2.3 From 4be45ed13bcd1dfab7455e09b2a9b96cd2064fca Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Wed, 7 Aug 2019 15:18:36 -0700 Subject: cleaning, remote outdated functions --- Source/Particles/MultiParticleContainer.cpp | 9 +-------- Source/Particles/PhysicalParticleContainer.H | 2 -- Source/Particles/PhysicalParticleContainer.cpp | 13 ------------- Source/Particles/WarpXParticleContainer.H | 3 --- 4 files changed, 1 insertion(+), 26 deletions(-) (limited to 'Source/Particles/MultiParticleContainer.cpp') diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 136d23f6c..48f685e8f 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -646,8 +646,6 @@ void MultiParticleContainer::doFieldIonization() { - Print()<<"in MultiParticleContainer::doFieldIonization\n"; - // Loop over all species. // Ionized particles in pc_source create particles in pc_product for (auto& pc_source : allcontainers){ @@ -698,16 +696,11 @@ MultiParticleContainer::doFieldIonization() { // Ionization mask: one element per source particles. // 0 if not ionized, 1 if ionized. - // Print()<<"here 1\n"; amrex::Gpu::ManagedVector is_ionized; - // Print()<<"here 2\n"; pc_source->buildIonizationMask(mfi, lev, is_ionized); - // Print()<<"here 3\n"; const int * const AMREX_RESTRICT p_is_ionized = is_ionized.dataPtr(); - // Print()<<"here 4\n"; createIonizedParticles(lev, mfi, pc_source, pc_product, p_is_ionized); - // Print()<<"here 5\n"; - } // MFIter + } } // lev } // pc_source } diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index 8dab1a2f7..c10dcc61b 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -108,8 +108,6 @@ public: virtual void PostRestart () final {} void SplitParticles(int lev); - - virtual int doFieldIonization (const int lev) override; virtual void buildIonizationMask (const amrex::MFIter& mfi, const int lev, amrex::Gpu::ManagedVector& ionization_mask) override; diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 75e975131..02602b844 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -80,12 +80,6 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core) : WarpXParticleContainer(amr_core, 0) { plasma_injector.reset(new PlasmaInjector()); -/* - if (do_field_ionization){ - Print()<<"AddRealComp\n"; - AddRealComp("ionization_level"); - } -*/ } void PhysicalParticleContainer::InitData() @@ -2121,13 +2115,6 @@ void PhysicalParticleContainer::InitIonizationModule() } } -int -PhysicalParticleContainer::doFieldIonization(const int lev) -{ - Print()<<"in PhysicalParticleContainer::doFieldIonization\n"; - return 0; -} - /* \brief create mask of ionized particles (1 if ionized, 0 otherwise) * * \param mfi: tile or grid diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index 017a857a7..82190b307 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -276,9 +276,6 @@ public: virtual void setIonizationProduct (int i_product) {}; - virtual int doFieldIonization ( - const int lev) { return 0; }; - virtual void buildIonizationMask (const amrex::MFIter& mfi, const int lev, amrex::Gpu::ManagedVector& ionization_mask) {}; -- cgit v1.2.3 From 30fa1002b9d106acbbbbfd0f01c7e80f47c72747 Mon Sep 17 00:00:00 2001 From: Maxence Thevenet Date: Wed, 7 Aug 2019 18:57:59 -0400 Subject: add profile data, remove typo in example --- Examples/Modules/ionization/inputs.2d | 1 - Source/Particles/MultiParticleContainer.cpp | 4 +++- Source/Particles/PhysicalParticleContainer.cpp | 1 + 3 files changed, 4 insertions(+), 2 deletions(-) (limited to 'Source/Particles/MultiParticleContainer.cpp') diff --git a/Examples/Modules/ionization/inputs.2d b/Examples/Modules/ionization/inputs.2d index 1f2b77a7e..ddf30763a 100644 --- a/Examples/Modules/ionization/inputs.2d +++ b/Examples/Modules/ionization/inputs.2d @@ -4,7 +4,6 @@ max_step = 500 warpx.verbose = 1 warpx.cfl = 1 -warpx.plot_rho = 1 warpx.do_pml = 1 warpx.use_filter = 1 diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 48f685e8f..dafe8fe93 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -521,6 +521,8 @@ namespace std::unique_ptr< WarpXParticleContainer>& pc_product, const int * const p_is_ionized) { + BL_PROFILE("createIonizedParticles"); + const int grid_id = mfi.index(); const int tile_id = mfi.LocalTileIndex(); @@ -645,7 +647,7 @@ namespace 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){ diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 02602b844..22b7df88f 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -2127,6 +2127,7 @@ void PhysicalParticleContainer::buildIonizationMask(const amrex::MFIter& mfi, const int lev, amrex::Gpu::ManagedVector& 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(); -- cgit v1.2.3 From a27dbd90b54e81652abe789a50b913c1bf9831fc Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Wed, 7 Aug 2019 17:49:33 -0700 Subject: use amrex::inclusive_scan for portability --- Source/Particles/MultiParticleContainer.cpp | 26 ++++++++++++++++---------- 1 file changed, 16 insertions(+), 10 deletions(-) (limited to 'Source/Particles/MultiParticleContainer.cpp') diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 48f685e8f..b83478e62 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -519,8 +519,10 @@ namespace int lev, const MFIter& mfi, std::unique_ptr< WarpXParticleContainer>& pc_source, std::unique_ptr< WarpXParticleContainer>& pc_product, - const int * const p_is_ionized) + amrex::Gpu::ManagedVector& is_ionized) { + const int * const AMREX_RESTRICT p_is_ionized = is_ionized.dataPtr(); + const int grid_id = mfi.index(); const int tile_id = mfi.LocalTileIndex(); @@ -544,17 +546,20 @@ namespace runtime_uold_source[2] = soa_source.GetRealData(PIdx::uz).data(); // Indices of product particle for each ionized source particle. - // i_product[i] is the location in product tile of product particle + // i_product[i]-1 is the location in product tile of product particle // from source particle i. amrex::Gpu::ManagedVector i_product; i_product.resize(np_source); // 0 is_ionized; pc_source->buildIonizationMask(mfi, lev, is_ionized); - const int * const AMREX_RESTRICT p_is_ionized = is_ionized.dataPtr(); - createIonizedParticles(lev, mfi, pc_source, pc_product, p_is_ionized); + // Create particles in pc_product + createIonizedParticles(lev, mfi, pc_source, pc_product, is_ionized); } } // lev } // pc_source -- cgit v1.2.3 From b2754b63e9a25c550cc849f84078503bc94a417f Mon Sep 17 00:00:00 2001 From: Maxence Thevenet Date: Wed, 7 Aug 2019 21:17:42 -0400 Subject: exclusive_scan requires ManagedDeviceVector --- Source/Particles/MultiParticleContainer.cpp | 6 +++--- Source/Particles/PhysicalParticleContainer.H | 2 +- Source/Particles/PhysicalParticleContainer.cpp | 2 +- Source/Particles/WarpXParticleContainer.H | 2 +- 4 files changed, 6 insertions(+), 6 deletions(-) (limited to 'Source/Particles/MultiParticleContainer.cpp') diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 9eb978d05..6d977e16c 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -519,7 +519,7 @@ namespace int lev, const MFIter& mfi, std::unique_ptr< WarpXParticleContainer>& pc_source, std::unique_ptr< WarpXParticleContainer>& pc_product, - amrex::Gpu::ManagedVector& is_ionized) + amrex::Gpu::ManagedDeviceVector& is_ionized) { BL_PROFILE("createIonizedParticles"); @@ -550,7 +550,7 @@ namespace // 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::ManagedVector i_product; + amrex::Gpu::ManagedDeviceVector i_product; i_product.resize(np_source); // 0 is_ionized; + amrex::Gpu::ManagedDeviceVector is_ionized; pc_source->buildIonizationMask(mfi, lev, is_ionized); // Create particles in pc_product createIonizedParticles(lev, mfi, pc_source, pc_product, is_ionized); diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index c10dcc61b..44c5030fa 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -110,7 +110,7 @@ public: void SplitParticles(int lev); virtual void buildIonizationMask (const amrex::MFIter& mfi, const int lev, - amrex::Gpu::ManagedVector& ionization_mask) override; + amrex::Gpu::ManagedDeviceVector& 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 22b7df88f..f1544d90c 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -2125,7 +2125,7 @@ void PhysicalParticleContainer::InitIonizationModule() */ void PhysicalParticleContainer::buildIonizationMask(const amrex::MFIter& mfi, const int lev, - amrex::Gpu::ManagedVector& ionization_mask) + amrex::Gpu::ManagedDeviceVector& ionization_mask) { BL_PROFILE("PPC::buildIonizationMask"); // Get pointers to ionization data from pc_source diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index 82190b307..334eb3282 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -277,7 +277,7 @@ public: virtual void setIonizationProduct (int i_product) {}; virtual void buildIonizationMask (const amrex::MFIter& mfi, const int lev, - amrex::Gpu::ManagedVector& ionization_mask) + amrex::Gpu::ManagedDeviceVector& ionization_mask) {}; std::map getParticleComps() { return particle_comps;} -- cgit v1.2.3 From 8dcbb28a52f429e22a6643ffa33e81c6224942e0 Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Thu, 8 Aug 2019 11:31:22 -0700 Subject: add some const, and remove unnecessary includes --- Source/Particles/MultiParticleContainer.H | 1 - Source/Particles/MultiParticleContainer.cpp | 5 ----- Source/Particles/WarpXParticleContainer.H | 2 +- Source/Particles/WarpXParticleContainer.cpp | 2 +- 4 files changed, 2 insertions(+), 8 deletions(-) (limited to 'Source/Particles/MultiParticleContainer.cpp') diff --git a/Source/Particles/MultiParticleContainer.H b/Source/Particles/MultiParticleContainer.H index 8f7fd56e3..fd9a33e86 100644 --- a/Source/Particles/MultiParticleContainer.H +++ b/Source/Particles/MultiParticleContainer.H @@ -7,7 +7,6 @@ #include #include -#include #include #include #include diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 6d977e16c..bff369edc 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -1,15 +1,10 @@ #include #include #include -#include -#include -#include -#include #include #include #include -#include using namespace amrex; diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index 2613be167..dd5a150c4 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -184,7 +184,7 @@ public: RealVector& uxp, RealVector& uyp, RealVector& uzp, - int* ion_lev, + const int * const ion_lev, amrex::MultiFab* jx, amrex::MultiFab* jy, amrex::MultiFab* jz, diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index a7c30764d..eba072572 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -428,7 +428,7 @@ void WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, RealVector& wp, RealVector& uxp, RealVector& uyp, RealVector& uzp, - int* ion_lev, + 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, -- cgit v1.2.3 From e3bb554cbc8df012084d652bc581690d87655add Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Thu, 8 Aug 2019 11:39:54 -0700 Subject: catch errors in inputs --- Source/Particles/MultiParticleContainer.cpp | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) (limited to 'Source/Particles/MultiParticleContainer.cpp') diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index bff369edc..215bc5bd0 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -479,6 +479,9 @@ MultiParticleContainer::mapSpeciesProduct() // 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; } } @@ -500,14 +503,12 @@ MultiParticleContainer::getSpeciesID(std::string product_str) i_product = i; } } - AMREX_ALWAYS_ASSERT_WITH_MESSAGE(found != 0, - "ERROR: could not find ID for species"); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + found != 0, + "ERROR: could not find product species ID for ionization. Wrong name?"); return i_product; } -//void -//MultiParticleContainer::InitIonizationModule() - namespace { static void createIonizedParticles( -- cgit v1.2.3 From 33c7e5860405f7785aec896200bda495d9c5b4a4 Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Thu, 8 Aug 2019 13:03:21 -0700 Subject: add few comments --- Source/Particles/MultiParticleContainer.cpp | 2 ++ 1 file changed, 2 insertions(+) (limited to 'Source/Particles/MultiParticleContainer.cpp') diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 215bc5bd0..24f9f2f11 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -511,6 +511,8 @@ MultiParticleContainer::getSpeciesID(std::string product_str) namespace { + // For particle i in mfi, if is_ionized[i]=1, copy particle from + // particle i container pc_source into pc_product static void createIonizedParticles( int lev, const MFIter& mfi, std::unique_ptr< WarpXParticleContainer>& pc_source, -- cgit v1.2.3 From 74f85ac51bea0f1c53d78d6a9ca87a3534696162 Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Thu, 8 Aug 2019 14:34:40 -0700 Subject: use Weiqun's convention: space before brackets for function declaration and definition --- Source/Particles/MultiParticleContainer.H | 2 +- Source/Particles/MultiParticleContainer.cpp | 10 +++++----- Source/Particles/PhysicalParticleContainer.H | 4 +--- Source/Particles/PhysicalParticleContainer.cpp | 12 +++++------- Source/Particles/WarpXParticleContainer.H | 2 +- 5 files changed, 13 insertions(+), 17 deletions(-) (limited to 'Source/Particles/MultiParticleContainer.cpp') diff --git a/Source/Particles/MultiParticleContainer.H b/Source/Particles/MultiParticleContainer.H index fd9a33e86..aaed96423 100644 --- a/Source/Particles/MultiParticleContainer.H +++ b/Source/Particles/MultiParticleContainer.H @@ -127,7 +127,7 @@ public: /// std::unique_ptr GetChargeDensity(int lev, bool local = false); - void doFieldIonization(); + void doFieldIonization (); void Checkpoint (const std::string& dir) const; diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 24f9f2f11..04f052086 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -453,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& pc_source, std::unique_ptr< WarpXParticleContainer>& pc_product, @@ -647,7 +647,7 @@ namespace } void -MultiParticleContainer::doFieldIonization() +MultiParticleContainer::doFieldIonization () { BL_PROFILE("MPC::doFieldIonization"); // Loop over all species. diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index 44c5030fa..c7494fbdf 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -22,9 +22,7 @@ public: virtual void InitData () override; - void InitIonizationLevel (); - - void InitIonizationModule(); + void InitIonizationModule (); #ifdef WARPX_DO_ELECTROSTATIC virtual void FieldGatherES(const amrex::Vector, 3> >& E, diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index a9b57aef6..b20ea8090 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -39,10 +39,6 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp pp.query("do_boosted_frame_diags", do_boosted_frame_diags); pp.query("do_field_ionization", do_field_ionization); - // If do_field_ionization, read initialization data from input file and - // read ionization energies from table. - //if (do_field_ionization) - // InitIonizationModule(); pp.query("plot_species", plot_species); int do_user_plot_vars; @@ -84,6 +80,8 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core) void PhysicalParticleContainer::InitData() { + // Init ionization module here instead of in the PhysicalParticleContainer + // constructor to have access to dt if (do_field_ionization) {InitIonizationModule();} AddParticles(0); // Note - add on level 0 Redistribute(); // We then redistribute @@ -1964,7 +1962,7 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti, } -void PhysicalParticleContainer::InitIonizationModule() +void PhysicalParticleContainer::InitIonizationModule () { if (!do_field_ionization) return; ParmParse pp(species_name); @@ -2017,8 +2015,8 @@ void PhysicalParticleContainer::InitIonizationModule() * depending on whether the particle is ionized or not. */ void -PhysicalParticleContainer::buildIonizationMask(const amrex::MFIter& mfi, const int lev, - amrex::Gpu::ManagedDeviceVector& ionization_mask) +PhysicalParticleContainer::buildIonizationMask (const amrex::MFIter& mfi, const int lev, + amrex::Gpu::ManagedDeviceVector& ionization_mask) { BL_PROFILE("PPC::buildIonizationMask"); // Get pointers to ionization data from pc_source diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index dd5a150c4..efa299c43 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -284,7 +284,7 @@ public: amrex::Gpu::ManagedDeviceVector& ionization_mask) {}; - std::map getParticleComps() { return particle_comps;} + std::map getParticleComps () { return particle_comps;} protected: -- cgit v1.2.3 From c87ee4b5e5d62f8e90a719a42ff164f99983a62c Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Thu, 8 Aug 2019 14:44:06 -0700 Subject: final minor cleaning --- Source/Particles/MultiParticleContainer.cpp | 6 +++--- Source/Particles/PhysicalParticleContainer.cpp | 6 +++--- Source/Utils/WarpXConst.H | 1 - 3 files changed, 6 insertions(+), 7 deletions(-) (limited to 'Source/Particles/MultiParticleContainer.cpp') diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 04f052086..e39930216 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -511,8 +511,8 @@ MultiParticleContainer::getSpeciesID (std::string product_str) namespace { - // For particle i in mfi, if is_ionized[i]=1, copy particle from - // particle i container pc_source into pc_product + // For particle i in mfi, if is_ionized[i]=1, copy particle + // particle i from container pc_source into pc_product static void createIonizedParticles ( int lev, const MFIter& mfi, std::unique_ptr< WarpXParticleContainer>& pc_source, @@ -663,7 +663,7 @@ MultiParticleContainer::doFieldIonization () for (int lev = 0; lev <= pc_source->finestLevel(); ++lev){ #ifdef _OPENMP - // Touch all tiles of source species in serial if additional arguments + // 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(); diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index b20ea8090..466446157 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -81,7 +81,7 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core) void PhysicalParticleContainer::InitData() { // Init ionization module here instead of in the PhysicalParticleContainer - // constructor to have access to dt + // constructor because dt is required if (do_field_ionization) {InitIonizationModule();} AddParticles(0); // Note - add on level 0 Redistribute(); // We then redistribute @@ -208,10 +208,10 @@ PhysicalParticleContainer::CheckAndAddParticle(Real x, Real y, Real z, auto& particle_tile = DefineAndReturnParticleTile(0, 0, 0); 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["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]); diff --git a/Source/Utils/WarpXConst.H b/Source/Utils/WarpXConst.H index bfb67b3d5..6901d4629 100644 --- a/Source/Utils/WarpXConst.H +++ b/Source/Utils/WarpXConst.H @@ -21,7 +21,6 @@ namespace PhysConst static constexpr amrex::Real hbar = 1.05457180010e-34; static constexpr amrex::Real alpha = mu0/(4*MathConst::pi)*q_e*q_e*c/hbar; static constexpr amrex::Real r_e = 1./(4*MathConst::pi*ep0) * q_e*q_e/(m_e*c*c); - // static constexpr amrex::Real alpha = 0.0072973525693; } #endif -- cgit v1.2.3 From 9795e23c12b3f3553c748e2e335b716f3df62e22 Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Wed, 21 Aug 2019 08:45:25 -0700 Subject: add comments and make sure charge=q_e for ionizable species --- Source/Particles/Deposition/CurrentDeposition.H | 12 ++++++++++++ Source/Particles/MultiParticleContainer.cpp | 9 ++++++++- Source/Particles/PhysicalParticleContainer.cpp | 9 +++++++-- Source/Utils/write_atomic_data_cpp.py | 6 ++++++ 4 files changed, 33 insertions(+), 3 deletions(-) (limited to 'Source/Particles/MultiParticleContainer.cpp') diff --git a/Source/Particles/Deposition/CurrentDeposition.H b/Source/Particles/Deposition/CurrentDeposition.H index 86cf89cd3..7dde90e96 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. @@ -37,6 +41,8 @@ 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]; @@ -166,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. @@ -195,6 +205,8 @@ 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; diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index e39930216..9e2cd097f 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -513,7 +513,7 @@ namespace { // For particle i in mfi, if is_ionized[i]=1, copy particle // particle i from container pc_source into pc_product - static void createIonizedParticles ( + void createIonizedParticles ( int lev, const MFIter& mfi, std::unique_ptr< WarpXParticleContainer>& pc_source, std::unique_ptr< WarpXParticleContainer>& pc_product, @@ -558,6 +558,9 @@ namespace // 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; @@ -662,6 +665,10 @@ MultiParticleContainer::doFieldIonization () 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) { diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 4096d1911..4562b5c0a 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1979,12 +1979,17 @@ 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 Real component for ionization level + // Add runtime integer component for ionization level AddIntComp("ionization_level"); - plot_flags.resize(PIdx::nattribs + 1, 1); + // plot_flags.resize(PIdx::nattribs + 1, 1); // 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]; diff --git a/Source/Utils/write_atomic_data_cpp.py b/Source/Utils/write_atomic_data_cpp.py index 21d61a075..b085d50eb 100644 --- a/Source/Utils/write_atomic_data_cpp.py +++ b/Source/Utils/write_atomic_data_cpp.py @@ -1,3 +1,9 @@ +''' +This python script reads ionization tables in atomic_data.txt (generated from +the NIST website) and extracts ionization levels into C++ file +IonizationEnergiesTable.H, which contains tables + metadata. +''' + import re, os import numpy as np from scipy.constants import e -- cgit v1.2.3