3 files changed, 103 insertions, 0 deletions
diff --git a/Source/Particles/Pusher/Make.package b/Source/Particles/Pusher/Make.package
index 8c8e77905..95a38fa2d 100644
--- a/Source/Particles/Pusher/Make.package
+++ b/Source/Particles/Pusher/Make.package
@@ -1,4 +1,6 @@
 CEXE_headers += GetAndSetPosition.H
 CEXE_headers += UpdatePosition.H
+CEXE_headers += UpdateMomentumBoris.H
+CEXE_headers += UpdateMomentumVay.H
 INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Particles/Pusher
 VPATH_LOCATIONS   += $(WARPX_HOME)/Source/Particles/Pusher
diff --git a/Source/Particles/Pusher/UpdateMomentumBoris.H b/Source/Particles/Pusher/UpdateMomentumBoris.H
new file mode 100644
index 000000000..71e9a8ed1
--- /dev/null
+++ b/Source/Particles/Pusher/UpdateMomentumBoris.H
@@ -0,0 +1,47 @@
+#ifndef WARPX_PARTICLES_PUSHER_UPDATEMOMENTUM_BORIS_H_
+#define WARPX_PARTICLES_PUSHER_UPDATEMOMENTUM_BORIS_H_
+
+#include <AMReX_REAL.H>
+
+/* \brief Push the particle's positions over one timestep,
+ *    given the value of its momenta `ux`, `uy`, `uz` */
+AMREX_GPU_HOST_DEVICE AMREX_INLINE
+void UpdateMomentumBoris(
+    amrex::Real& ux, amrex::Real& uy, amrex::Real& uz, amrex::Real& gaminv,
+    const amrex::Real Ex, const amrex::Real Ey, const amrex::Real Ez,
+    const amrex::Real Bx, const amrex::Real By, const amrex::Real Bz,
+    const amrex::Real q, const amrex::Real m, const amrex::Real dt )
+{
+    const amrex::Real econst = 0.5*q*dt/m;
+
+    // First half-push for E
+    ux += econst*Ex;
+    uy += econst*Ey;
+    uz += econst*Ez;
+    // Compute temporary gamma factor
+    constexpr amrex::Real inv_c2 = 1./(PhysConst::c*PhysConst::c);
+    const amrex::Real inv_gamma = 1./std::sqrt(1. + (ux*ux + uy*uy + uz*uz)*inv_c2);
+    // Magnetic rotation
+    // - Compute temporary variables
+    const amrex::Real tx = econst*inv_gamma*Bx;
+    const amrex::Real ty = econst*inv_gamma*By;
+    const amrex::Real tz = econst*inv_gamma*Bz;
+    const amrex::Real tsqi = 2./(1. + tx*tx + ty*ty + tz*tz);
+    const amrex::Real sx = tx*tsqi;
+    const amrex::Real sy = ty*tsqi;
+    const amrex::Real sz = tz*tsqi;
+    const amrex::Real ux_p = ux + uy*tz - uz*ty;
+    const amrex::Real uy_p = uy + uz*tx - ux*tz;
+    const amrex::Real uz_p = uz + ux*ty - uy*tx;
+    // - Update momentum
+    ux += uy_p*sz - uz_p*sy;
+    uy += uz_p*sx - ux_p*sz;
+    uz += ux_p*sy - uy_p*sx;
+    // Second half-push for E
+    ux += econst*Ex;
+    uy += econst*Ey;
+    uz += econst*Ez;
+    gaminv = 1./std::sqrt(1. + (ux*ux + uy*uy + uz*uz)*inv_c2);
+}
+
+#endif // WARPX_PARTICLES_PUSHER_UPDATEMOMENTUM_BORIS_H_
diff --git a/Source/Particles/Pusher/UpdateMomentumVay.H b/Source/Particles/Pusher/UpdateMomentumVay.H
new file mode 100644
index 000000000..044297e22
--- /dev/null
+++ b/Source/Particles/Pusher/UpdateMomentumVay.H
@@ -0,0 +1,54 @@
+#ifndef WARPX_PARTICLES_PUSHER_UPDATEMOMENTUM_VAY_H_
+#define WARPX_PARTICLES_PUSHER_UPDATEMOMENTUM_VAY_H_
+
+#include <AMReX_FArrayBox.H>
+#include <WarpXConst.H>
+#include <AMReX_REAL.H>
+
+/* \brief Push the particle's positions over one timestep,
+ *    given the value of its momenta `ux`, `uy`, `uz` */
+AMREX_GPU_HOST_DEVICE AMREX_INLINE
+void UpdateMomentumVay(
+    amrex::Real& ux, amrex::Real& uy, amrex::Real& uz, amrex::Real& gaminv,
+    const amrex::Real Ex, const amrex::Real Ey, const amrex::Real Ez,
+    const amrex::Real Bx, const amrex::Real By, const amrex::Real Bz,
+    const amrex::Real q, const amrex::Real m, const amrex::Real dt )
+{
+    // Constants
+    const amrex::Real econst = q*dt/m;
+    const amrex::Real bconst = 0.5*q*dt/m;
+    constexpr amrex::Real invclight = 1./PhysConst::c;
+    constexpr amrex::Real invclightsq = 1./(PhysConst::c*PhysConst::c);
+    // Compute initial gamma
+    const amrex::Real inv_gamma = 1./std::sqrt(1. + (ux*ux + uy*uy + uz*uz)*invclightsq);
+    // Get tau
+    const amrex::Real taux = bconst*Bx;
+    const amrex::Real tauy = bconst*By;
+    const amrex::Real tauz = bconst*Bz;
+    const amrex::Real tausq = taux*taux+tauy*tauy+tauz*tauz;
+    // Get U', gamma'^2
+    const amrex::Real uxpr = ux + econst*Ex + (uy*tauz-uz*tauy)*inv_gamma;
+    const amrex::Real uypr = uy + econst*Ey + (uz*taux-ux*tauz)*inv_gamma;
+    const amrex::Real uzpr = uz + econst*Ez + (ux*tauy-uy*taux)*inv_gamma;
+    const amrex::Real gprsq = (1. + (uxpr*uxpr + uypr*uypr + uzpr*uzpr)*invclightsq);
+    // Get u*
+    const amrex::Real ust = (uxpr*taux + uypr*tauy + uzpr*tauz)*invclight;
+    // Get new gamma
+    const amrex::Real sigma = gprsq-tausq;
+    const amrex::Real gisq = 2./(sigma + std::sqrt(sigma*sigma + 4.*(tausq + ust*ust)) );
+    // Get t, s
+    const amrex::Real bg = bconst*std::sqrt(gisq);
+    const amrex::Real tx = bg*Bx;
+    const amrex::Real ty = bg*By;
+    const amrex::Real tz = bg*Bz;
+    const amrex::Real s = 1./(1.+tausq*gisq);
+    // Get t.u'
+    const amrex::Real tu = tx*uxpr + ty*uypr + tz*uzpr;
+    // Get new U
+    ux = s*(uxpr+tx*tu+uypr*tz-uzpr*ty);
+    uy = s*(uypr+ty*tu+uzpr*tx-uxpr*tz);
+    uz = s*(uzpr+tz*tu+uxpr*ty-uypr*tx);
+    gaminv = 1./std::sqrt(1. + (ux*ux + uy*uy + uz*uz)*invclightsq);
+}
+
+#endif // WARPX_PARTICLES_PUSHER_UPDATEMOMENTUM_VAY_H_