diff options
Diffstat (limited to 'Source/Particles')
-rw-r--r-- | Source/Particles/MultiParticleContainer.H | 1 | ||||
-rw-r--r-- | Source/Particles/PhotonParticleContainer.H | 14 | ||||
-rw-r--r-- | Source/Particles/PhotonParticleContainer.cpp | 51 | ||||
-rw-r--r-- | Source/Particles/PhysicalParticleContainer.H | 38 | ||||
-rw-r--r-- | Source/Particles/PhysicalParticleContainer.cpp | 119 | ||||
-rw-r--r-- | Source/Particles/WarpXParticleContainer.H | 3 |
6 files changed, 208 insertions, 18 deletions
diff --git a/Source/Particles/MultiParticleContainer.H b/Source/Particles/MultiParticleContainer.H index f9a0e51d7..06898a5c1 100644 --- a/Source/Particles/MultiParticleContainer.H +++ b/Source/Particles/MultiParticleContainer.H @@ -11,6 +11,7 @@ #include <AMReX_Particles.H> #ifdef WARPX_QED + #include <QedChiFunctions.H> #include <BreitWheelerEngineWrapper.H> #include <QuantumSyncEngineWrapper.H> #endif diff --git a/Source/Particles/PhotonParticleContainer.H b/Source/Particles/PhotonParticleContainer.H index 9132cf4a9..6a474a92a 100644 --- a/Source/Particles/PhotonParticleContainer.H +++ b/Source/Particles/PhotonParticleContainer.H @@ -46,7 +46,7 @@ public: amrex::Cuda::ManagedDeviceVector<amrex::ParticleReal>& zp, amrex::Real dt, DtType a_dt_type=DtType::Full) override; - // Don't push momenta for photons + // Do nothing virtual void PushP (int lev, amrex::Real dt, const amrex::MultiFab& Ex, @@ -79,6 +79,18 @@ public: return false; }; +#ifdef WARPX_QED + /** + * This function evolves the optical depth of the photons if QED effects + * are enabled. + * @param[in,out] pti particle iterator (optical depth will be modified) + * @param[in] dt temporal step + */ + virtual void EvolveOpticalDepth(WarpXParIter& pti, + amrex::Real dt) override; + +#endif + }; #endif // #ifndef WARPX_PhotonParticleContainer_H_ diff --git a/Source/Particles/PhotonParticleContainer.cpp b/Source/Particles/PhotonParticleContainer.cpp index 7e52b52e1..7d63fff11 100644 --- a/Source/Particles/PhotonParticleContainer.cpp +++ b/Source/Particles/PhotonParticleContainer.cpp @@ -78,8 +78,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) { @@ -88,7 +86,6 @@ PhotonParticleContainer::PushPX(WarpXParIter& pti, ux[i], uy[i], uz[i], dt ); } ); - } void @@ -116,3 +113,51 @@ PhotonParticleContainer::Evolve (int lev, t, dt); } + +#ifdef WARPX_QED + +void +PhotonParticleContainer::EvolveOpticalDepth( + WarpXParIter& pti,amrex::Real dt) +{ + //m_shr_p_bw_engine->are_lookup_tables_initialized() is necessary here if we want + //to perform just initialization tests of the optical depth without actually + //enabling QED effects (this requires lookup tables). + if(!has_breit_wheeler() || !m_shr_p_bw_engine->are_lookup_tables_initialized()) + return; + + 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 = + m_shr_p_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) { + const ParticleReal px = me * ux[i]; + const ParticleReal py = me * uy[i]; + const ParticleReal pz = me * uz[i]; + + bool has_event_happened = evolve_opt( + px, py, pz, + Ex[i], Ey[i], Ez[i], + Bx[i], By[i], Bz[i], + dt, p_tau[i]); + } + ); +} +#endif diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index 17a504719..4701c559d 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -187,13 +187,47 @@ public: amrex::FArrayBox const * & byfab, amrex::FArrayBox const * & bzfab); #ifdef WARPX_QED + //Functions decleared in WarpXParticleContainer.H + //containers for which QED processes could be relevant + //are expected to override these functions + + /** + * Tells if this PhysicalParticleContainer has Quantum + * Synchrotron process enabled + * @return true if process is enabled + */ bool has_quantum_sync() override; + + /** + * Tells if this PhysicalParticleContainer has Breit + * Wheeler process enabled + * @return true if process is enabled + */ bool has_breit_wheeler() override; + /** + * Acquires a shared smart pointer to a BreitWheelerEngine + * @param[in] ptr the pointer + */ void set_breit_wheeler_engine_ptr - (std::shared_ptr<BreitWheelerEngine>) override; + (std::shared_ptr<BreitWheelerEngine> ptr) override; + + /** + * Acquires a shared smart pointer to a QuantumSynchrotronEngine + * @param[in] ptr the pointer + */ void set_quantum_sync_engine_ptr - (std::shared_ptr<QuantumSynchrotronEngine>) override; + (std::shared_ptr<QuantumSynchrotronEngine> ptr) override; + //__________ + + /** + * This function evolves the optical depth of the particles if QED effects + * are enabled. + * @param[in,out] pti particle iterator (optical depth will be modified) + * @param[in] dt temporal step + */ + virtual void EvolveOpticalDepth(WarpXParIter& pti, + amrex::Real dt); #endif protected: diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 51690d659..54cc1a880 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1062,6 +1062,7 @@ PhysicalParticleContainer::Evolve (int lev, BL_PROFILE("PPC::Evolve()"); BL_PROFILE_VAR_NS("PPC::Evolve::Copy", blp_copy); BL_PROFILE_VAR_NS("PPC::FieldGather", blp_fg); + BL_PROFILE_VAR_NS("PPC::EvolveOpticalDepth", blp_ppc_qed_ev); BL_PROFILE_VAR_NS("PPC::ParticlePush", blp_ppc_pp); const std::array<Real,3>& dx = WarpX::CellSize(lev); @@ -1253,6 +1254,15 @@ PhysicalParticleContainer::Evolve (int lev, BL_PROFILE_VAR_STOP(blp_fg); +#ifdef WARPX_QED + // + //Evolve Optical Depth + // + BL_PROFILE_VAR_START(blp_ppc_qed_ev); + EvolveOpticalDepth(pti, dt); + BL_PROFILE_VAR_STOP(blp_ppc_qed_ev); +#endif + // // Particle Push // @@ -1600,8 +1610,45 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti, const Real q = this->charge; const Real m = this-> mass; +#ifdef WARPX_QED - //Assumes that all consistency checks have been done at initialization + auto t_chi_max = m_shr_p_qs_engine->get_ref_ctrl().chi_part_min; + + if(do_classical_radiation_reaction){ + if(m_do_qed_quantum_sync){ + amrex::ParallelFor( + pti.numParticles(), + [=] AMREX_GPU_DEVICE (long i) { + auto chi = QedUtils::chi_lepton(m*ux[i], m*uy[i], m*uz[i], + Ex[i], Ey[i], Ez[i], + Bx[i], By[i], Bz[i]); + if(chi < t_chi_max){ + UpdateMomentumBorisWithRadiationReaction( ux[i], uy[i], uz[i], + Ex[i], Ey[i], Ez[i], Bx[i], + By[i], Bz[i], q, m, dt); + } + else{ + UpdateMomentumBoris( ux[i], uy[i], uz[i], + Ex[i], Ey[i], Ez[i], Bx[i], + By[i], Bz[i], q, m, dt); + } + UpdatePosition( x[i], y[i], z[i], + ux[i], uy[i], uz[i], dt ); + } + ); + }else{ + amrex::ParallelFor( + pti.numParticles(), + [=] AMREX_GPU_DEVICE (long i) { + UpdateMomentumBorisWithRadiationReaction( ux[i], uy[i], uz[i], + Ex[i], Ey[i], Ez[i], Bx[i], + By[i], Bz[i], q, m, dt); + UpdatePosition( x[i], y[i], z[i], + ux[i], uy[i], uz[i], dt ); + } + ); + } +#else if(do_classical_radiation_reaction){ amrex::ParallelFor( pti.numParticles(), @@ -1615,6 +1662,7 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti, ux[i], uy[i], uz[i], dt ); } ); +#endif } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Boris){ amrex::ParallelFor( pti.numParticles(), @@ -1659,6 +1707,52 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti, }; } +#ifdef WARPX_QED +void PhysicalParticleContainer::EvolveOpticalDepth( + WarpXParIter& pti, amrex::Real dt) +{ + //m_shr_p_qs_engine->are_lookup_tables_initialized() is necessary here if we want + //to perform just initialization tests of the optical depth without actually + //enabling QED effects (this requires lookup tables). + if(!has_quantum_sync() || !m_shr_p_qs_engine->are_lookup_tables_initialized()) + return; + + QuantumSynchrotronEvolveOpticalDepth evolve_opt = + m_shr_p_qs_engine->build_evolve_functor(); + + auto& attribs = pti.GetAttribs(); + const ParticleReal* const AMREX_RESTRICT ux = attribs[PIdx::ux].dataPtr(); + const ParticleReal* const AMREX_RESTRICT uy = attribs[PIdx::uy].dataPtr(); + const 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(); + + ParticleReal* const AMREX_RESTRICT p_tau = + pti.GetAttribs(particle_comps["tau"]).dataPtr(); + + const ParticleReal m = this->mass; + + amrex::ParallelFor(pti.numParticles(), + [=] AMREX_GPU_DEVICE (long i) { + const ParticleReal px = m * ux[i]; + const ParticleReal py = m * uy[i]; + const ParticleReal pz = m * uz[i]; + + bool has_event_happened = evolve_opt( + px, py, pz, + Ex[i], Ey[i], Ez[i], + Bx[i], By[i], Bz[i], + dt, p_tau[i]); + } + ); + +} +#endif + void PhysicalParticleContainer::PushP (int lev, Real dt, const MultiFab& Ex, const MultiFab& Ey, const MultiFab& Ez, @@ -1742,12 +1836,10 @@ PhysicalParticleContainer::PushP (int lev, Real dt, ion_lev = pti.GetiAttribs(particle_icomps["ionization_level"]).dataPtr(); } - //Assumes that all consistency checks have been done at initialization if(do_classical_radiation_reaction){ - amrex::ParallelFor( - pti.numParticles(), - [=] AMREX_GPU_DEVICE (long i) { + amrex::ParallelFor(pti.numParticles(), + [=] AMREX_GPU_DEVICE (long i){ Real qp = q; if (ion_lev){ qp *= ion_lev[i]; } UpdateMomentumBorisWithRadiationReaction( @@ -1758,7 +1850,7 @@ PhysicalParticleContainer::PushP (int lev, Real dt, } ); } else 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]; } @@ -1769,8 +1861,8 @@ PhysicalParticleContainer::PushP (int lev, Real dt, qp, m, dt); } ); - } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Vay) { - amrex::ParallelFor( pti.numParticles(), + } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Vay){ + amrex::ParallelFor(pti.numParticles(), [=] AMREX_GPU_DEVICE (long i) { Real qp = q; if (ion_lev){ qp *= ion_lev[i]; } @@ -1781,16 +1873,19 @@ PhysicalParticleContainer::PushP (int lev, Real dt, qp, m, dt); } ); - } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::HigueraCary) { - amrex::ParallelFor( pti.numParticles(), + } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::HigueraCary){ + amrex::ParallelFor(pti.numParticles(), [=] AMREX_GPU_DEVICE (long i) { UpdateMomentumHigueraCary( ux[i], uy[i], uz[i], - Expp[i], Eypp[i], Ezpp[i], Bxpp[i], Bypp[i], Bzpp[i], q, m, dt); + Expp[i], Eypp[i], Ezpp[i], + Bxpp[i], Bypp[i], Bzpp[i], + q, m, dt); } ); } else { amrex::Abort("Unknown particle pusher"); - }; + } + } } } diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index dbd913c5b..7eaf7c277 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -322,9 +322,12 @@ protected: #ifdef WARPX_QED bool m_do_qed = false; + //Species for which QED effects are relevant should override these methods virtual bool has_quantum_sync(){return false;}; virtual bool has_breit_wheeler(){return false;}; + //Species can receive a shared pointer to a QED engine (species for + //which this is relevant should override these functions) virtual void set_breit_wheeler_engine_ptr(std::shared_ptr<BreitWheelerEngine>){}; virtual void |