Merge branch 'dev' into pml_exchange_gpu

1 files changed, 184 insertions, 168 deletions

diff --git a/Source/Laser/LaserParticleContainer.cpp b/Source/Laser/LaserParticleContainer.cpp
index f6ed12fed..15e82e940 100644
--- a/Source/Laser/LaserParticleContainer.cpp
+++ b/Source/Laser/LaserParticleContainer.cpp
@@ -27,128 +27,128 @@ LaserParticleContainer::LaserParticleContainer (AmrCore* amr_core, int ispecies,
     charge = 1.0;
     mass = std::numeric_limits<Real>::max();
     do_boosted_frame_diags = 0;
-        
+
     ParmParse pp(laser_name);
 
-	// Parse the type of laser profile and set the corresponding flag `profile`
-	std::string laser_type_s;
-	pp.get("profile", laser_type_s);
-	std::transform(laser_type_s.begin(), laser_type_s.end(), laser_type_s.begin(), ::tolower);
-	if (laser_type_s == "gaussian") {
-	    profile = laser_t::Gaussian;
+    // Parse the type of laser profile and set the corresponding flag `profile`
+    std::string laser_type_s;
+    pp.get("profile", laser_type_s);
+    std::transform(laser_type_s.begin(), laser_type_s.end(), laser_type_s.begin(), ::tolower);
+    if (laser_type_s == "gaussian") {
+        profile = laser_t::Gaussian;
     } else if(laser_type_s == "harris") {
         profile = laser_t::Harris;
     } else if(laser_type_s == "parse_field_function") {
         profile = laser_t::parse_field_function;
-	} else {
-	    amrex::Abort("Unknown laser type");
-	}
-
-	// Parse the properties of the antenna
-	pp.getarr("position", position);
-	pp.getarr("direction", nvec);
-	pp.getarr("polarization", p_X);
-	pp.query("pusher_algo", pusher_algo);
-	pp.get("wavelength", wavelength);
-	pp.get("e_max", e_max);
-    pp.query("do_continuous_injection", do_continuous_injection);
-
-	if ( profile == laser_t::Gaussian ) {
-	    // Parse the properties of the Gaussian profile
-        pp.get("profile_waist", profile_waist);
-        pp.get("profile_duration", profile_duration);
-        pp.get("profile_t_peak", profile_t_peak);
-        pp.get("profile_focal_distance", profile_focal_distance);
-        stc_direction = p_X;
-        pp.queryarr("stc_direction", stc_direction);
-        pp.query("zeta", zeta);
-        pp.query("beta", beta);
-        pp.query("phi2", phi2);
-	}
-
-    if ( profile == laser_t::Harris ) {
-        // Parse the properties of the Harris profile
-        pp.get("profile_waist", profile_waist);
-        pp.get("profile_duration", profile_duration);
-        pp.get("profile_focal_distance", profile_focal_distance);
-    }
+        } else {
+            amrex::Abort("Unknown laser type");
+        }
 
-    if ( profile == laser_t::parse_field_function ) {
-        // Parse the properties of the parse_field_function profile
-        pp.get("field_function(X,Y,t)", field_function);
-        parser.define(field_function);
-        parser.registerVariables({"X","Y","t"});
-
-        ParmParse ppc("my_constants");
-        std::set<std::string> symbols = parser.symbols();
-        symbols.erase("X");
-        symbols.erase("Y");
-        symbols.erase("t"); // after removing variables, we are left with constants
-        for (auto it = symbols.begin(); it != symbols.end(); ) {
-            Real v;
-            if (ppc.query(it->c_str(), v)) {
-                parser.setConstant(*it, v);
-                it = symbols.erase(it);
-            } else {
-                ++it;
-            }
+        // Parse the properties of the antenna
+        pp.getarr("position", position);
+        pp.getarr("direction", nvec);
+        pp.getarr("polarization", p_X);
+        pp.query("pusher_algo", pusher_algo);
+        pp.get("wavelength", wavelength);
+        pp.get("e_max", e_max);
+        pp.query("do_continuous_injection", do_continuous_injection);
+
+        if ( profile == laser_t::Gaussian ) {
+            // Parse the properties of the Gaussian profile
+            pp.get("profile_waist", profile_waist);
+            pp.get("profile_duration", profile_duration);
+            pp.get("profile_t_peak", profile_t_peak);
+            pp.get("profile_focal_distance", profile_focal_distance);
+            stc_direction = p_X;
+            pp.queryarr("stc_direction", stc_direction);
+            pp.query("zeta", zeta);
+            pp.query("beta", beta);
+            pp.query("phi2", phi2);
         }
-        for (auto const& s : symbols) { // make sure there no unknown symbols
-            amrex::Abort("Laser Profile: Unknown symbol "+s);
+
+        if ( profile == laser_t::Harris ) {
+            // Parse the properties of the Harris profile
+            pp.get("profile_waist", profile_waist);
+            pp.get("profile_duration", profile_duration);
+            pp.get("profile_focal_distance", profile_focal_distance);
         }
-    }
 
-	// Plane normal
-	Real s = 1.0/std::sqrt(nvec[0]*nvec[0] + nvec[1]*nvec[1] + nvec[2]*nvec[2]);
-	nvec = { nvec[0]*s, nvec[1]*s, nvec[2]*s };
-
-    if (WarpX::gamma_boost > 1.) {
-        // Check that the laser direction is equal to the boost direction
-        AMREX_ALWAYS_ASSERT_WITH_MESSAGE(  nvec[0]*WarpX::boost_direction[0]
-                                         + nvec[1]*WarpX::boost_direction[1]
-                                         + nvec[2]*WarpX::boost_direction[2] - 1. < 1.e-12,
-                                           "The Lorentz boost should be in the same direction as the laser propagation");
-        // Get the position of the plane, along the boost direction, in the lab frame
-        // and convert the position of the antenna to the boosted frame
-        Z0_lab = nvec[0]*position[0] + nvec[1]*position[1] + nvec[2]*position[2];
-        Real Z0_boost = Z0_lab/WarpX::gamma_boost;
-        position[0] += (Z0_boost-Z0_lab)*nvec[0];
-        position[1] += (Z0_boost-Z0_lab)*nvec[1];
-        position[2] += (Z0_boost-Z0_lab)*nvec[2];
-    }
+        if ( profile == laser_t::parse_field_function ) {
+            // Parse the properties of the parse_field_function profile
+            pp.get("field_function(X,Y,t)", field_function);
+            parser.define(field_function);
+            parser.registerVariables({"X","Y","t"});
+
+            ParmParse ppc("my_constants");
+            std::set<std::string> symbols = parser.symbols();
+            symbols.erase("X");
+            symbols.erase("Y");
+            symbols.erase("t"); // after removing variables, we are left with constants
+            for (auto it = symbols.begin(); it != symbols.end(); ) {
+                Real v;
+                if (ppc.query(it->c_str(), v)) {
+                    parser.setConstant(*it, v);
+                    it = symbols.erase(it);
+                } else {
+                    ++it;
+                }
+            }
+            for (auto const& s : symbols) { // make sure there no unknown symbols
+                amrex::Abort("Laser Profile: Unknown symbol "+s);
+            }
+        }
+
+        // Plane normal
+        Real s = 1.0/std::sqrt(nvec[0]*nvec[0] + nvec[1]*nvec[1] + nvec[2]*nvec[2]);
+        nvec = { nvec[0]*s, nvec[1]*s, nvec[2]*s };
+
+        if (WarpX::gamma_boost > 1.) {
+            // Check that the laser direction is equal to the boost direction
+            AMREX_ALWAYS_ASSERT_WITH_MESSAGE(  nvec[0]*WarpX::boost_direction[0]
+                                             + nvec[1]*WarpX::boost_direction[1]
+                                             + nvec[2]*WarpX::boost_direction[2] - 1. < 1.e-12,
+                                               "The Lorentz boost should be in the same direction as the laser propagation");
+            // Get the position of the plane, along the boost direction, in the lab frame
+            // and convert the position of the antenna to the boosted frame
+            Z0_lab = nvec[0]*position[0] + nvec[1]*position[1] + nvec[2]*position[2];
+            Real Z0_boost = Z0_lab/WarpX::gamma_boost;
+            position[0] += (Z0_boost-Z0_lab)*nvec[0];
+            position[1] += (Z0_boost-Z0_lab)*nvec[1];
+            position[2] += (Z0_boost-Z0_lab)*nvec[2];
+        }
 
-	// The first polarization vector
-	s = 1.0/std::sqrt(p_X[0]*p_X[0] + p_X[1]*p_X[1] + p_X[2]*p_X[2]);
-	p_X = { p_X[0]*s, p_X[1]*s, p_X[2]*s };
+        // The first polarization vector
+        s = 1.0/std::sqrt(p_X[0]*p_X[0] + p_X[1]*p_X[1] + p_X[2]*p_X[2]);
+        p_X = { p_X[0]*s, p_X[1]*s, p_X[2]*s };
 
-	Real dp = std::inner_product(nvec.begin(), nvec.end(), p_X.begin(), 0.0);
-    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(std::abs(dp) < 1.0e-14,
-                                     "Laser plane vector is not perpendicular to the main polarization vector");
+        Real dp = std::inner_product(nvec.begin(), nvec.end(), p_X.begin(), 0.0);
+        AMREX_ALWAYS_ASSERT_WITH_MESSAGE(std::abs(dp) < 1.0e-14,
+            "Laser plane vector is not perpendicular to the main polarization vector");
 
-	p_Y = CrossProduct(nvec, p_X);   // The second polarization vector
+        p_Y = CrossProduct(nvec, p_X);   // The second polarization vector
 
-	s = 1.0/std::sqrt(stc_direction[0]*stc_direction[0] + stc_direction[1]*stc_direction[1] + stc_direction[2]*stc_direction[2]);
-	stc_direction = { stc_direction[0]*s, stc_direction[1]*s, stc_direction[2]*s };
-	dp = std::inner_product(nvec.begin(), nvec.end(), stc_direction.begin(), 0.0);
-    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(std::abs(dp) < 1.0e-14,
-                                     "stc_direction is not perpendicular to the laser plane vector");
+        s = 1.0/std::sqrt(stc_direction[0]*stc_direction[0] + stc_direction[1]*stc_direction[1] + stc_direction[2]*stc_direction[2]);
+        stc_direction = { stc_direction[0]*s, stc_direction[1]*s, stc_direction[2]*s };
+        dp = std::inner_product(nvec.begin(), nvec.end(), stc_direction.begin(), 0.0);
+        AMREX_ALWAYS_ASSERT_WITH_MESSAGE(std::abs(dp) < 1.0e-14,
+            "stc_direction is not perpendicular to the laser plane vector");
 
-	// Get angle between p_X and stc_direction
-	// in 2d, stcs are in the simulation plane
+        // Get angle between p_X and stc_direction
+        // in 2d, stcs are in the simulation plane
 #if AMREX_SPACEDIM == 3
-	theta_stc = acos(stc_direction[0]*p_X[0] +
+        theta_stc = acos(stc_direction[0]*p_X[0] +
                      stc_direction[1]*p_X[1] +
                      stc_direction[2]*p_X[2]);
 #else
-	theta_stc = 0.;
+        theta_stc = 0.;
 #endif
 
 #if AMREX_SPACEDIM == 3
-	u_X = p_X;
-	u_Y = p_Y;
+        u_X = p_X;
+        u_Y = p_Y;
 #else
-	u_X = CrossProduct({0., 1., 0.}, nvec);
-	u_Y = {0., 1., 0.};
+        u_X = CrossProduct({0., 1., 0.}, nvec);
+        u_Y = {0., 1., 0.};
 #endif
 
     laser_injection_box= Geom(0).ProbDomain();
@@ -166,7 +166,7 @@ LaserParticleContainer::LaserParticleContainer (AmrCore* amr_core, int ispecies,
         // If laser antenna initially outside of the box, store its theoretical
         // position in z_antenna_th
         updated_position = position;
-        
+
         // Sanity checks
         int dir = WarpX::moving_window_dir;
         std::vector<Real> windir(3, 0.0);
@@ -175,16 +175,16 @@ LaserParticleContainer::LaserParticleContainer (AmrCore* amr_core, int ispecies,
 #else
         windir[dir] = 1.0;
 #endif
-        AMREX_ALWAYS_ASSERT_WITH_MESSAGE(  
-            (nvec[0]-windir[0]) + (nvec[1]-windir[1]) + (nvec[2]-windir[2]) 
+        AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
+            (nvec[0]-windir[0]) + (nvec[1]-windir[1]) + (nvec[2]-windir[2])
             < 1.e-12, "do_continous_injection for laser particle only works" +
             " if moving window direction and laser propagation direction are the same");
         if ( WarpX::gamma_boost>1 ){
             AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
-                (WarpX::boost_direction[0]-0)*(WarpX::boost_direction[0]-0) + 
-                (WarpX::boost_direction[1]-0)*(WarpX::boost_direction[1]-0) + 
+                (WarpX::boost_direction[0]-0)*(WarpX::boost_direction[0]-0) +
+                (WarpX::boost_direction[1]-0)*(WarpX::boost_direction[1]-0) +
                 (WarpX::boost_direction[2]-1)*(WarpX::boost_direction[2]-1) < 1.e-12,
-                "do_continous_injection for laser particle only works if " + 
+                "do_continous_injection for laser particle only works if " +
                 "warpx.boost_direction = z. TODO: all directions.");
         }
     }
@@ -192,13 +192,13 @@ LaserParticleContainer::LaserParticleContainer (AmrCore* amr_core, int ispecies,
 
 /* \brief Check if laser particles enter the box, and inject if necessary.
  * \param injection_box: a RealBox where particles should be injected.
- */ 
+ */
 void
 LaserParticleContainer::ContinuousInjection (const RealBox& injection_box)
 {
     // Input parameter injection_box contains small box where injection
     // should occur.
-    // So far, LaserParticleContainer::laser_injection_box contains the 
+    // So far, LaserParticleContainer::laser_injection_box contains the
     // outdated full problem domain at t=0.
 
     // Convert updated_position to Real* to use RealBox::contains().
@@ -258,7 +258,7 @@ LaserParticleContainer::InitData (int lev)
     ComputeSpacing(lev, S_X, S_Y);
     ComputeWeightMobility(S_X, S_Y);
 
-    // LaserParticleContainer::position contains the initial position of the 
+    // LaserParticleContainer::position contains the initial position of the
     // laser antenna. In the boosted frame, the antenna is moving.
     // Update its position with updated_position.
     if (do_continuous_injection){
@@ -327,19 +327,19 @@ LaserParticleContainer::InitData (int lev)
                          IntVect(plane_hi[0],plane_hi[1],0)};
     BoxArray plane_ba {plane_box};
     {
-	IntVect chunk(plane_box.size());
-	const int min_size = 8;
-	while (plane_ba.size() < nprocs && chunk[0] > min_size && chunk[1] > min_size)
-	{
-	    for (int j = 1; j >= 0 ; j--)
-	    {
-		chunk[j] /= 2;
-
-		if (plane_ba.size() < nprocs) {
-		    plane_ba.maxSize(chunk);
-		}
-	    }
-	}
+        IntVect chunk(plane_box.size());
+        const int min_size = 8;
+        while (plane_ba.size() < nprocs && chunk[0] > min_size && chunk[1] > min_size)
+        {
+            for (int j = 1; j >= 0 ; j--)
+            {
+                chunk[j] /= 2;
+
+                if (plane_ba.size() < nprocs) {
+                    plane_ba.maxSize(chunk);
+                }
+            }
+        }
     }
 #else
     BoxArray plane_ba { Box {IntVect(plane_lo[0],0), IntVect(plane_hi[0],0)} };
@@ -351,29 +351,29 @@ LaserParticleContainer::InitData (int lev)
     const Vector<int>& procmap = plane_dm.ProcessorMap();
     for (int i = 0, n = plane_ba.size(); i < n; ++i)
     {
-	if (procmap[i] == myproc)
-	{
-	    const Box& bx = plane_ba[i];
-	    for (IntVect cell = bx.smallEnd(); cell <= bx.bigEnd(); bx.next(cell))
-	    {
-            const Vector<Real>& pos = Transform(cell[0], cell[1]);
+        if (procmap[i] == myproc)
+        {
+            const Box& bx = plane_ba[i];
+            for (IntVect cell = bx.smallEnd(); cell <= bx.bigEnd(); bx.next(cell))
+            {
+                const Vector<Real>& pos = Transform(cell[0], cell[1]);
 #if (AMREX_SPACEDIM == 3)
-            const Real* x = pos.dataPtr();
+                const Real* x = pos.dataPtr();
 #else
-            const Real x[2] = {pos[0], pos[2]};
+                const Real x[2] = {pos[0], pos[2]};
 #endif
-            if (laser_injection_box.contains(x))
-            {
-                for (int k = 0; k<2; ++k) {
-                    particle_x.push_back(pos[0]);
-                    particle_y.push_back(pos[1]);
-                    particle_z.push_back(pos[2]);
+                if (laser_injection_box.contains(x))
+                {
+                    for (int k = 0; k<2; ++k) {
+                        particle_x.push_back(pos[0]);
+                        particle_y.push_back(pos[1]);
+                        particle_z.push_back(pos[2]);
+                    }
+                    particle_w.push_back( weight);
+                    particle_w.push_back(-weight);
                 }
-                particle_w.push_back( weight);
-                particle_w.push_back(-weight);
             }
-	    }
-	}
+        }
     }
     const int np = particle_z.size();
     RealVector particle_ux(np, 0.0);
@@ -459,7 +459,7 @@ LaserParticleContainer::Evolve (int lev,
 
             if (rho) {
                 int* AMREX_RESTRICT ion_lev = nullptr;
-                DepositCharge(pti, wp, ion_lev, rho, 0, 0, 
+                DepositCharge(pti, wp, ion_lev, rho, 0, 0,
                               np_current, thread_num, lev, lev);
                 if (crho) {
                     DepositCharge(pti, wp, ion_lev, crho, 0, np_current,
@@ -493,6 +493,7 @@ LaserParticleContainer::Evolve (int lev,
                     amplitude_E[i] = parser.eval(plane_Xp[i], plane_Yp[i], t);
                 }
             }
+
             // Calculate the corresponding momentum and position for the particles
             update_laser_particle(np, uxp.dataPtr(), uyp.dataPtr(),
                                   uzp.dataPtr(), wp.dataPtr(),
@@ -507,6 +508,7 @@ LaserParticleContainer::Evolve (int lev,
             DepositCurrent(pti, wp, uxp, uyp, uzp, ion_lev, &jx, &jy, &jz,
                            0, np_current, thread_num,
                            lev, lev, dt);
+
             bool has_buffer = cjx;
             if (has_buffer){
                 // Deposit in buffers
@@ -604,36 +606,45 @@ LaserParticleContainer::PushP (int lev, Real dt,
  * \param np: number of laser particles
  * \param thread_num: thread number
  * \param pplane_Xp, pplane_Yp: pointers to arrays of particle positions
- * in laser plane coordinate. 
+ * in laser plane coordinate.
  */
 void
 LaserParticleContainer::calculate_laser_plane_coordinates (
     const int np, const int thread_num,
-    Real * const pplane_Xp,
-    Real * const pplane_Yp)
+    Real * AMREX_RESTRICT const pplane_Xp,
+    Real * AMREX_RESTRICT const pplane_Yp)
 {
     Real const * const AMREX_RESTRICT xp = m_xp[thread_num].dataPtr();
     Real const * const AMREX_RESTRICT yp = m_yp[thread_num].dataPtr();
     Real const * const AMREX_RESTRICT zp = m_zp[thread_num].dataPtr();
-    Real const * const AMREX_RESTRICT tmp_u_X = u_X.dataPtr();
-    Real const * const AMREX_RESTRICT tmp_u_Y = u_Y.dataPtr();
-    Real const * const AMREX_RESTRICT tmp_position = position.dataPtr();
+    Real tmp_u_X_0 = u_X[0];
+    Real tmp_u_X_2 = u_X[2];
+    Real tmp_position_0 = position[0];
+    Real tmp_position_2 = position[2];
+#if (AMREX_SPACEDIM == 3)
+    Real tmp_u_X_1 = u_X[1];
+    Real tmp_u_Y_0 = u_Y[0];
+    Real tmp_u_Y_1 = u_Y[1];
+    Real tmp_u_Y_2 = u_Y[2];
+    Real tmp_position_1 = position[1];
+#endif
+
     amrex::ParallelFor(
         np,
         [=] AMREX_GPU_DEVICE (int i) {
 #if (AMREX_SPACEDIM == 3)
-            pplane_Xp[i] = 
-                tmp_u_X[0] * (xp[i] - tmp_position[0]) +
-                tmp_u_X[1] * (yp[i] - tmp_position[1]) +
-                tmp_u_X[2] * (zp[i] - tmp_position[2]);
-            pplane_Yp[i] = 
-                tmp_u_Y[0] * (xp[i] - tmp_position[0]) +
-                tmp_u_Y[1] * (yp[i] - tmp_position[1]) +
-                tmp_u_Y[2] * (zp[i] - tmp_position[2]);
+            pplane_Xp[i] =
+                tmp_u_X_0 * (xp[i] - tmp_position_0) +
+                tmp_u_X_1 * (yp[i] - tmp_position_1) +
+                tmp_u_X_2 * (zp[i] - tmp_position_2);
+            pplane_Yp[i] =
+                tmp_u_Y_0 * (xp[i] - tmp_position_0) +
+                tmp_u_Y_1 * (yp[i] - tmp_position_1) +
+                tmp_u_Y_2 * (zp[i] - tmp_position_2);
 #elif (AMREX_SPACEDIM == 2)
-            pplane_Xp[i] = 
-                tmp_u_X[0] * (xp[i] - tmp_position[0]) +
-                tmp_u_X[2] * (zp[i] - tmp_position[2]);
+            pplane_Xp[i] =
+                tmp_u_X_0 * (xp[i] - tmp_position_0) +
+                tmp_u_X_2 * (zp[i] - tmp_position_2);
             pplane_Yp[i] = 0.;
 #endif
         }
@@ -651,17 +662,22 @@ LaserParticleContainer::calculate_laser_plane_coordinates (
  */
 void
 LaserParticleContainer::update_laser_particle(
-    const int np, Real * const puxp, Real * const puyp, Real * const puzp,
-    Real * const pwp, Real const * const amplitude, const Real dt,
+    const int np, Real * AMREX_RESTRICT const puxp, Real * AMREX_RESTRICT const puyp,
+    Real * AMREX_RESTRICT const puzp, Real const * AMREX_RESTRICT const pwp,
+    Real const * AMREX_RESTRICT const amplitude, const Real dt,
     const int thread_num)
 {
     Real * const AMREX_RESTRICT xp = m_xp[thread_num].dataPtr();
     Real * const AMREX_RESTRICT yp = m_yp[thread_num].dataPtr();
     Real * const AMREX_RESTRICT zp = m_zp[thread_num].dataPtr();
     Real * const AMREX_RESTRICT giv = m_giv[thread_num].dataPtr();
-    Real const * const AMREX_RESTRICT tmp_p_X = p_X.dataPtr();
-    Real const * const AMREX_RESTRICT tmp_nvec = nvec.dataPtr();
-    
+    Real tmp_p_X_0 = p_X[0];
+    Real tmp_p_X_1 = p_X[1];
+    Real tmp_p_X_2 = p_X[2];
+    Real tmp_nvec_0 = nvec[0];
+    Real tmp_nvec_1 = nvec[1];
+    Real tmp_nvec_2 = nvec[2];
+
     // Copy member variables to tmp copies for GPU runs.
     Real tmp_mobility = mobility;
     Real gamma_boost = WarpX::gamma_boost;
@@ -673,14 +689,14 @@ LaserParticleContainer::update_laser_particle(
             const Real sign_charge = (pwp[i]>0) ? 1 : -1;
             const Real v_over_c = sign_charge * tmp_mobility * amplitude[i];
             // The velocity is along the laser polarization p_X
-            Real vx = PhysConst::c * v_over_c * tmp_p_X[0];
-            Real vy = PhysConst::c * v_over_c * tmp_p_X[1];
-            Real vz = PhysConst::c * v_over_c * tmp_p_X[2];
+            Real vx = PhysConst::c * v_over_c * tmp_p_X_0;
+            Real vy = PhysConst::c * v_over_c * tmp_p_X_1;
+            Real vz = PhysConst::c * v_over_c * tmp_p_X_2;
             // When running in the boosted-frame, their is additional velocity along nvec
             if (gamma_boost > 1.){
-                vx -= PhysConst::c * beta_boost * tmp_nvec[0];
-                vy -= PhysConst::c * beta_boost * tmp_nvec[1];
-                vz -= PhysConst::c * beta_boost * tmp_nvec[2];
+                vx -= PhysConst::c * beta_boost * tmp_nvec_0;
+                vy -= PhysConst::c * beta_boost * tmp_nvec_1;
+                vz -= PhysConst::c * beta_boost * tmp_nvec_2;
             }
             // Get the corresponding momenta
             const Real gamma = gamma_boost/std::sqrt(1. - v_over_c*v_over_c);