diff options
Diffstat (limited to 'Source')
5 files changed, 201 insertions, 11 deletions
diff --git a/Source/Particles/PhotonParticleContainer.H b/Source/Particles/PhotonParticleContainer.H index 407cf26f3..9132cf4a9 100644 --- a/Source/Particles/PhotonParticleContainer.H +++ b/Source/Particles/PhotonParticleContainer.H @@ -73,6 +73,12 @@ public: int lev, int depos_lev, amrex::Real dt) override {}; + //Photons are not leptons + virtual bool AmIALepton () override + { + return false; + }; + }; #endif // #ifndef WARPX_PhotonParticleContainer_H_ diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index 08493edcd..ed3a698e2 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -210,6 +210,14 @@ protected: // Inject particles during the whole simulation void ContinuousInjection (const amrex::RealBox& injection_box) override; + //This function return true if the PhysicalParticleContainer contains electrons + //or positrons, false otherwise + virtual bool AmIALepton (); + + //When true PhysicalParticleContainer tries to use a pusher including + //radiation reaction + bool do_classical_radiation_reaction = false; + #ifdef WARPX_QED // A flag to enable quantum_synchrotron process for leptons bool do_qed_quantum_sync = false; diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 66331db26..d12a4dbff 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -17,6 +17,7 @@ #include <UpdatePosition.H> #include <UpdateMomentumBoris.H> #include <UpdateMomentumVay.H> +#include <UpdateMomentumBorisWithRadiationReaction.H> #include <UpdateMomentumHigueraCary.H> using namespace amrex; @@ -45,6 +46,23 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp pp.query("do_boosted_frame_diags", do_boosted_frame_diags); pp.query("do_field_ionization", do_field_ionization); + + //check if Radiation Reaction is enabled and do consistency checks + pp.query("do_classical_radiation_reaction", do_classical_radiation_reaction); + //if the species is not a lepton, do_classical_radiation_reaction + //should be false + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + !(do_classical_radiation_reaction && !AmIALepton()), + "Can't enable classical radiation reaction for non lepton species. " ); + + //Only Boris pusher is compatible with radiation reaction + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + !(do_classical_radiation_reaction && + WarpX::particle_pusher_algo != ParticlePusherAlgo::Boris), + "Radiation reaction can be enabled only if Boris pusher is used"); + //_____________________________ + + #ifdef AMREX_USE_GPU AMREX_ALWAYS_ASSERT_WITH_MESSAGE( do_field_ionization == 0, @@ -1580,7 +1598,23 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti, // Loop over the particles and update their momentum const Real q = this->charge; const Real m = this-> mass; - if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Boris){ + + + //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], + Ex[i], Ey[i], Ez[i], Bx[i], + By[i], Bz[i], qp, 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( pti.numParticles(), [=] AMREX_GPU_DEVICE (long i) { @@ -1701,18 +1735,49 @@ PhysicalParticleContainer::PushP (int lev, Real dt, // Loop over the particles and update their momentum const Real q = this->charge; const Real m = this-> mass; - if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Boris){ + + int* AMREX_RESTRICT ion_lev = nullptr; + if (do_field_ionization){ + 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) { + 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) { - UpdateMomentumBoris( ux[i], uy[i], uz[i], - Expp[i], Eypp[i], Ezpp[i], Bxpp[i], Bypp[i], Bzpp[i], q, m, dt); + 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) { - UpdateMomentumVay( ux[i], uy[i], uz[i], - Expp[i], Eypp[i], Ezpp[i], Bxpp[i], Bypp[i], Bzpp[i], q, m, dt); + 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) { @@ -2163,6 +2228,13 @@ PhysicalParticleContainer::buildIonizationMask (const amrex::MFIter& mfi, const ); } +//This function return true if the PhysicalParticleContainer contains electrons +//or positrons, false otherwise +bool +PhysicalParticleContainer::AmIALepton(){ + return (this-> mass == PhysConst::m_e); +} + #ifdef WARPX_QED bool PhysicalParticleContainer::has_quantum_sync() diff --git a/Source/Particles/Pusher/UpdateMomentumBorisWithRadiationReaction.H b/Source/Particles/Pusher/UpdateMomentumBorisWithRadiationReaction.H new file mode 100644 index 000000000..0abedd258 --- /dev/null +++ b/Source/Particles/Pusher/UpdateMomentumBorisWithRadiationReaction.H @@ -0,0 +1,84 @@ + +#ifndef WARPX_PARTICLES_PUSHER_UPDATEMOMENTUM_BORIS_WITHRR_H_ +#define WARPX_PARTICLES_PUSHER_UPDATEMOMENTUM_BORIS_WITHRR_H_ + +#include <AMReX_REAL.H> +#include <UpdateMomentumBoris.H> + +/** + * Push the particle's positions over one timestep, + * given the value of its momenta `ux`, `uy`, `uz`. + * Includes Radiation Reaction according to + * https://doi.org/10.1088/1367-2630/12/12/123005 + */ +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void UpdateMomentumBorisWithRadiationReaction( + amrex::ParticleReal& ux, amrex::ParticleReal& uy, amrex::ParticleReal& uz, + const amrex::ParticleReal Ex, const amrex::ParticleReal Ey, const amrex::ParticleReal Ez, + const amrex::ParticleReal Bx, const amrex::ParticleReal By, const amrex::ParticleReal Bz, + const amrex::Real q, const amrex::Real m, const amrex::Real dt ) +{ + //RR algorithm needs to store old value of the normalized momentum + const amrex::Real ux_old = ux; + const amrex::Real uy_old = uy; + const amrex::Real uz_old = uz; + + //Useful constants + const amrex::Real inv_dt = 1.0/dt; + constexpr amrex::Real inv_c2 = 1./(PhysConst::c*PhysConst::c); + + //Call to regular Boris pusher + UpdateMomentumBoris( + ux, uy, uz, + Ex, Ey, Ez, + Bx, By, Bz, + q, m, dt ); + + //Estimation of the normalized momentum at intermediate (integer) time + const amrex::Real ux_n = (ux+ux_old)*0.5; + const amrex::Real uy_n = (uy+uy_old)*0.5; + const amrex::Real uz_n = (uz+uz_old)*0.5; + + // Compute Lorentz factor (and inverse) at intermediate (integer) time + const amrex::Real gamma_n = std::sqrt( 1. + + (ux_n*ux_n + uy_n*uy_n + uz_n*uz_n)*inv_c2); + const amrex::Real inv_gamma_n = 1.0/gamma_n; + + //Estimation of the velocity at intermediate (integer) time + const amrex::Real vx_n = ux_n*inv_gamma_n; + const amrex::Real vy_n = uy_n*inv_gamma_n; + const amrex::Real vz_n = uz_n*inv_gamma_n; + const amrex::Real bx_n = vx_n/PhysConst::c; + const amrex::Real by_n = vy_n/PhysConst::c; + const amrex::Real bz_n = vz_n/PhysConst::c; + + //Lorentz force over charge + const amrex::Real flx_q = (Ex + vy_n*Bz - vz_n*By); + const amrex::Real fly_q = (Ey + vz_n*Bx - vx_n*Bz); + const amrex::Real flz_q = (Ez + vx_n*By - vy_n*Bx); + const amrex::Real fl_q2 = flx_q*flx_q + fly_q*fly_q + flz_q*flz_q; + + //Calculation of auxiliary quantities + const amrex::Real bdotE = (bx_n*Ex + by_n*Ey + bz_n*Ez); + const amrex::Real bdotE2 = bdotE*bdotE; + const amrex::Real coeff = gamma_n*gamma_n*(fl_q2-bdotE2); + + //Radiation reaction constant + const amrex::Real RRcoeff = 2.0*PhysConst::r_e*q*q/ + (3.0*m*m*PhysConst::c*PhysConst::c); + + //Compute the components of the RR force + const amrex::Real frx = + RRcoeff*(PhysConst::c*(fly_q*Bz - flz_q*By) + bdotE*Ex - coeff*bx_n); + const amrex::Real fry = + RRcoeff*(PhysConst::c*(flz_q*Bx - flx_q*Bz) + bdotE*Ey - coeff*by_n); + const amrex::Real frz = + RRcoeff*(PhysConst::c*(flx_q*By - fly_q*Bx) + bdotE*Ez - coeff*bz_n); + + //Update momentum using the RR force + ux += frx*dt; + uy += fry*dt; + uz += frz*dt; +} + +#endif // WARPX_PARTICLES_PUSHER_UPDATEMOMENTUM_BORIS_WITHRR_H_ diff --git a/Source/Particles/RigidInjectedParticleContainer.cpp b/Source/Particles/RigidInjectedParticleContainer.cpp index 2ef833151..535ffec6f 100644 --- a/Source/Particles/RigidInjectedParticleContainer.cpp +++ b/Source/Particles/RigidInjectedParticleContainer.cpp @@ -13,6 +13,7 @@ #include <WarpXAlgorithmSelection.H> #include <UpdateMomentumBoris.H> #include <UpdateMomentumVay.H> +#include <UpdateMomentumBorisWithRadiationReaction.H> #include <UpdateMomentumHigueraCary.H> using namespace amrex; @@ -430,18 +431,37 @@ RigidInjectedParticleContainer::PushP (int lev, Real dt, // Loop over the particles and update their momentum const Real q = this->charge; const Real m = this->mass; - if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Boris){ + + //Assumes that all consistency checks have been done at initialization + if(do_classical_radiation_reaction){ + amrex::ParallelFor( + pti.numParticles(), + [=] AMREX_GPU_DEVICE (long i) { + UpdateMomentumBorisWithRadiationReaction( + uxpp[i], uypp[i], uzpp[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( pti.numParticles(), [=] AMREX_GPU_DEVICE (long i) { - UpdateMomentumBoris( uxpp[i], uypp[i], uzpp[i], - Expp[i], Eypp[i], Ezpp[i], Bxpp[i], Bypp[i], Bzpp[i], q, m, dt); + UpdateMomentumBoris( + uxpp[i], uypp[i], uzpp[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( pti.numParticles(), [=] AMREX_GPU_DEVICE (long i) { - UpdateMomentumVay( uxpp[i], uypp[i], uzpp[i], - Expp[i], Eypp[i], Ezpp[i], Bxpp[i], Bypp[i], Bzpp[i], q, m, dt); + UpdateMomentumVay( + uxpp[i], uypp[i], uzpp[i], + Expp[i], Eypp[i], Ezpp[i], + Bxpp[i], Bypp[i], Bzpp[i], + q, m, dt); } ); } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::HigueraCary) { |