diff options
Diffstat (limited to 'Source/Particles/PhysicalParticleContainer.cpp')
| -rw-r--r-- | Source/Particles/PhysicalParticleContainer.cpp | 454 |
1 files changed, 389 insertions, 65 deletions
diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index fed7266e1..ee0711372 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -74,13 +74,13 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp #ifdef WARPX_QED //Add real component if QED is enabled - pp.query("do_qed", do_qed); - if(do_qed) + pp.query("do_qed", m_do_qed); + if(m_do_qed) AddRealComp("tau"); //IF do_qed is enabled, find out if Quantum Synchrotron process is enabled - if(do_qed) - pp.query("do_qed_quantum_sync", do_qed_quantum_sync); + if(m_do_qed) + pp.query("do_qed_quantum_sync", m_do_qed_quantum_sync); //TODO: SHOULD CHECK IF SPECIES IS EITHER ELECTRONS OR POSITRONS!! #endif @@ -91,7 +91,7 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp plot_flag_size += 6; #ifdef WARPX_QED - if(do_qed){ + if(m_do_qed){ // plot_flag will have an entry for the optical depth plot_flag_size++; } @@ -123,7 +123,7 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp } #ifdef WARPX_QED - if(do_qed){ + if(m_do_qed){ //Optical depths is always plotted if QED is on plot_flags[plot_flag_size-1] = 1; } @@ -544,11 +544,11 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) BreitWheelerGetOpticalDepth breit_wheeler_get_opt; if(loc_has_quantum_sync){ quantum_sync_get_opt = - shr_ptr_qs_engine->build_optical_depth_functor(); + m_shr_p_qs_engine->build_optical_depth_functor(); } if(loc_has_breit_wheeler){ breit_wheeler_get_opt = - shr_ptr_bw_engine->build_optical_depth_functor(); + m_shr_p_bw_engine->build_optical_depth_functor(); } #endif @@ -1600,11 +1600,48 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti, const Real q = this->charge; const Real m = this-> mass; +#ifdef WARPX_QED + if(has_quantum_sync()){ + Real* AMREX_RESTRICT p_tau = + pti.GetAttribs(particle_comps["tau"]).dataPtr(); + PushPX_QedQuantumSynchrotron(x, y, z, ux, uy, uz, + Ex, Ey, Ez, Bx, By, Bz, + q, m, p_tau, dt, pti.numParticles()); + } + else { + PushPX_classical(x, y, z, ux, uy, uz, + Ex, Ey, Ez, Bx, By, Bz, + q, m, ion_lev, dt, pti.numParticles()); + } +#else + PushPX_classical(x, y, z, ux, uy, uz, + Ex, Ey, Ez, Bx, By, Bz, + q, m, ion_lev, dt, pti.numParticles()); +#endif +} + +void PhysicalParticleContainer::PushPX_classical( + ParticleReal* const AMREX_RESTRICT x, + ParticleReal* const AMREX_RESTRICT y, + ParticleReal* const AMREX_RESTRICT z, + ParticleReal* const AMREX_RESTRICT ux, + ParticleReal* const AMREX_RESTRICT uy, + ParticleReal* const AMREX_RESTRICT uz, + const ParticleReal* const AMREX_RESTRICT Ex, + const ParticleReal* const AMREX_RESTRICT Ey, + const ParticleReal* const AMREX_RESTRICT Ez, + const ParticleReal* const AMREX_RESTRICT Bx, + const ParticleReal* const AMREX_RESTRICT By, + const ParticleReal* const AMREX_RESTRICT Bz, + Real q, Real m, int* AMREX_RESTRICT ion_lev, + amrex::Real dt, long num_particles +) +{ //Assumes that all consistency checks have been done at initialization if(do_classical_radiation_reaction){ amrex::ParallelFor( - pti.numParticles(), + num_particles, [=] AMREX_GPU_DEVICE (long i) { Real qp = q; if (ion_lev){ qp *= ion_lev[i]; } @@ -1617,7 +1654,7 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti, ); } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Boris){ amrex::ParallelFor( - pti.numParticles(), + num_particles, [=] AMREX_GPU_DEVICE (long i) { Real qp = q; if (ion_lev){ qp *= ion_lev[i]; } @@ -1630,7 +1667,7 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti, ); } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Vay) { amrex::ParallelFor( - pti.numParticles(), + num_particles, [=] AMREX_GPU_DEVICE (long i) { Real qp = q; if (ion_lev){ qp *= ion_lev[i]; } @@ -1643,7 +1680,7 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti, ); } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::HigueraCary) { amrex::ParallelFor( - pti.numParticles(), + num_particles, [=] AMREX_GPU_DEVICE (long i) { Real qp = q; if (ion_lev){ qp *= ion_lev[i]; } @@ -1656,9 +1693,141 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti, ); } else { amrex::Abort("Unknown particle pusher"); - }; + } } +#ifdef WARPX_QED +void PhysicalParticleContainer::PushPX_QedQuantumSynchrotron( + ParticleReal* const AMREX_RESTRICT x, + ParticleReal* const AMREX_RESTRICT y, + ParticleReal* const AMREX_RESTRICT z, + ParticleReal* const AMREX_RESTRICT ux, + ParticleReal* const AMREX_RESTRICT uy, + ParticleReal* const AMREX_RESTRICT uz, + const ParticleReal* const AMREX_RESTRICT Ex, + const ParticleReal* const AMREX_RESTRICT Ey, + const ParticleReal* const AMREX_RESTRICT Ez, + const ParticleReal* const AMREX_RESTRICT Bx, + const ParticleReal* const AMREX_RESTRICT By, + const ParticleReal* const AMREX_RESTRICT Bz, + Real q, amrex::Real m, + ParticleReal* const AMREX_RESTRICT tau, + Real dt, long num_particles +) +{ + QuantumSynchrotronEvolveOpticalDepth evolve_opt = + m_shr_p_qs_engine->build_evolve_functor(); + + const auto chi_min = m_shr_p_qs_engine->get_ref_ctrl().chi_part_min; + + //Assumes that all consistency checks have been done at initialization + if(do_classical_radiation_reaction){ + amrex::ParallelFor(num_particles, + [=] AMREX_GPU_DEVICE (long i) { + const ParticleReal px = m * ux[i]; + const ParticleReal py = m * uy[i]; + const ParticleReal pz = m * uz[i]; + + const ParticleReal chi = QedUtils::chi_lepton(px, py, pz, + Ex[i], Ey[i], Ez[i], + Bx[i], By[i], Bz[i]); + + if(chi < chi_min){ + UpdateMomentumBorisWithRadiationReaction( + ux[i], uy[i], uz[i], + Ex[i], Ey[i], Ez[i], + Bx[i], By[i], Bz[i], + q, m, dt); + } + else + { + bool has_event_happened = evolve_opt( + px, py, pz, + Ex[i], Ey[i], Ez[i], + Bx[i], By[i], Bz[i], + dt, tau[i]); + + 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 if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Boris){ + amrex::ParallelFor(num_particles, + [=] 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, tau[i]); + + 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 if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Vay) { + amrex::ParallelFor(num_particles, + [=] 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, tau[i]); + + UpdateMomentumVay( 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 (WarpX::particle_pusher_algo == ParticlePusherAlgo::HigueraCary) { + amrex::ParallelFor( + num_particles, + [=] 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, tau[i]); + + UpdateMomentumHigueraCary( 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::Abort("Unknown particle pusher"); + } +} +#endif + void PhysicalParticleContainer::PushP (int lev, Real dt, const MultiFab& Ex, const MultiFab& Ey, const MultiFab& Ez, @@ -1742,59 +1911,214 @@ PhysicalParticleContainer::PushP (int lev, Real dt, ion_lev = pti.GetiAttribs(particle_icomps["ionization_level"]).dataPtr(); } +#ifdef WARPX_QED + if(has_quantum_sync()){ + Real* AMREX_RESTRICT p_tau = + pti.GetAttribs(particle_comps["tau"]).dataPtr(); + PushP_QedQuantumSynchrotron( + ux, uy, uz, Expp, Eypp, Ezpp, Bxpp, Bypp, Bzpp, + q, m, p_tau, dt, pti.numParticles()); + } + else{ + PushP_classical( + ux, uy, uz, Expp, Eypp, Ezpp, Bxpp, Bypp, Bzpp, + q, m, ion_lev, dt, pti.numParticles()); + } +#else + PushP_classical( + ux, uy, uz, Expp, Eypp, Ezpp, Bxpp, Bypp, Bzpp, + q, m, ion_lev, dt, pti.numParticles()); +#endif + - //Assumes that all consistency checks have been done at initialization - if(do_classical_radiation_reaction){ - amrex::ParallelFor( - pti.numParticles(), - [=] AMREX_GPU_DEVICE (long i) { - Real qp = q; - if (ion_lev){ qp *= ion_lev[i]; } - UpdateMomentumBorisWithRadiationReaction( - ux[i], uy[i], uz[i], - Expp[i], Eypp[i], Ezpp[i], - Bxpp[i], Bypp[i], Bzpp[i], - qp, m, dt); - } - ); - } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Boris){ - 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], - Expp[i], Eypp[i], Ezpp[i], - Bxpp[i], Bypp[i], Bzpp[i], - qp, m, dt); - } - ); - } 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]; } - UpdateMomentumVay( - ux[i], uy[i], uz[i], - Expp[i], Eypp[i], Ezpp[i], - Bxpp[i], Bypp[i], Bzpp[i], - qp, m, dt); - } - ); - } 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); - } - ); - } else { - amrex::Abort("Unknown particle pusher"); - }; } } } +void PhysicalParticleContainer::PushP_classical( + ParticleReal* const AMREX_RESTRICT ux, + ParticleReal* const AMREX_RESTRICT uy, + ParticleReal* const AMREX_RESTRICT uz, + const ParticleReal* const AMREX_RESTRICT Expp, + const ParticleReal* const AMREX_RESTRICT Eypp, + const ParticleReal* const AMREX_RESTRICT Ezpp, + const ParticleReal* const AMREX_RESTRICT Bxpp, + const ParticleReal* const AMREX_RESTRICT Bypp, + const ParticleReal* const AMREX_RESTRICT Bzpp, + Real q, Real m, int* AMREX_RESTRICT ion_lev, + Real dt, long num_particles +) +{ + //Assumes that all consistency checks have been done at initialization + if(do_classical_radiation_reaction){ + amrex::ParallelFor(num_particles, + [=] AMREX_GPU_DEVICE (long i){ + Real qp = q; + if (ion_lev){ qp *= ion_lev[i]; } + UpdateMomentumBorisWithRadiationReaction( + ux[i], uy[i], uz[i], + Expp[i], Eypp[i], Ezpp[i], + Bxpp[i], Bypp[i], Bzpp[i], + qp, m, dt); + } + ); + } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Boris){ + amrex::ParallelFor(num_particles, + [=] AMREX_GPU_DEVICE (long i) { + Real qp = q; + if (ion_lev){ qp *= ion_lev[i]; } + UpdateMomentumBoris( + ux[i], uy[i], uz[i], + Expp[i], Eypp[i], Ezpp[i], + Bxpp[i], Bypp[i], Bzpp[i], + qp, m, dt); + } + ); + } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Vay){ + amrex::ParallelFor(num_particles, + [=] AMREX_GPU_DEVICE (long i) { + Real qp = q; + if (ion_lev){ qp *= ion_lev[i]; } + UpdateMomentumVay( + ux[i], uy[i], uz[i], + Expp[i], Eypp[i], Ezpp[i], + Bxpp[i], Bypp[i], Bzpp[i], + qp, m, dt); + } + ); + } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::HigueraCary){ + amrex::ParallelFor(num_particles, + [=] 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); + } + ); + } else { + amrex::Abort("Unknown particle pusher"); + } +} + +#ifdef WARPX_QED + +void PhysicalParticleContainer::PushP_QedQuantumSynchrotron( + ParticleReal* const AMREX_RESTRICT ux, + ParticleReal* const AMREX_RESTRICT uy, + ParticleReal* const AMREX_RESTRICT uz, + const ParticleReal* const AMREX_RESTRICT Expp, + const ParticleReal* const AMREX_RESTRICT Eypp, + const ParticleReal* const AMREX_RESTRICT Ezpp, + const ParticleReal* const AMREX_RESTRICT Bxpp, + const ParticleReal* const AMREX_RESTRICT Bypp, + const ParticleReal* const AMREX_RESTRICT Bzpp, + Real q, Real m, + ParticleReal* const AMREX_RESTRICT tau, + Real dt, long num_particles +) +{ + QuantumSynchrotronEvolveOpticalDepth evolve_opt = + m_shr_p_qs_engine->build_evolve_functor(); + + const auto chi_min = m_shr_p_qs_engine->get_ref_ctrl().chi_part_min; + + //Assumes that all consistency checks have been done at initialization + if(do_classical_radiation_reaction){ + amrex::ParallelFor(num_particles, + [=] AMREX_GPU_DEVICE (long i) { + const ParticleReal px = m * ux[i]; + const ParticleReal py = m * uy[i]; + const ParticleReal pz = m * uz[i]; + + const ParticleReal chi = QedUtils::chi_lepton(px, py, pz, + Expp[i], Eypp[i], Ezpp[i], + Bxpp[i], Bypp[i], Bzpp[i]); + + if(chi < chi_min){ + UpdateMomentumBorisWithRadiationReaction( + ux[i], uy[i], uz[i], + Expp[i], Eypp[i], Ezpp[i], + Bxpp[i], Bypp[i], Bzpp[i], + q, m, dt); + } + else + { + bool has_event_happened = evolve_opt( + px, py, pz, + Expp[i], Eypp[i], Ezpp[i], + Bxpp[i], Bypp[i], Bzpp[i], + dt, tau[i]); + + UpdateMomentumBoris( + ux[i], uy[i], uz[i], + Expp[i], Eypp[i], Ezpp[i], + Bxpp[i], Bypp[i], Bzpp[i], + q, m, dt); + } + } + ); + } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Boris){ + amrex::ParallelFor(num_particles, + [=] 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, + Expp[i], Eypp[i], Ezpp[i], + Bxpp[i], Bypp[i], Bzpp[i], + dt, tau[i]); + + UpdateMomentumBoris( + ux[i], uy[i], uz[i], + Expp[i], Eypp[i], Ezpp[i], + Bxpp[i], Bypp[i], Bzpp[i], + q, m, dt); + } + ); + } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Vay){ + amrex::ParallelFor(num_particles, + [=] 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, + Expp[i], Eypp[i], Ezpp[i], + Bxpp[i], Bypp[i], Bzpp[i], + dt, tau[i]); + + UpdateMomentumVay( + ux[i], uy[i], uz[i], + Expp[i], Eypp[i], Ezpp[i], + Bxpp[i], Bypp[i], Bzpp[i], + q, m, dt); + } + ); + } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::HigueraCary){ + amrex::ParallelFor(num_particles, + [=] 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, + Expp[i], Eypp[i], Ezpp[i], + Bxpp[i], Bypp[i], Bzpp[i], + dt, tau[i]); + + UpdateMomentumHigueraCary( ux[i], uy[i], uz[i], + Expp[i], Eypp[i], Ezpp[i], Bxpp[i], Bypp[i], Bzpp[i], q, m, dt); + } + ); + } else { + amrex::Abort("Unknown particle pusher"); + } +} + +#endif + void PhysicalParticleContainer::copy_attribs(WarpXParIter& pti,const ParticleReal* xp, const ParticleReal* yp, const ParticleReal* zp) { @@ -2238,25 +2562,25 @@ PhysicalParticleContainer::AmIALepton(){ bool PhysicalParticleContainer::has_quantum_sync() { - return do_qed_quantum_sync; + return m_do_qed_quantum_sync; } bool PhysicalParticleContainer::has_breit_wheeler() { - return do_qed_breit_wheeler; + return m_do_qed_breit_wheeler; } void PhysicalParticleContainer:: set_breit_wheeler_engine_ptr(std::shared_ptr<BreitWheelerEngine> ptr) { - shr_ptr_bw_engine = ptr; + m_shr_p_bw_engine = ptr; } void PhysicalParticleContainer:: set_quantum_sync_engine_ptr(std::shared_ptr<QuantumSynchrotronEngine> ptr) { - shr_ptr_qs_engine = ptr; + m_shr_p_qs_engine = ptr; } #endif |