Merge branch 'classical_radiation_reaction' of github.com:lucafedeli88/WarpX into classical_radiation_reaction

1 files changed, 29 insertions, 13 deletions

diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp
index 95f2243b3..3d727af84 100644
--- a/Source/Particles/PhysicalParticleContainer.cpp
+++ b/Source/Particles/PhysicalParticleContainer.cpp
@@ -943,7 +943,7 @@ PhysicalParticleContainer::Evolve (int lev,
                                    MultiFab* rho, MultiFab* crho,
                                    const MultiFab* cEx, const MultiFab* cEy, const MultiFab* cEz,
                                    const MultiFab* cBx, const MultiFab* cBy, const MultiFab* cBz,
-                                   Real t, Real dt)
+                                   Real t, Real dt, DtType a_dt_type)
 {
     BL_PROFILE("PPC::Evolve()");
     BL_PROFILE_VAR_NS("PPC::Evolve::Copy", blp_copy);
@@ -1302,7 +1302,7 @@ PhysicalParticleContainer::Evolve (int lev,
                 // Particle Push
                 //
                 BL_PROFILE_VAR_START(blp_ppc_pp);
-                PushPX(pti, m_xp[thread_num], m_yp[thread_num], m_zp[thread_num], dt);
+                PushPX(pti, m_xp[thread_num], m_yp[thread_num], m_zp[thread_num], dt, a_dt_type);
                 BL_PROFILE_VAR_STOP(blp_ppc_pp);
 
                 //
@@ -1364,8 +1364,16 @@ PhysicalParticleContainer::Evolve (int lev,
             }
         }
     }
-    // Split particles
-    if (do_splitting){ SplitParticles(lev); }
+    // Split particles at the end of the timestep.
+    // When subcycling is ON, the splitting is done on the last call to
+    // PhysicalParticleContainer::Evolve on the finest level, i.e., at the
+    // end of the large timestep. Otherwise, the pushes on different levels
+    // are not consistent, and the call to Redistribute (inside
+    // SplitParticles) may result in split particles to deposit twice on the
+    // coarse level.
+    if (do_splitting && (a_dt_type == DtType::SecondHalf || a_dt_type == DtType::Full) ){
+        SplitParticles(lev);
+    }
 }
 
 // Loop over all particles in the particle container and
@@ -1409,7 +1417,7 @@ PhysicalParticleContainer::SplitParticles(int lev)
                 np_split_to_add += np_split;
 #if (AMREX_SPACEDIM==2)
                 if (split_type==0){
-                    // Split particle in two along each axis
+                    // Split particle in two along each diagonals
                     // 4 particles in 2d
                     for (int ishift = -1; ishift < 2; ishift +=2 ){
                         for (int kshift = -1; kshift < 2; kshift +=2 ){
@@ -1424,7 +1432,7 @@ PhysicalParticleContainer::SplitParticles(int lev)
                         }
                     }
                 } else {
-                    // Split particle in two along each diagonal
+                    // Split particle in two along each axis
                     // 4 particles in 2d
                     for (int ishift = -1; ishift < 2; ishift +=2 ){
                         // Add one particle with offset in x
@@ -1447,15 +1455,15 @@ PhysicalParticleContainer::SplitParticles(int lev)
                 }
 #elif (AMREX_SPACEDIM==3)
                 if (split_type==0){
-                    // Split particle in two along each axis
-                    // 6 particles in 2d
+                    // Split particle in two along each diagonals
+                    // 8 particles in 3d
                     for (int ishift = -1; ishift < 2; ishift +=2 ){
                         for (int jshift = -1; jshift < 2; jshift +=2 ){
                             for (int kshift = -1; kshift < 2; kshift +=2 ){
                                 // Add one particle with offset in x, y and z
                                 psplit_x.push_back( xp[i] + ishift*dx[0]/2 );
                                 psplit_y.push_back( yp[i] + jshift*dx[1]/2 );
-                                psplit_z.push_back( zp[i] + jshift*dx[2]/2 );
+                                psplit_z.push_back( zp[i] + kshift*dx[2]/2 );
                                 psplit_ux.push_back( uxp[i] );
                                 psplit_uy.push_back( uyp[i] );
                                 psplit_uz.push_back( uzp[i] );
@@ -1464,8 +1472,8 @@ PhysicalParticleContainer::SplitParticles(int lev)
                         }
                     }
                 } else {
-                    // Split particle in two along each diagonal
-                    // 8 particles in 3d
+                    // Split particle in two along each axis
+                    // 6 particles in 3d
                     for (int ishift = -1; ishift < 2; ishift +=2 ){
                         // Add one particle with offset in x
                         psplit_x.push_back( xp[i] + ishift*dx[0]/2 );
@@ -1527,7 +1535,7 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti,
                                   Cuda::ManagedDeviceVector<Real>& xp,
                                   Cuda::ManagedDeviceVector<Real>& yp,
                                   Cuda::ManagedDeviceVector<Real>& zp,
-                                  Real dt)
+                                  Real dt, DtType a_dt_type)
 {
 
     // This wraps the momentum and position advance so that inheritors can modify the call.
@@ -1546,7 +1554,7 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti,
     const Real* const AMREX_RESTRICT By = attribs[PIdx::By].dataPtr();
     const Real* const AMREX_RESTRICT Bz = attribs[PIdx::Bz].dataPtr();
 
-    if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags)
+    if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags && (a_dt_type!=DtType::SecondHalf))
     {
         copy_attribs(pti, x, y, z);
     }
@@ -1559,6 +1567,7 @@ 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){
         amrex::ParallelFor(
             pti.numParticles(),
@@ -1693,6 +1702,13 @@ PhysicalParticleContainer::PushP (int lev, Real dt,
                               Expp[i], Eypp[i], Ezpp[i], Bxpp[i], Bypp[i], Bzpp[i], q, m, dt);
                     }
                 );
+            } else if (WarpX::particle_pusher_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");
             };