1 files changed, 45 insertions, 4 deletions

diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp
index 02dee1967..89953120d 100644
--- a/Source/Particles/PhysicalParticleContainer.cpp
+++ b/Source/Particles/PhysicalParticleContainer.cpp
@@ -15,6 +15,7 @@
 #include <UpdatePosition.H>
 #include <UpdateMomentumBoris.H>
 #include <UpdateMomentumVay.H>
+#include <UpdateMomentumBorisWithRadiationReaction.H>
 
 using namespace amrex;
 
@@ -1486,7 +1487,14 @@ 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){
+
+
+    // Boris + RR is enabled only for leptons
+    auto algo = WarpX::particle_pusher_algo;
+    if (!AmIALepton() && algo == ParticlePusherAlgo::BorisRR)
+        algo = ParticlePusherAlgo::Boris;
+
+    if (algo == ParticlePusherAlgo::Boris){
         amrex::ParallelFor(
             pti.numParticles(),
             [=] AMREX_GPU_DEVICE (long i) {
@@ -1499,7 +1507,7 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti,
                       ux[i], uy[i], uz[i], dt );
             }
         );
-    } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Vay) {
+    } else if (algo == ParticlePusherAlgo::Vay) {
         amrex::ParallelFor(
             pti.numParticles(),
             [=] AMREX_GPU_DEVICE (long i) {
@@ -1512,6 +1520,19 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti,
                                 ux[i], uy[i], uz[i], dt );
             }
         );
+    } else if (algo == ParticlePusherAlgo::BorisRR) {
+        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 {
       amrex::Abort("Unknown particle pusher");
     };
@@ -1593,20 +1614,33 @@ 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){
+
+            // Boris + RR is enabled only for leptons
+            auto algo = WarpX::particle_pusher_algo;
+            if (!AmIALepton() && algo == ParticlePusherAlgo::BorisRR)
+                algo = ParticlePusherAlgo::Boris;
+
+            if (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);
                     }
                 );
-            } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Vay) {
+            } else if (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);
                     }
                 );
+            } else if (algo == ParticlePusherAlgo::BorisRR) {
+                amrex::ParallelFor( pti.numParticles(),
+                    [=] AMREX_GPU_DEVICE (long i) {
+                        UpdateMomentumBorisWithRadiationReaction( 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");
             };
@@ -2047,3 +2081,10 @@ 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);
+}