diff options
Diffstat (limited to 'Source/Particles/PhysicalParticleContainer.cpp')
| -rw-r--r-- | Source/Particles/PhysicalParticleContainer.cpp | 430 |
1 files changed, 377 insertions, 53 deletions
diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 51690d659..a02ac8593 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -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 = + shr_ptr_qs_engine->build_evolve_functor(); + + const auto chi_min = shr_ptr_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 = + shr_ptr_qs_engine->build_evolve_functor(); + + const auto chi_min = shr_ptr_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) { |