aboutsummaryrefslogtreecommitdiff
path: root/Source/Particles
diff options
context:
space:
mode:
Diffstat (limited to 'Source/Particles')
-rw-r--r--Source/Particles/MultiParticleContainer.H1
-rw-r--r--Source/Particles/PhotonParticleContainer.H14
-rw-r--r--Source/Particles/PhotonParticleContainer.cpp51
-rw-r--r--Source/Particles/PhysicalParticleContainer.H38
-rw-r--r--Source/Particles/PhysicalParticleContainer.cpp119
-rw-r--r--Source/Particles/WarpXParticleContainer.H3
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