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.cpp | 155 +++++++++++++++++++++++++++- 1 file changed, 153 insertions(+), 2 deletions(-) (limited to 'Source/Particles/MultiParticleContainer.cpp') 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 + } + } +} -- cgit v1.2.3