add RR force implicitly

RevathiJambunathan · Dec 22, 2023 · b862443 · b862443
1 parent 2874189
commit b862443
Showing 1 changed file with 17 additions and 7 deletions.
diff --git a/Source/Particles/Pusher/UpdateMomentumBorisWithRadiationReaction.H b/Source/Particles/Pusher/UpdateMomentumBorisWithRadiationReaction.H
@@ -105,20 +105,30 @@ void UpdateMomentumBorisWithRadiationReaction(
     const amrex::ParticleReal RRcoeff = (2.0_prt/3.0_prt)*PhysConst::r_e*q_over_mc*q_over_mc;
 #endif
 
+
+#ifdef PULSAR
+    //Update momentum using the RR force
+    if (r_p < scaledRR_rmax) {
+        amrex::Real lambda_n = gamma_n*(fl_q2-bdotE2) / (PhysConst::c);
+        if (lambda_n > 0) {
+            amrex::Real RRfactor = 1._prt/(1._prt + RRcoeff*lambda_n*dt);
+            ux = RRfactor * ux;
+            uy = RRfactor * uy;
+            uz = RRfactor * uz;
+        }    
+    }
+#else
     //Compute the components of the RR force
     const amrex::ParticleReal frx =
         RRcoeff*(PhysConst::c*(fly_q*Bz - flz_q*By) + bdotE*Ex - coeff*bx_n);
     const amrex::ParticleReal fry =
         RRcoeff*(PhysConst::c*(flz_q*Bx - flx_q*Bz) + bdotE*Ey - coeff*by_n);
     const amrex::ParticleReal frz =
         RRcoeff*(PhysConst::c*(flx_q*By - fly_q*Bx) + bdotE*Ez - coeff*bz_n);
-
-    //Update momentum using the RR force
-    if (r_p < scaledRR_rmax) {
-        ux += frx*dt;
-        uy += fry*dt;
-        uz += frz*dt;
-    }
+    ux += frx*dt;
+    uy += fry*dt;
+    uz += frz*dt;
+#endif
 }
 
 #endif // WARPX_PARTICLES_PUSHER_UPDATEMOMENTUM_BORIS_WITHRR_H_