diff options
Diffstat (limited to 'Source/Particles/PhotonParticleContainer.cpp')
| -rw-r--r-- | Source/Particles/PhotonParticleContainer.cpp | 158 |
1 files changed, 156 insertions, 2 deletions
diff --git a/Source/Particles/PhotonParticleContainer.cpp b/Source/Particles/PhotonParticleContainer.cpp index 3e678fb26..ca4dcc81b 100644 --- a/Source/Particles/PhotonParticleContainer.cpp +++ b/Source/Particles/PhotonParticleContainer.cpp @@ -79,8 +79,6 @@ PhotonParticleContainer::PushPX(WarpXParIter& pti, copy_attribs(pti, x, y, z); } - //No need to update momentum for photons (for now) - amrex::ParallelFor( pti.numParticles(), [=] AMREX_GPU_DEVICE (long i) { @@ -90,6 +88,28 @@ PhotonParticleContainer::PushPX(WarpXParIter& pti, } ); +#ifdef WARPX_QED + DoBreitWheelerPti(pti, dt); +#endif + +} + +void +PhotonParticleContainer::PushP (int lev, + amrex::Real dt, + const amrex::MultiFab& Ex, + const amrex::MultiFab& Ey, + const amrex::MultiFab& Ez, + const amrex::MultiFab& Bx, + const amrex::MultiFab& By, + const amrex::MultiFab& Bz) +{ + BL_PROFILE("PhotonParticleContainer::PushP"); + if (do_not_push) return; + +#ifdef WARPX_QED + if(has_breit_wheeler()) DoBreitWheeler(lev,dt, Ex,Ey,Ez,Bx,By,Bz); +#endif } void @@ -117,3 +137,137 @@ PhotonParticleContainer::Evolve (int lev, t, dt); } + +#ifdef WARPX_QED +void +PhotonParticleContainer::DoBreitWheeler(int lev, + amrex::Real dt, + const amrex::MultiFab& Ex, + const amrex::MultiFab& Ey, + const amrex::MultiFab& Ez, + const amrex::MultiFab& Bx, + const amrex::MultiFab& By, + const amrex::MultiFab& Bz) +{ +#ifdef _OPENMP + #pragma omp parallel +#endif + { +#ifdef _OPENMP + int thread_num = omp_get_thread_num(); +#else + int thread_num = 0; +#endif + + for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) + { + const Box& box = pti.validbox(); + + auto& attribs = pti.GetAttribs(); + + auto& Exp = attribs[PIdx::Ex]; + auto& Eyp = attribs[PIdx::Ey]; + auto& Ezp = attribs[PIdx::Ez]; + auto& Bxp = attribs[PIdx::Bx]; + auto& Byp = attribs[PIdx::By]; + auto& Bzp = attribs[PIdx::Bz]; + + const long np = pti.numParticles(); + + // Data on the grid + const FArrayBox& exfab = Ex[pti]; + const FArrayBox& eyfab = Ey[pti]; + const FArrayBox& ezfab = Ez[pti]; + const FArrayBox& bxfab = Bx[pti]; + const FArrayBox& byfab = By[pti]; + const FArrayBox& bzfab = Bz[pti]; + + Exp.assign(np,WarpX::E_external[0]); + Eyp.assign(np,WarpX::E_external[1]); + Ezp.assign(np,WarpX::E_external[2]); + + Bxp.assign(np,WarpX::B_external[0]); + Byp.assign(np,WarpX::B_external[1]); + Bzp.assign(np,WarpX::B_external[2]); + + // + // copy data from particle container to temp arrays + // + pti.GetPosition(m_xp[thread_num], m_yp[thread_num], m_zp[thread_num]); + + int e_is_nodal = Ex.is_nodal() and Ey.is_nodal() and Ez.is_nodal(); + FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp, + &exfab, &eyfab, &ezfab, &bxfab, &byfab, &bzfab, + Ex.nGrow(), e_is_nodal, + 0, np, thread_num, lev, lev); + + // This wraps the momentum advance so that inheritors can modify the call. + // Extract pointers to the different particle quantities + ParticleReal* const AMREX_RESTRICT ux = attribs[PIdx::ux].dataPtr(); + ParticleReal* const AMREX_RESTRICT uy = attribs[PIdx::uy].dataPtr(); + ParticleReal* const AMREX_RESTRICT uz = attribs[PIdx::uz].dataPtr(); + const ParticleReal* const AMREX_RESTRICT Expp = Exp.dataPtr(); + const ParticleReal* const AMREX_RESTRICT Eypp = Eyp.dataPtr(); + const ParticleReal* const AMREX_RESTRICT Ezpp = Ezp.dataPtr(); + const ParticleReal* const AMREX_RESTRICT Bxpp = Bxp.dataPtr(); + const ParticleReal* const AMREX_RESTRICT Bypp = Byp.dataPtr(); + const ParticleReal* const AMREX_RESTRICT Bzpp = Bzp.dataPtr(); + + BreitWheelerEvolveOpticalDepth evolve_opt = + shr_ptr_bw_engine->build_evolve_functor(); + + amrex::Real* AMREX_RESTRICT p_tau = + pti.GetAttribs(particle_comps["tau"]).dataPtr(); + + const auto me = PhysConst::m_e; + + // Loop over the particles and update their optical_depth + amrex::ParallelFor( + pti.numParticles(), + [=] AMREX_GPU_DEVICE (long i) { + evolve_opt( + me*ux[i], me*uy[i], me*uz[i], + Expp[i], Eypp[i], Ezpp[i], + Bxpp[i], Bypp[i], Bzpp[i], + dt, p_tau[i]); + } + ); + } + } +} + +void +PhotonParticleContainer::DoBreitWheelerPti(WarpXParIter& pti, + amrex::Real dt) +{ + auto& attribs = pti.GetAttribs(); + ParticleReal* const AMREX_RESTRICT ux = attribs[PIdx::ux].dataPtr(); + ParticleReal* const AMREX_RESTRICT uy = attribs[PIdx::uy].dataPtr(); + ParticleReal* const AMREX_RESTRICT uz = attribs[PIdx::uz].dataPtr(); + const ParticleReal* const AMREX_RESTRICT Ex = attribs[PIdx::Ex].dataPtr(); + const ParticleReal* const AMREX_RESTRICT Ey = attribs[PIdx::Ey].dataPtr(); + const ParticleReal* const AMREX_RESTRICT Ez = attribs[PIdx::Ez].dataPtr(); + const ParticleReal* const AMREX_RESTRICT Bx = attribs[PIdx::Bx].dataPtr(); + const ParticleReal* const AMREX_RESTRICT By = attribs[PIdx::By].dataPtr(); + const ParticleReal* const AMREX_RESTRICT Bz = attribs[PIdx::Bz].dataPtr(); + + BreitWheelerEvolveOpticalDepth evolve_opt = + shr_ptr_bw_engine->build_evolve_functor(); + + amrex::Real* AMREX_RESTRICT p_tau = + pti.GetAttribs(particle_comps["tau"]).dataPtr(); + + const auto me = PhysConst::m_e; + + amrex::ParallelFor( + pti.numParticles(), + [=] AMREX_GPU_DEVICE (long i) { + evolve_opt( + me*ux[i], me*uy[i], me*uz[i], + Ex[i], Ey[i], Ez[i], + Bx[i], By[i], Bz[i], + dt, p_tau[i]); + } + ); +} +#endif |