From 7195cc8eecf81896f3f2cd57eaf6161dd6834f5c Mon Sep 17 00:00:00 2001 From: Dave Grote Date: Thu, 16 May 2019 17:56:07 -0700 Subject: Full implementation of the multimode RZ solver --- 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 6d618c096..f39a2a36d 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -299,7 +299,7 @@ MultiParticleContainer::GetChargeDensity (int lev, bool local) std::unique_ptr rho = allcontainers[0]->GetChargeDensity(lev, true); for (unsigned i = 1, n = allcontainers.size(); i < n; ++i) { std::unique_ptr rhoi = allcontainers[i]->GetChargeDensity(lev, true); - MultiFab::Add(*rho, *rhoi, 0, 0, 1, rho->nGrow()); + MultiFab::Add(*rho, *rhoi, 0, 0, rho->nComp(), rho->nGrow()); } if (!local) { const Geometry& gm = allcontainers[0]->Geom(lev); -- cgit v1.2.3 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/Evolve/WarpXEvolveEM.cpp | 2 +- Source/Particles/MultiParticleContainer.H | 5 ++ Source/Particles/MultiParticleContainer.cpp | 51 ++++++++++++++++++++ Source/Particles/PhysicalParticleContainer.H | 22 +++++++++ Source/Particles/PhysicalParticleContainer.cpp | 66 ++++++++++++++++++++++---- Source/Particles/WarpXParticleContainer.H | 14 +++++- 6 files changed, 149 insertions(+), 11 deletions(-) (limited to 'Source/Particles/MultiParticleContainer.cpp') diff --git a/Source/Evolve/WarpXEvolveEM.cpp b/Source/Evolve/WarpXEvolveEM.cpp index aa180cc95..4bd9d4566 100644 --- a/Source/Evolve/WarpXEvolveEM.cpp +++ b/Source/Evolve/WarpXEvolveEM.cpp @@ -278,7 +278,7 @@ WarpX::OneStep_nosub (Real cur_time) { // Loop over species. For each ionizable species, create particles in // product species. - // mypc->doFieldIonization(); + mypc->doFieldIonization(); // Push particle from x^{n} to x^{n+1} // from p^{n-1/2} to p^{n+1/2} // Deposit current j^{n+1/2} diff --git a/Source/Particles/MultiParticleContainer.H b/Source/Particles/MultiParticleContainer.H index 7c9ede411..4ea9c2763 100644 --- a/Source/Particles/MultiParticleContainer.H +++ b/Source/Particles/MultiParticleContainer.H @@ -128,6 +128,8 @@ public: /// std::unique_ptr GetChargeDensity(int lev, bool local = false); + void doFieldIonization(); + void Checkpoint (const std::string& dir) const; void WritePlotFile (const std::string& dir) const; @@ -207,6 +209,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..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, 3> >& E, const amrex::Vector > > >& masks) override; @@ -105,6 +109,24 @@ public: void SplitParticles(int lev); + virtual void copyParticles(int lev) override; + + virtual int doFieldIonization( + const int lev, + amrex::Cuda::ManagedDeviceVector& x_buf, + amrex::Cuda::ManagedDeviceVector& y_buf, + amrex::Cuda::ManagedDeviceVector& z_buf, + amrex::Cuda::ManagedDeviceVector& ux_buf, + amrex::Cuda::ManagedDeviceVector& uy_buf, + amrex::Cuda::ManagedDeviceVector& uz_buf, + amrex::Cuda::ManagedDeviceVector& ex_buf, + amrex::Cuda::ManagedDeviceVector& ey_buf, + amrex::Cuda::ManagedDeviceVector& ez_buf, + amrex::Cuda::ManagedDeviceVector& bx_buf, + amrex::Cuda::ManagedDeviceVector& by_buf, + amrex::Cuda::ManagedDeviceVector& bz_buf, + amrex::Cuda::ManagedDeviceVector& w_buf) 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 d0db0fbdc..d2fcc93e0 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -6,6 +6,7 @@ #include #include #include +#include #include #include @@ -38,10 +39,10 @@ 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){ - AddRealComp("ionization_level"); - plot_flags.resize(PIdx::nattribs + 1, 1); - } + // 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; @@ -79,10 +80,12 @@ 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() @@ -1604,25 +1607,38 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti, copy_attribs(pti, x, y, z); } + Real* AMREX_RESTRICT ion_lev = nullptr; + if (do_field_ionization){ + ion_lev = pti.GetAttribs(particle_comps["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 { @@ -2058,3 +2074,35 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti, } } } + + +void PhysicalParticleContainer::InitIonizationModule() +{ + ParmParse pp(species_name); + pp.query("species_ionization_level", species_ionization_level); + pp.get("ionization_product", ionization_product_name); + pp.get("physical_element", physical_element); + // Add Real component for ionization level + AddRealComp("ionization_level"); + 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]; + ionization_energies.resize(ion_atomic_number); + int offset = ion_energy_offsets[ion_element_id]; + for(int i=0; i particle_comps; @@ -290,8 +292,18 @@ protected: // This is currently required because continuous injection does not // support all features allowed by direct injection. int do_continuous_injection = 0; - int do_field_ionization = 0; + int do_field_ionization = 0; + int ionization_product; + std::string ionization_product_name; + int ion_atomic_number; + int species_ionization_level = 0; + amrex::Gpu::ManagedVector ionization_energies; + amrex::Gpu::ManagedVector adk_power; + amrex::Gpu::ManagedVector adk_prefactor; + amrex::Gpu::ManagedVector adk_exp_prefactor; + std::string physical_element; + int do_boosted_frame_diags = 1; amrex::Vector local_rho; -- cgit v1.2.3 From a0f579f27969c3f86e54874b5a3c6a732dff9f53 Mon Sep 17 00:00:00 2001 From: MaxThevenet 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 From eb4080b05342edd63c09091f393638c9053b7b19 Mon Sep 17 00:00:00 2001 From: Andrew Myers Date: Wed, 21 Aug 2019 11:31:10 -0700 Subject: reimplement the boosted frame diagnostics without the runtime components --- Source/Diagnostics/ParticleIO.cpp | 11 ---- Source/Evolve/WarpXEvolveEM.cpp | 16 +++--- Source/Particles/MultiParticleContainer.cpp | 20 ------- Source/Particles/PhysicalParticleContainer.cpp | 65 ++++++---------------- .../Particles/RigidInjectedParticleContainer.cpp | 43 +++++--------- Source/Particles/WarpXParticleContainer.H | 2 + Source/Particles/WarpXParticleContainer.cpp | 27 --------- 7 files changed, 41 insertions(+), 143 deletions(-) (limited to 'Source/Particles/MultiParticleContainer.cpp') diff --git a/Source/Diagnostics/ParticleIO.cpp b/Source/Diagnostics/ParticleIO.cpp index f159e5302..4ac5cb6a1 100644 --- a/Source/Diagnostics/ParticleIO.cpp +++ b/Source/Diagnostics/ParticleIO.cpp @@ -102,17 +102,6 @@ MultiParticleContainer::WritePlotFile (const std::string& dir) const real_names.push_back("theta"); #endif - if (WarpX::do_boosted_frame_diagnostic && pc->DoBoostedFrameDiags()) - { - real_names.push_back("xold"); - real_names.push_back("yold"); - real_names.push_back("zold"); - - real_names.push_back("uxold"); - real_names.push_back("uyold"); - real_names.push_back("uzold"); - } - // Convert momentum to SI pc->ConvertUnits(ConvertDirection::WarpX_to_SI); // real_names contains a list of all particle attributes. diff --git a/Source/Evolve/WarpXEvolveEM.cpp b/Source/Evolve/WarpXEvolveEM.cpp index 16b5905d1..ff55724f8 100644 --- a/Source/Evolve/WarpXEvolveEM.cpp +++ b/Source/Evolve/WarpXEvolveEM.cpp @@ -140,6 +140,14 @@ WarpX::EvolveEM (int numsteps) bool do_insitu = ((step+1) >= insitu_start) && (insitu_int > 0) && ((step+1) % insitu_int == 0); + if (do_boosted_frame_diagnostic) { + std::unique_ptr cell_centered_data = nullptr; + if (WarpX::do_boosted_frame_fields) { + cell_centered_data = GetCellCenteredData(); + } + myBFD->writeLabFrameData(cell_centered_data.get(), *mypc, geom[0], cur_time, dt[0]); + } + bool move_j = is_synchronized || to_make_plot || do_insitu; // If is_synchronized we need to shift j too so that next step we can evolve E by dt/2. // We might need to move j because we are going to make a plotfile. @@ -173,14 +181,6 @@ WarpX::EvolveEM (int numsteps) t_new[i] = cur_time; } - if (do_boosted_frame_diagnostic) { - std::unique_ptr cell_centered_data = nullptr; - if (WarpX::do_boosted_frame_fields) { - cell_centered_data = GetCellCenteredData(); - } - myBFD->writeLabFrameData(cell_centered_data.get(), *mypc, geom[0], cur_time, dt[0]); - } - // slice gen // if (to_make_plot || do_insitu || to_make_slice_plot) { diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 982e04e39..4980bc62a 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -42,26 +42,6 @@ MultiParticleContainer::MultiParticleContainer (AmrCore* amr_core) nspecies_boosted_frame_diags += 1; } } - - if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) - { - for (int i = 0; i < nspecies_boosted_frame_diags; ++i) - { - int is = map_species_boosted_frame_diags[i]; - allcontainers[is]->AddRealComp("xold"); - allcontainers[is]->AddRealComp("yold"); - allcontainers[is]->AddRealComp("zold"); - allcontainers[is]->AddRealComp("uxold"); - allcontainers[is]->AddRealComp("uyold"); - allcontainers[is]->AddRealComp("uzold"); - } - pc_tmp->AddRealComp("xold"); - pc_tmp->AddRealComp("yold"); - pc_tmp->AddRealComp("zold"); - pc_tmp->AddRealComp("uxold"); - pc_tmp->AddRealComp("uyold"); - pc_tmp->AddRealComp("uzold"); - } } void diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index d10390204..14807d0d9 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -196,18 +196,6 @@ 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) - { - // 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]); - } // add particle AddOneParticle(0, 0, 0, x, y, z, attribs); } @@ -430,15 +418,6 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) for (int ia = 0; ia < PIdx::nattribs; ++ia) { pa[ia] = soa.GetRealData(ia).data() + old_size; } - GpuArray pb; - if (do_boosted) { - pb[0] = soa.GetRealData(particle_comps[ "xold"]).data() + old_size; - pb[1] = soa.GetRealData(particle_comps[ "yold"]).data() + old_size; - pb[2] = soa.GetRealData(particle_comps[ "zold"]).data() + old_size; - pb[3] = soa.GetRealData(particle_comps["uxold"]).data() + old_size; - pb[4] = soa.GetRealData(particle_comps["uyold"]).data() + old_size; - pb[5] = soa.GetRealData(particle_comps["uzold"]).data() + old_size; - } const GpuArray overlap_corner {AMREX_D_DECL(overlap_realbox.lo(0), @@ -589,15 +568,6 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) pa[PIdx::uy][ip] = u.y; pa[PIdx::uz][ip] = u.z; - if (do_boosted) { - pb[0][ip] = x; - pb[1][ip] = y; - pb[2][ip] = z; - pb[3][ip] = u.x; - pb[4][ip] = u.y; - pb[5][ip] = u.z; - } - #if (AMREX_SPACEDIM == 3) p.pos(0) = x; p.pos(1) = y; @@ -1633,21 +1603,22 @@ PhysicalParticleContainer::PushP (int lev, Real dt, void PhysicalParticleContainer::copy_attribs(WarpXParIter& pti,const Real* xp, const Real* yp, const Real* zp) { - auto& attribs = pti.GetAttribs(); - Real* AMREX_RESTRICT uxp = attribs[PIdx::ux].dataPtr(); Real* AMREX_RESTRICT uyp = attribs[PIdx::uy].dataPtr(); Real* AMREX_RESTRICT uzp = attribs[PIdx::uz].dataPtr(); - Real* AMREX_RESTRICT xpold = pti.GetAttribs(particle_comps["xold"]).dataPtr(); - Real* AMREX_RESTRICT ypold = pti.GetAttribs(particle_comps["yold"]).dataPtr(); - Real* AMREX_RESTRICT zpold = pti.GetAttribs(particle_comps["zold"]).dataPtr(); - Real* AMREX_RESTRICT uxpold = pti.GetAttribs(particle_comps["uxold"]).dataPtr(); - Real* AMREX_RESTRICT uypold = pti.GetAttribs(particle_comps["uyold"]).dataPtr(); - Real* AMREX_RESTRICT uzpold = pti.GetAttribs(particle_comps["uzold"]).dataPtr(); - - const long np = pti.numParticles(); + const auto np = pti.numParticles(); + const auto lev = pti.GetLevel(); + const auto index = pti.GetPairIndex(); + tmp_particle_data.resize(finestLevel()+1); + for (int i = 0; i < 6; ++i) tmp_particle_data[lev][index][i].resize(np); + Real* AMREX_RESTRICT xpold = tmp_particle_data[lev][index][0].dataPtr(); + Real* AMREX_RESTRICT ypold = tmp_particle_data[lev][index][1].dataPtr(); + Real* AMREX_RESTRICT zpold = tmp_particle_data[lev][index][2].dataPtr(); + Real* AMREX_RESTRICT uxpold = tmp_particle_data[lev][index][3].dataPtr(); + Real* AMREX_RESTRICT uypold = tmp_particle_data[lev][index][4].dataPtr(); + Real* AMREX_RESTRICT uzpold = tmp_particle_data[lev][index][5].dataPtr(); ParallelFor( np, [=] AMREX_GPU_DEVICE (long i) { @@ -1733,15 +1704,15 @@ void PhysicalParticleContainer::GetParticleSlice(const int direction, const Real auto& uyp_new = attribs[PIdx::uy ]; auto& uzp_new = attribs[PIdx::uz ]; - auto& xp_old = pti.GetAttribs(particle_comps["xold"]); - auto& yp_old = pti.GetAttribs(particle_comps["yold"]); - auto& zp_old = pti.GetAttribs(particle_comps["zold"]); - auto& uxp_old = pti.GetAttribs(particle_comps["uxold"]); - auto& uyp_old = pti.GetAttribs(particle_comps["uyold"]); - auto& uzp_old = pti.GetAttribs(particle_comps["uzold"]); + auto& xp_old = tmp_particle_data[lev][index][0]; + auto& yp_old = tmp_particle_data[lev][index][1]; + auto& zp_old = tmp_particle_data[lev][index][2]; + auto& uxp_old = tmp_particle_data[lev][index][3]; + auto& uyp_old = tmp_particle_data[lev][index][4]; + auto& uzp_old = tmp_particle_data[lev][index][5]; const long np = pti.numParticles(); - + Real uzfrm = -WarpX::gamma_boost*WarpX::beta_boost*PhysConst::c; Real inv_c2 = 1.0/PhysConst::c/PhysConst::c; diff --git a/Source/Particles/RigidInjectedParticleContainer.cpp b/Source/Particles/RigidInjectedParticleContainer.cpp index 36cb9d224..0c94d1e33 100644 --- a/Source/Particles/RigidInjectedParticleContainer.cpp +++ b/Source/Particles/RigidInjectedParticleContainer.cpp @@ -239,15 +239,13 @@ RigidInjectedParticleContainer::PushPX(WarpXParIter& pti, Real* const AMREX_RESTRICT Bzp = attribs[PIdx::Bz].dataPtr(); if (!done_injecting_lev) { - if (!(WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags)) { - // If the old values are not already saved, create copies here. - xp_save = xp; - yp_save = yp; - zp_save = zp; - uxp_save = uxp; - uyp_save = uyp; - uzp_save = uzp; - } + // If the old values are not already saved, create copies here. + xp_save = xp; + yp_save = yp; + zp_save = zp; + uxp_save = uxp; + uyp_save = uyp; + uzp_save = uzp; // Scale the fields of particles about to cross the injection plane. // This only approximates what should be happening. The particles @@ -275,27 +273,12 @@ RigidInjectedParticleContainer::PushPX(WarpXParIter& pti, if (!done_injecting_lev) { - Real* AMREX_RESTRICT x_save; - Real* AMREX_RESTRICT y_save; - Real* AMREX_RESTRICT z_save; - Real* AMREX_RESTRICT ux_save; - Real* AMREX_RESTRICT uy_save; - Real* AMREX_RESTRICT uz_save; - if (!(WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags)) { - x_save = xp_save.dataPtr(); - y_save = yp_save.dataPtr(); - z_save = zp_save.dataPtr(); - ux_save = uxp_save.dataPtr(); - uy_save = uyp_save.dataPtr(); - uz_save = uzp_save.dataPtr(); - } else { - x_save = pti.GetAttribs(particle_comps["xold"]).dataPtr(); - y_save = pti.GetAttribs(particle_comps["yold"]).dataPtr(); - z_save = pti.GetAttribs(particle_comps["zold"]).dataPtr(); - ux_save = pti.GetAttribs(particle_comps["uxold"]).dataPtr(); - uy_save = pti.GetAttribs(particle_comps["uyold"]).dataPtr(); - uz_save = pti.GetAttribs(particle_comps["uzold"]).dataPtr(); - } + Real* AMREX_RESTRICT x_save = xp_save.dataPtr(); + Real* AMREX_RESTRICT y_save = yp_save.dataPtr(); + Real* AMREX_RESTRICT z_save = zp_save.dataPtr(); + Real* AMREX_RESTRICT ux_save = uxp_save.dataPtr(); + Real* AMREX_RESTRICT uy_save = uyp_save.dataPtr(); + Real* AMREX_RESTRICT uz_save = uzp_save.dataPtr(); // Undo the push for particles not injected yet. // The zp are advanced a fixed amount. diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index ac5b47ada..cbdf0cdbc 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -308,6 +308,8 @@ protected: // list of names of attributes to dump. amrex::Vector plot_vars; + amrex::Vector, std::array, 6> > > tmp_particle_data; + private: virtual void particlePostLocate(ParticleType& p, const amrex::ParticleLocData& pld, const int lev) override; diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index befa5cfed..db0f683c7 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -85,17 +85,6 @@ WarpXParticleContainer::WarpXParticleContainer (AmrCore* amr_core, int ispecies) particle_comps["theta"] = PIdx::theta; #endif - if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) - { - particle_comps["xold"] = PIdx::nattribs; - particle_comps["yold"] = PIdx::nattribs+1; - particle_comps["zold"] = PIdx::nattribs+2; - particle_comps["uxold"] = PIdx::nattribs+3; - particle_comps["uyold"] = PIdx::nattribs+4; - particle_comps["uzold"] = PIdx::nattribs+5; - - } - // Initialize temporary local arrays for charge/current deposition int num_threads = 1; #ifdef _OPENMP @@ -238,14 +227,6 @@ 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]); - } - particle_tile.push_back(p); } @@ -256,14 +237,6 @@ 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); - } - for (int comp = PIdx::uz+1; comp < PIdx::nattribs; ++comp) { #ifdef WARPX_DIM_RZ -- cgit v1.2.3