Merge branch 'dev' into generalize_nodal_deposition

1 files changed, 267 insertions, 343 deletions

diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp
index 2dc25e6fa..02dee1967 100644
--- a/Source/Particles/PhysicalParticleContainer.cpp
+++ b/Source/Particles/PhysicalParticleContainer.cpp
@@ -31,10 +31,13 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp
 
     pp.query("boost_adjust_transverse_positions", boost_adjust_transverse_positions);
     pp.query("do_backward_propagation", do_backward_propagation);
+
+    // Initialize splitting
     pp.query("do_splitting", do_splitting);
     pp.query("split_type", split_type);
+
     pp.query("do_continuous_injection", do_continuous_injection);
-    // Whether to plot back-transformed (lab-frame) diagnostics 
+    // Whether to plot back-transformed (lab-frame) diagnostics
     // for this species.
     pp.query("do_boosted_frame_diags", do_boosted_frame_diags);
 
@@ -51,7 +54,7 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp
     do_user_plot_vars = pp.queryarr("plot_vars", plot_vars);
     if (not do_user_plot_vars){
         // By default, all particle variables are dumped to plotfiles,
-        // including {x,y,z,ux,uy,uz}old variables when running in a 
+        // including {x,y,z,ux,uy,uz}old variables when running in a
         // boosted frame
         if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags){
             plot_flags.resize(PIdx::nattribs + 6, 1);
@@ -68,9 +71,9 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp
         // If not none, set plot_flags values to 1 for elements in plot_vars.
         if (plot_vars[0] != "none"){
             for (const auto& var : plot_vars){
-                // Return error if var not in PIdx. 
-                AMREX_ALWAYS_ASSERT_WITH_MESSAGE( 
-                    ParticleStringNames::to_index.count(var), 
+                // Return error if var not in PIdx.
+                AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
+                    ParticleStringNames::to_index.count(var),
                     "plot_vars argument not in ParticleStringNames");
                 plot_flags[ParticleStringNames::to_index.at(var)] = 1;
             }
@@ -142,7 +145,7 @@ void PhysicalParticleContainer::MapParticletoBoostedFrame(Real& x, Real& y, Real
 void
 PhysicalParticleContainer::AddGaussianBeam(Real x_m, Real y_m, Real z_m,
                                            Real x_rms, Real y_rms, Real z_rms,
-                                           Real q_tot, long npart, 
+                                           Real q_tot, long npart,
                                            int do_symmetrize) {
 
     const Geometry& geom     = m_gdb->Geom(0);
@@ -154,7 +157,7 @@ PhysicalParticleContainer::AddGaussianBeam(Real x_m, Real y_m, Real z_m,
     std::normal_distribution<double> distz(z_m, z_rms);
 
     if (ParallelDescriptor::IOProcessor()) {
-        // If do_symmetrize, create 4x fewer particles, and 
+        // If do_symmetrize, create 4x fewer particles, and
         // Replicate each particle 4 times (x,y) (-x,y) (x,-y) (-x,-y)
         if (do_symmetrize){
             npart /= 4;
@@ -196,7 +199,7 @@ PhysicalParticleContainer::CheckAndAddParticle(Real x, Real y, Real z,
                                                std::array<Real, 3> u,
                                                Real weight)
 {
-    std::array<Real,PIdx::nattribs> attribs;
+    std::array<ParticleReal,PIdx::nattribs> attribs;
     attribs.fill(0.0);
 
     // update attribs with input arguments
@@ -361,13 +364,13 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox)
         for (int dir=0; dir<AMREX_SPACEDIM; dir++) {
             if ( tile_realbox.lo(dir) <= part_realbox.hi(dir) ) {
                 Real ncells_adjust = std::floor( (tile_realbox.lo(dir) - part_realbox.lo(dir))/dx[dir] );
-                overlap_realbox.setLo( dir, part_realbox.lo(dir) + std::max(ncells_adjust, 0.) * dx[dir]);
+                overlap_realbox.setLo( dir, part_realbox.lo(dir) + std::max(ncells_adjust, Real(0.)) * dx[dir]);
             } else {
                 no_overlap = true; break;
             }
             if ( tile_realbox.hi(dir) >= part_realbox.lo(dir) ) {
                 Real ncells_adjust = std::floor( (part_realbox.hi(dir) - tile_realbox.hi(dir))/dx[dir] );
-                overlap_realbox.setHi( dir, part_realbox.hi(dir) - std::max(ncells_adjust, 0.) * dx[dir]);
+                overlap_realbox.setHi( dir, part_realbox.hi(dir) - std::max(ncells_adjust, Real(0.)) * dx[dir]);
             } else {
                 no_overlap = true; break;
             }
@@ -387,11 +390,11 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox)
         const int tile_id = mfi.LocalTileIndex();
 
         // Max number of new particles, if particles are created in the whole
-        // overlap_box. All of them are created, and invalid ones are then 
+        // overlap_box. All of them are created, and invalid ones are then
         // discaded
         int max_new_particles = overlap_box.numPts() * num_ppc;
 
-        // If refine injection, build pointer dp_cellid that holds pointer to 
+        // If refine injection, build pointer dp_cellid that holds pointer to
         // array of refined cell IDs.
         Vector<int> cellid_v;
         if (refine_injection and lev == 0)
@@ -437,7 +440,7 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox)
 
         ParticleType* pp = particle_tile.GetArrayOfStructs()().data() + old_size;
         auto& soa = particle_tile.GetStructOfArrays();
-        GpuArray<Real*,PIdx::nattribs> pa;
+        GpuArray<ParticleReal*,PIdx::nattribs> pa;
         for (int ia = 0; ia < PIdx::nattribs; ++ia) {
             pa[ia] = soa.GetRealData(ia).data() + old_size;
         }
@@ -446,7 +449,7 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox)
         if (do_field_ionization) {
             pi = soa.GetIntData(particle_icomps["ionization_level"]).data() + old_size;
         }
-        
+
         const GpuArray<Real,AMREX_SPACEDIM> overlap_corner
             {AMREX_D_DECL(overlap_realbox.lo(0),
                           overlap_realbox.lo(1),
@@ -455,12 +458,12 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox)
         std::size_t shared_mem_bytes = plasma_injector->sharedMemoryNeeded();
         int lrrfac = rrfac;
 
-	bool loc_do_field_ionization = do_field_ionization;
-	int loc_ionization_initial_level = ionization_initial_level;
+        bool loc_do_field_ionization = do_field_ionization;
+        int loc_ionization_initial_level = ionization_initial_level;
 
-        // Loop over all new particles and inject them (creates too many 
+        // Loop over all new particles and inject them (creates too many
         // particles, in particular does not consider xmin, xmax etc.).
-        // The invalid ones are given negative ID and are deleted during the 
+        // The invalid ones are given negative ID and are deleted during the
         // next redistribute.
         amrex::For(max_new_particles, [=] AMREX_GPU_DEVICE (int ip) noexcept
         {
@@ -627,7 +630,7 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox)
             p.pos(1) = z;
 #endif
         }, shared_mem_bytes);
-    			 
+
         if (cost) {
             wt = (amrex::second() - wt) / tile_box.d_numPts();
             Array4<Real> const& costarr = cost->array(mfi);
@@ -865,7 +868,7 @@ PhysicalParticleContainer::FieldGather (int lev,
     MultiFab* cost = WarpX::getCosts(lev);
 
 #ifdef _OPENMP
-#pragma omp parallel 
+#pragma omp parallel
 #endif
     {
 #ifdef _OPENMP
@@ -915,7 +918,7 @@ PhysicalParticleContainer::FieldGather (int lev,
             //
             int e_is_nodal = Ex.is_nodal() and Ey.is_nodal() and Ez.is_nodal();
             FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
-                        &exfab, &eyfab, &ezfab, &bxfab, &byfab, &bzfab, 
+                        &exfab, &eyfab, &ezfab, &bxfab, &byfab, &bzfab,
                         Ex.nGrow(), e_is_nodal,
                         0, np, thread_num, lev, lev);
 
@@ -942,25 +945,23 @@ 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);
     BL_PROFILE_VAR_NS("PPC::FieldGather", blp_fg);
     BL_PROFILE_VAR_NS("PPC::ParticlePush", blp_ppc_pp);
-    BL_PROFILE_VAR_NS("PPC::Evolve::partition", blp_partition);
-    
+
     const std::array<Real,3>& dx = WarpX::CellSize(lev);
     const std::array<Real,3>& cdx = WarpX::CellSize(std::max(lev-1,0));
 
-    // Get instances of NCI Godfrey filters 
+    // Get instances of NCI Godfrey filters
     const auto& nci_godfrey_filter_exeybz = WarpX::GetInstance().nci_godfrey_filter_exeybz;
     const auto& nci_godfrey_filter_bxbyez = WarpX::GetInstance().nci_godfrey_filter_bxbyez;
 
     BL_ASSERT(OnSameGrids(lev,jx));
 
     MultiFab* cost = WarpX::getCosts(lev);
-
     const iMultiFab* current_masks = WarpX::CurrentBufferMasks(lev);
     const iMultiFab* gather_masks = WarpX::GatherBufferMasks(lev);
 
@@ -978,9 +979,9 @@ PhysicalParticleContainer::Evolve (int lev,
                 tmp_particle_data[lev][index][i].resize(np);
         }
     }
-    
+
 #ifdef _OPENMP
-#pragma omp parallel 
+#pragma omp parallel
 #endif
     {
 #ifdef _OPENMP
@@ -991,10 +992,6 @@ PhysicalParticleContainer::Evolve (int lev,
 
         FArrayBox filtered_Ex, filtered_Ey, filtered_Ez;
         FArrayBox filtered_Bx, filtered_By, filtered_Bz;
-        std::vector<bool> inexflag;
-        Vector<long> pid;
-        RealVector tmp;
-        ParticleVector particle_tmp;
 
         for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti)
         {
@@ -1026,56 +1023,18 @@ PhysicalParticleContainer::Evolve (int lev,
             FArrayBox const* bzfab = &(Bz[pti]);
 
             Elixir exeli, eyeli, ezeli, bxeli, byeli, bzeli;
+
             if (WarpX::use_fdtd_nci_corr)
             {
-#if (AMREX_SPACEDIM == 2)
-                const Box& tbox = amrex::grow(pti.tilebox(),{static_cast<int>(WarpX::nox),
-                            static_cast<int>(WarpX::noz)});
-#else
-                const Box& tbox = amrex::grow(pti.tilebox(),{static_cast<int>(WarpX::nox),
-                            static_cast<int>(WarpX::noy),
-                            static_cast<int>(WarpX::noz)});
-#endif
-
-                // Filter Ex (Both 2D and 3D)
-                filtered_Ex.resize(amrex::convert(tbox,WarpX::Ex_nodal_flag));
-                // Safeguard for GPU
-                exeli = filtered_Ex.elixir();
-                // Apply filter on Ex, result stored in filtered_Ex
-                nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Ex, Ex[pti], filtered_Ex.box());
-                // Update exfab reference
-                exfab = &filtered_Ex;
-
-                // Filter Ez
-                filtered_Ez.resize(amrex::convert(tbox,WarpX::Ez_nodal_flag));
-                ezeli = filtered_Ez.elixir();
-                nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_Ez, Ez[pti], filtered_Ez.box());
-                ezfab = &filtered_Ez;
-
-                // Filter By
-                filtered_By.resize(amrex::convert(tbox,WarpX::By_nodal_flag));
-                byeli = filtered_By.elixir();
-                nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_By, By[pti], filtered_By.box());
-                byfab = &filtered_By;
-#if (AMREX_SPACEDIM == 3)
-                // Filter Ey
-                filtered_Ey.resize(amrex::convert(tbox,WarpX::Ey_nodal_flag));
-                eyeli = filtered_Ey.elixir();
-                nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Ey, Ey[pti], filtered_Ey.box());
-                eyfab = &filtered_Ey;
-
-                // Filter Bx
-                filtered_Bx.resize(amrex::convert(tbox,WarpX::Bx_nodal_flag));
-                bxeli = filtered_Bx.elixir();
-                nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_Bx, Bx[pti], filtered_Bx.box());
-                bxfab = &filtered_Bx;
-
-                // Filter Bz
-                filtered_Bz.resize(amrex::convert(tbox,WarpX::Bz_nodal_flag));
-                bzeli = filtered_Bz.elixir();
-                nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Bz, Bz[pti], filtered_Bz.box());
-                bzfab = &filtered_Bz;
-#endif
+                // Filter arrays Ex[pti], store the result in
+                // filtered_Ex and update pointer exfab so that it
+                // points to filtered_Ex (and do the same for all
+                // components of E and B).
+                applyNCIFilter(lev, pti.tilebox(), exeli, eyeli, ezeli, bxeli, byeli, bzeli,
+                               filtered_Ex, filtered_Ey, filtered_Ez,
+                               filtered_Bx, filtered_By, filtered_Bz,
+                               Ex[pti], Ey[pti], Ez[pti], Bx[pti], By[pti], Bz[pti],
+                               exfab, eyfab, ezfab, bxfab, byfab, bzfab);
             }
 
             Exp.assign(np,0.0);
@@ -1085,101 +1044,25 @@ PhysicalParticleContainer::Evolve (int lev,
             Byp.assign(np,WarpX::B_external[1]);
             Bzp.assign(np,WarpX::B_external[2]);
 
-            m_giv[thread_num].resize(np);
-
+            // Determine which particles deposit/gather in the buffer, and
+            // which particles deposit/gather in the fine patch
             long nfine_current = np;
             long nfine_gather = np;
-            if (has_buffer && !do_not_push)
-            {
-                BL_PROFILE_VAR_START(blp_partition);
-                inexflag.resize(np);
-                auto& aos = pti.GetArrayOfStructs();
-                // We need to partition the large buffer first
-                iMultiFab const* bmasks = (WarpX::n_field_gather_buffer >= WarpX::n_current_deposition_buffer) ?
-                    gather_masks : current_masks;
-                int i = 0;
-                const auto& msk = (*bmasks)[pti];
-                for (const auto& p : aos) {
-                    const IntVect& iv = Index(p, lev);
-                    inexflag[i++] = msk(iv);
-                }
-
-                pid.resize(np);
-                std::iota(pid.begin(), pid.end(), 0L);
-
-                auto sep = std::stable_partition(pid.begin(), pid.end(),
-                                                 [&inexflag](long id) { return inexflag[id]; });
-
-                if (WarpX::n_current_deposition_buffer == WarpX::n_field_gather_buffer) {
-                    nfine_current = nfine_gather = std::distance(pid.begin(), sep);
-                } else if (sep != pid.end()) {
-                    int n_buf;
-                    if (bmasks == gather_masks) {
-                        nfine_gather = std::distance(pid.begin(), sep);
-                        bmasks = current_masks;
-                        n_buf = WarpX::n_current_deposition_buffer;
-                    } else {
-                        nfine_current = std::distance(pid.begin(), sep);
-                        bmasks = gather_masks;
-                        n_buf = WarpX::n_field_gather_buffer;
-                    }
-                    if (n_buf > 0)
-                    {
-                        const auto& msk2 = (*bmasks)[pti];
-                        for (auto it = sep; it != pid.end(); ++it) {
-                            const long id = *it;
-                            const IntVect& iv = Index(aos[id], lev);
-                            inexflag[id] = msk2(iv);
-                        }
-
-                        auto sep2 = std::stable_partition(sep, pid.end(),
-                                                          [&inexflag](long id) {return inexflag[id];});
-                        if (bmasks == gather_masks) {
-                            nfine_gather = std::distance(pid.begin(), sep2);
-                        } else {
-                            nfine_current = std::distance(pid.begin(), sep2);
-                        }
-                    }
-                }
-
-                if (deposit_on_main_grid && lev > 0) {
-                    nfine_current = 0;
-                }
-
-                if (nfine_current != np || nfine_gather != np)
-                {
-                    particle_tmp.resize(np);
-                    for (long ip = 0; ip < np; ++ip) {
-                        particle_tmp[ip] = aos[pid[ip]];
-                    }
-                    std::swap(aos(), particle_tmp);
-
-                    tmp.resize(np);
-                    for (long ip = 0; ip < np; ++ip) {
-                        tmp[ip] = wp[pid[ip]];
-                    }
-                    std::swap(wp, tmp);
-
-                    for (long ip = 0; ip < np; ++ip) {
-                        tmp[ip] = uxp[pid[ip]];
-                    }
-                    std::swap(uxp, tmp);
-
-                    for (long ip = 0; ip < np; ++ip) {
-                        tmp[ip] = uyp[pid[ip]];
-                    }
-                    std::swap(uyp, tmp);
-
-                    for (long ip = 0; ip < np; ++ip) {
-                        tmp[ip] = uzp[pid[ip]];
-                    }
-                    std::swap(uzp, tmp);
-                }
-                BL_PROFILE_VAR_STOP(blp_partition);
+            if (has_buffer && !do_not_push) {
+                // - Modify `nfine_current` and `nfine_gather` (in place)
+                //    so that they correspond to the number of particles
+                //    that deposit/gather in the fine patch respectively.
+                // - Reorder the particle arrays,
+                //    so that the `nfine_current`/`nfine_gather` first particles
+                //    deposit/gather in the fine patch
+                //    and (thus) the `np-nfine_current`/`np-nfine_gather` last particles
+                //    deposit/gather in the buffer
+                PartitionParticlesInBuffers( nfine_current, nfine_gather, np,
+                    pti, lev, current_masks, gather_masks, uxp, uyp, uzp, wp );
             }
 
             const long np_current = (cjx) ? nfine_current : np;
-            
+
             //
             // copy data from particle container to temp arrays
             //
@@ -1195,14 +1078,14 @@ PhysicalParticleContainer::Evolve (int lev,
                 } else {
                     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 (has_buffer){
                     DepositCharge(pti, wp, ion_lev, crho, 0, np_current,
                                   np-np_current, thread_num, lev, lev-1);
                 }
             }
-            
+
             if (! do_not_push)
             {
                 const long np_gather = (cEx) ? nfine_gather : np;
@@ -1214,7 +1097,7 @@ PhysicalParticleContainer::Evolve (int lev,
                 //
                 BL_PROFILE_VAR_START(blp_fg);
                 FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
-                            exfab, eyfab, ezfab, bxfab, byfab, bzfab, 
+                            exfab, eyfab, ezfab, bxfab, byfab, bzfab,
                             Ex.nGrow(), e_is_nodal,
                             0, np_gather, thread_num, lev, lev);
 
@@ -1230,66 +1113,28 @@ PhysicalParticleContainer::Evolve (int lev,
                     FArrayBox const* cbxfab = &(*cBx)[pti];
                     FArrayBox const* cbyfab = &(*cBy)[pti];
                     FArrayBox const* cbzfab = &(*cBz)[pti];
-                    
+
                     if (WarpX::use_fdtd_nci_corr)
                     {
-#if (AMREX_SPACEDIM == 2)
-                        const Box& tbox = amrex::grow(cbox,{static_cast<int>(WarpX::nox),
-                                    static_cast<int>(WarpX::noz)});
-#else
-                        const Box& tbox = amrex::grow(cbox,{static_cast<int>(WarpX::nox),
-                                    static_cast<int>(WarpX::noy),
-                                    static_cast<int>(WarpX::noz)});
-#endif
-
-                        // Filter Ex (both 2D and 3D)
-                        filtered_Ex.resize(amrex::convert(tbox,WarpX::Ex_nodal_flag));
-                        // Safeguard for GPU
-                        exeli = filtered_Ex.elixir();
-                        // Apply filter on Ex, result stored in filtered_Ex
-                        nci_godfrey_filter_exeybz[lev-1]->ApplyStencil(filtered_Ex, (*cEx)[pti], filtered_Ex.box());
-                        // Update exfab reference
-                        cexfab = &filtered_Ex;
-
-                        // Filter Ez
-                        filtered_Ez.resize(amrex::convert(tbox,WarpX::Ez_nodal_flag));
-                        ezeli = filtered_Ez.elixir();
-                        nci_godfrey_filter_bxbyez[lev-1]->ApplyStencil(filtered_Ez, (*cEz)[pti], filtered_Ez.box());
-                        cezfab = &filtered_Ez;
-
-                        // Filter By
-                        filtered_By.resize(amrex::convert(tbox,WarpX::By_nodal_flag));
-                        byeli = filtered_By.elixir();
-                        nci_godfrey_filter_bxbyez[lev-1]->ApplyStencil(filtered_By, (*cBy)[pti], filtered_By.box());
-                        cbyfab = &filtered_By;
-#if (AMREX_SPACEDIM == 3)
-                        // Filter Ey
-                        filtered_Ey.resize(amrex::convert(tbox,WarpX::Ey_nodal_flag));
-                        eyeli = filtered_Ey.elixir();
-                        nci_godfrey_filter_exeybz[lev-1]->ApplyStencil(filtered_Ey, (*cEy)[pti], filtered_Ey.box());
-                        ceyfab = &filtered_Ey;
-                        
-                        // Filter Bx
-                        filtered_Bx.resize(amrex::convert(tbox,WarpX::Bx_nodal_flag));
-                        bxeli = filtered_Bx.elixir();
-                        nci_godfrey_filter_bxbyez[lev-1]->ApplyStencil(filtered_Bx, (*cBx)[pti], filtered_Bx.box());
-                        cbxfab = &filtered_Bx;
-                        
-                        // Filter Bz
-                        filtered_Bz.resize(amrex::convert(tbox,WarpX::Bz_nodal_flag));
-                        bzeli = filtered_Bz.elixir();
-                        nci_godfrey_filter_exeybz[lev-1]->ApplyStencil(filtered_Bz, (*cBz)[pti], filtered_Bz.box());
-                        cbzfab = &filtered_Bz;
-#endif
+                        // Filter arrays (*cEx)[pti], store the result in
+                        // filtered_Ex and update pointer cexfab so that it
+                        // points to filtered_Ex (and do the same for all
+                        // components of E and B)
+                        applyNCIFilter(lev-1, cbox, exeli, eyeli, ezeli, bxeli, byeli, bzeli,
+                                       filtered_Ex, filtered_Ey, filtered_Ez,
+                                       filtered_Bx, filtered_By, filtered_Bz,
+                                       (*cEx)[pti], (*cEy)[pti], (*cEz)[pti],
+                                       (*cBx)[pti], (*cBy)[pti], (*cBz)[pti],
+                                       cexfab, ceyfab, cezfab, cbxfab, cbyfab, cbzfab);
                     }
-                    
+
                     // Field gather for particles in gather buffers
                     e_is_nodal = cEx->is_nodal() and cEy->is_nodal() and cEz->is_nodal();
-                    FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp, 
+                    FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
                                 cexfab, ceyfab, cezfab,
                                 cbxfab, cbyfab, cbzfab,
-                                cEx->nGrow(), e_is_nodal, 
-                                nfine_gather, np-nfine_gather, 
+                                cEx->nGrow(), e_is_nodal,
+                                nfine_gather, np-nfine_gather,
                                 thread_num, lev, lev-1);
                 }
 
@@ -1299,8 +1144,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], 
-                       m_giv[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);
 
                 //
@@ -1313,7 +1157,7 @@ PhysicalParticleContainer::Evolve (int lev,
                 } else {
                     ion_lev = nullptr;
                 }
-                
+
                 // Deposit inside domains
                 DepositCurrent(pti, wp, uxp, uyp, uzp, ion_lev, &jx, &jy, &jz,
                                0, np_current, thread_num,
@@ -1333,7 +1177,7 @@ PhysicalParticleContainer::Evolve (int lev,
                 pti.SetPosition(m_xp[thread_num], m_yp[thread_num], m_zp[thread_num]);
                 BL_PROFILE_VAR_STOP(blp_copy);
             }
-            
+
             if (rho) {
                 // Deposit charge after particle push, in component 1 of MultiFab rho.
                 int* AMREX_RESTRICT ion_lev;
@@ -1362,8 +1206,84 @@ 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);
+    }
+}
+
+void
+PhysicalParticleContainer::applyNCIFilter (
+    int lev, const Box& box,
+    Elixir& exeli, Elixir& eyeli, Elixir& ezeli,
+    Elixir& bxeli, Elixir& byeli, Elixir& bzeli,
+    FArrayBox& filtered_Ex, FArrayBox& filtered_Ey, FArrayBox& filtered_Ez,
+    FArrayBox& filtered_Bx, FArrayBox& filtered_By, FArrayBox& filtered_Bz,
+    const FArrayBox& Ex, const FArrayBox& Ey, const FArrayBox& Ez,
+    const FArrayBox& Bx, const FArrayBox& By, const FArrayBox& Bz,
+    FArrayBox const * & ex_ptr, FArrayBox const * & ey_ptr,
+    FArrayBox const * & ez_ptr, FArrayBox const * & bx_ptr,
+    FArrayBox const * & by_ptr, FArrayBox const * & bz_ptr)
+{
+
+    // Get instances of NCI Godfrey filters
+    const auto& nci_godfrey_filter_exeybz = WarpX::GetInstance().nci_godfrey_filter_exeybz;
+    const auto& nci_godfrey_filter_bxbyez = WarpX::GetInstance().nci_godfrey_filter_bxbyez;
+
+#if (AMREX_SPACEDIM == 2)
+    const Box& tbox = amrex::grow(box,{static_cast<int>(WarpX::nox),
+                static_cast<int>(WarpX::noz)});
+#else
+    const Box& tbox = amrex::grow(box,{static_cast<int>(WarpX::nox),
+                static_cast<int>(WarpX::noy),
+                static_cast<int>(WarpX::noz)});
+#endif
+
+    // Filter Ex (Both 2D and 3D)
+    filtered_Ex.resize(amrex::convert(tbox,WarpX::Ex_nodal_flag));
+    // Safeguard for GPU
+    exeli = filtered_Ex.elixir();
+    // Apply filter on Ex, result stored in filtered_Ex
+    nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Ex, Ex, filtered_Ex.box());
+    // Update ex_ptr reference
+    ex_ptr = &filtered_Ex;
+
+    // Filter Ez
+    filtered_Ez.resize(amrex::convert(tbox,WarpX::Ez_nodal_flag));
+    ezeli = filtered_Ez.elixir();
+    nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_Ez, Ez, filtered_Ez.box());
+    ez_ptr = &filtered_Ez;
+
+    // Filter By
+    filtered_By.resize(amrex::convert(tbox,WarpX::By_nodal_flag));
+    byeli = filtered_By.elixir();
+    nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_By, By, filtered_By.box());
+    by_ptr = &filtered_By;
+#if (AMREX_SPACEDIM == 3)
+    // Filter Ey
+    filtered_Ey.resize(amrex::convert(tbox,WarpX::Ey_nodal_flag));
+    eyeli = filtered_Ey.elixir();
+    nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Ey, Ey, filtered_Ey.box());
+    ey_ptr = &filtered_Ey;
+
+    // Filter Bx
+    filtered_Bx.resize(amrex::convert(tbox,WarpX::Bx_nodal_flag));
+    bxeli = filtered_Bx.elixir();
+    nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_Bx, Bx, filtered_Bx.box());
+    bx_ptr = &filtered_Bx;
+
+    // Filter Bz
+    filtered_Bz.resize(amrex::convert(tbox,WarpX::Bz_nodal_flag));
+    bzeli = filtered_Bz.elixir();
+    nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Bz, Bz, filtered_Bz.box());
+    bz_ptr = &filtered_Bz;
+#endif
 }
 
 // Loop over all particles in the particle container and
@@ -1373,7 +1293,7 @@ PhysicalParticleContainer::SplitParticles(int lev)
 {
     auto& mypc = WarpX::GetInstance().GetPartContainer();
     auto& pctmp_split = mypc.GetPCtmp();
-    Cuda::ManagedDeviceVector<Real> xp, yp, zp;
+    Cuda::ManagedDeviceVector<ParticleReal> xp, yp, zp;
     RealVector psplit_x, psplit_y, psplit_z, psplit_w;
     RealVector psplit_ux, psplit_uy, psplit_uz;
     long np_split_to_add = 0;
@@ -1389,7 +1309,16 @@ PhysicalParticleContainer::SplitParticles(int lev)
     for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti)
     {
         pti.GetPosition(xp, yp, zp);
+
+        // offset for split particles is computed as a function of cell size
+        // and number of particles per cell, so that a uniform distribution
+        // before splitting results in a uniform distribution after splitting
+        const amrex::Vector<int> ppc_nd = plasma_injector->num_particles_per_cell_each_dim;
         const std::array<Real,3>& dx = WarpX::CellSize(lev);
+        amrex::Vector<amrex::Real> split_offset = {dx[0]/2./ppc_nd[0],
+                                                   dx[1]/2./ppc_nd[1],
+                                                   dx[2]/2./ppc_nd[2]};
+
         // particle Array Of Structs data
         auto& particles = pti.GetArrayOfStructs();
         // particle Struct Of Arrays data
@@ -1402,19 +1331,19 @@ PhysicalParticleContainer::SplitParticles(int lev)
         for(int i=0; i<np; i++){
             auto& p = particles[i];
             if (p.id() == DoSplitParticleID){
-                // If particle is tagged, split it and put the 
+                // If particle is tagged, split it and put the
                 // split particles in local arrays psplit_x etc.
                 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 ){
                             // Add one particle with offset in x and z
-                            psplit_x.push_back( xp[i] + ishift*dx[0]/2 );
+                            psplit_x.push_back( xp[i] + ishift*split_offset[0] );
                             psplit_y.push_back( yp[i] );
-                            psplit_z.push_back( zp[i] + kshift*dx[2]/2 );
+                            psplit_z.push_back( zp[i] + kshift*split_offset[2] );
                             psplit_ux.push_back( uxp[i] );
                             psplit_uy.push_back( uyp[i] );
                             psplit_uz.push_back( uzp[i] );
@@ -1422,11 +1351,11 @@ 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
-                        psplit_x.push_back( xp[i] + ishift*dx[0]/2 );
+                        psplit_x.push_back( xp[i] + ishift*split_offset[0] );
                         psplit_y.push_back( yp[i] );
                         psplit_z.push_back( zp[i] );
                         psplit_ux.push_back( uxp[i] );
@@ -1436,7 +1365,7 @@ PhysicalParticleContainer::SplitParticles(int lev)
                         // Add one particle with offset in z
                         psplit_x.push_back( xp[i] );
                         psplit_y.push_back( yp[i] );
-                        psplit_z.push_back( zp[i] + ishift*dx[2]/2 );
+                        psplit_z.push_back( zp[i] + ishift*split_offset[2] );
                         psplit_ux.push_back( uxp[i] );
                         psplit_uy.push_back( uyp[i] );
                         psplit_uz.push_back( uzp[i] );
@@ -1445,15 +1374,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_x.push_back( xp[i] + ishift*split_offset[0] );
+                                psplit_y.push_back( yp[i] + jshift*split_offset[1] );
+                                psplit_z.push_back( zp[i] + kshift*split_offset[2] );
                                 psplit_ux.push_back( uxp[i] );
                                 psplit_uy.push_back( uyp[i] );
                                 psplit_uz.push_back( uzp[i] );
@@ -1462,11 +1391,11 @@ 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 );
+                        psplit_x.push_back( xp[i] + ishift*split_offset[0] );
                         psplit_y.push_back( yp[i] );
                         psplit_z.push_back( zp[i] );
                         psplit_ux.push_back( uxp[i] );
@@ -1475,7 +1404,7 @@ PhysicalParticleContainer::SplitParticles(int lev)
                         psplit_w.push_back( wp[i]/np_split );
                         // Add one particle with offset in y
                         psplit_x.push_back( xp[i] );
-                        psplit_y.push_back( yp[i] + ishift*dx[1]/2 );
+                        psplit_y.push_back( yp[i] + ishift*split_offset[1] );
                         psplit_z.push_back( zp[i] );
                         psplit_ux.push_back( uxp[i] );
                         psplit_uy.push_back( uyp[i] );
@@ -1484,7 +1413,7 @@ PhysicalParticleContainer::SplitParticles(int lev)
                         // Add one particle with offset in z
                         psplit_x.push_back( xp[i] );
                         psplit_y.push_back( yp[i] );
-                        psplit_z.push_back( zp[i] + ishift*dx[2]/2 );
+                        psplit_z.push_back( zp[i] + ishift*split_offset[2] );
                         psplit_ux.push_back( uxp[i] );
                         psplit_uy.push_back( uyp[i] );
                         psplit_uz.push_back( uzp[i] );
@@ -1497,13 +1426,13 @@ PhysicalParticleContainer::SplitParticles(int lev)
             }
         }
     }
-	// Add local arrays psplit_x etc. to the temporary
-	// particle container pctmp_split. Split particles
-	// are tagged with p.id()=NoSplitParticleID so that 
-	// they are not re-split when entering a higher level
-	// AddNParticles calls Redistribute, so that particles
-	// in pctmp_split are in the proper grids and tiles
-	pctmp_split.AddNParticles(lev, 
+    // Add local arrays psplit_x etc. to the temporary
+    // particle container pctmp_split. Split particles
+    // are tagged with p.id()=NoSplitParticleID so that
+    // they are not re-split when entering a higher level
+    // AddNParticles calls Redistribute, so that particles
+    // in pctmp_split are in the proper grids and tiles
+    pctmp_split.AddNParticles(lev,
                               np_split_to_add,
                               psplit_x.dataPtr(),
                               psplit_y.dataPtr(),
@@ -1514,39 +1443,37 @@ PhysicalParticleContainer::SplitParticles(int lev)
                               1,
                               psplit_w.dataPtr(),
                               1, NoSplitParticleID);
-	// Copy particles from tmp to current particle container
+    // Copy particles from tmp to current particle container
     addParticles(pctmp_split,1);
-	// Clear tmp container
+    // Clear tmp container
     pctmp_split.clearParticles();
 }
 
 void
 PhysicalParticleContainer::PushPX(WarpXParIter& pti,
-                                  Cuda::ManagedDeviceVector<Real>& xp,
-                                  Cuda::ManagedDeviceVector<Real>& yp,
-                                  Cuda::ManagedDeviceVector<Real>& zp,
-                                  Cuda::ManagedDeviceVector<Real>& giv,
-                                  Real dt)
+                                  Cuda::ManagedDeviceVector<ParticleReal>& xp,
+                                  Cuda::ManagedDeviceVector<ParticleReal>& yp,
+                                  Cuda::ManagedDeviceVector<ParticleReal>& zp,
+                                  Real dt, DtType a_dt_type)
 {
 
     // This wraps the momentum and position advance so that inheritors can modify the call.
     auto& attribs = pti.GetAttribs();
     // Extract pointers to the different particle quantities
-    Real* const AMREX_RESTRICT x = xp.dataPtr();
-    Real* const AMREX_RESTRICT y = yp.dataPtr();
-    Real* const AMREX_RESTRICT z = zp.dataPtr();
-    Real* const AMREX_RESTRICT gi = giv.dataPtr();
-    Real* const AMREX_RESTRICT ux = attribs[PIdx::ux].dataPtr();
-    Real* const AMREX_RESTRICT uy = attribs[PIdx::uy].dataPtr();
-    Real* const AMREX_RESTRICT uz = attribs[PIdx::uz].dataPtr();
-    const Real* const AMREX_RESTRICT Ex = attribs[PIdx::Ex].dataPtr();
-    const Real* const AMREX_RESTRICT Ey = attribs[PIdx::Ey].dataPtr();
-    const Real* const AMREX_RESTRICT Ez = attribs[PIdx::Ez].dataPtr();
-    const Real* const AMREX_RESTRICT Bx = attribs[PIdx::Bx].dataPtr();
-    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)
+    ParticleReal* const AMREX_RESTRICT x = xp.dataPtr();
+    ParticleReal* const AMREX_RESTRICT y = yp.dataPtr();
+    ParticleReal* const AMREX_RESTRICT z = zp.dataPtr();
+    ParticleReal* const AMREX_RESTRICT ux = attribs[PIdx::ux].dataPtr();
+    ParticleReal* const AMREX_RESTRICT uy = attribs[PIdx::uy].dataPtr();
+    ParticleReal* const AMREX_RESTRICT uz = attribs[PIdx::uz].dataPtr();
+    const ParticleReal* const AMREX_RESTRICT Ex = attribs[PIdx::Ex].dataPtr();
+    const ParticleReal* const AMREX_RESTRICT Ey = attribs[PIdx::Ey].dataPtr();
+    const ParticleReal* const AMREX_RESTRICT Ez = attribs[PIdx::Ez].dataPtr();
+    const ParticleReal* const AMREX_RESTRICT Bx = attribs[PIdx::Bx].dataPtr();
+    const ParticleReal* const AMREX_RESTRICT By = attribs[PIdx::By].dataPtr();
+    const ParticleReal* const AMREX_RESTRICT Bz = attribs[PIdx::Bz].dataPtr();
+
+    if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags && (a_dt_type!=DtType::SecondHalf))
     {
         copy_attribs(pti, x, y, z);
     }
@@ -1555,17 +1482,17 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti,
     if (do_field_ionization){
         ion_lev = pti.GetiAttribs(particle_icomps["ionization_level"]).dataPtr();
     }
-    
+
     // 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( 
+        amrex::ParallelFor(
             pti.numParticles(),
             [=] AMREX_GPU_DEVICE (long i) {
                 Real qp = q;
                 if (ion_lev){ qp *= ion_lev[i]; }
-                UpdateMomentumBoris( ux[i], uy[i], uz[i], gi[i],
+                UpdateMomentumBoris( 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],
@@ -1578,7 +1505,7 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti,
             [=] AMREX_GPU_DEVICE (long i) {
                 Real qp = q;
                 if (ion_lev){ qp *= ion_lev[i]; }
-                UpdateMomentumVay( ux[i], uy[i], uz[i], gi[i],
+                UpdateMomentumVay( 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],
@@ -1609,7 +1536,7 @@ PhysicalParticleContainer::PushP (int lev, Real dt,
         int thread_num = omp_get_thread_num();
 #else
         int thread_num = 0;
-#endif      
+#endif
         for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti)
         {
             const Box& box = pti.validbox();
@@ -1640,8 +1567,6 @@ PhysicalParticleContainer::PushP (int lev, Real dt,
             Byp.assign(np,WarpX::B_external[1]);
             Bzp.assign(np,WarpX::B_external[2]);
 
-            m_giv[thread_num].resize(np);
-
             //
             // copy data from particle container to temp arrays
             //
@@ -1649,22 +1574,21 @@ PhysicalParticleContainer::PushP (int lev, Real dt,
 
             int e_is_nodal = Ex.is_nodal() and Ey.is_nodal() and Ez.is_nodal();
             FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
-                        &exfab, &eyfab, &ezfab, &bxfab, &byfab, &bzfab, 
+                        &exfab, &eyfab, &ezfab, &bxfab, &byfab, &bzfab,
                         Ex.nGrow(), e_is_nodal,
                         0, np, thread_num, lev, lev);
 
             // This wraps the momentum advance so that inheritors can modify the call.
             // Extract pointers to the different particle quantities
-            Real* const AMREX_RESTRICT gi = m_giv[thread_num].dataPtr();
-            Real* const AMREX_RESTRICT ux = attribs[PIdx::ux].dataPtr();
-            Real* const AMREX_RESTRICT uy = attribs[PIdx::uy].dataPtr();
-            Real* const AMREX_RESTRICT uz = attribs[PIdx::uz].dataPtr();
-            const Real* const AMREX_RESTRICT Expp = Exp.dataPtr();
-            const Real* const AMREX_RESTRICT Eypp = Eyp.dataPtr();
-            const Real* const AMREX_RESTRICT Ezpp = Ezp.dataPtr();
-            const Real* const AMREX_RESTRICT Bxpp = Bxp.dataPtr();
-            const Real* const AMREX_RESTRICT Bypp = Byp.dataPtr();
-            const Real* const AMREX_RESTRICT Bzpp = Bzp.dataPtr();
+            ParticleReal* const AMREX_RESTRICT ux = attribs[PIdx::ux].dataPtr();
+            ParticleReal* const AMREX_RESTRICT uy = attribs[PIdx::uy].dataPtr();
+            ParticleReal* const AMREX_RESTRICT uz = attribs[PIdx::uz].dataPtr();
+            const ParticleReal* const AMREX_RESTRICT Expp = Exp.dataPtr();
+            const ParticleReal* const AMREX_RESTRICT Eypp = Eyp.dataPtr();
+            const ParticleReal* const AMREX_RESTRICT Ezpp = Ezp.dataPtr();
+            const ParticleReal* const AMREX_RESTRICT Bxpp = Bxp.dataPtr();
+            const ParticleReal* const AMREX_RESTRICT Bypp = Byp.dataPtr();
+            const ParticleReal* const AMREX_RESTRICT Bzpp = Bzp.dataPtr();
 
             // Loop over the particles and update their momentum
             const Real q = this->charge;
@@ -1672,14 +1596,14 @@ PhysicalParticleContainer::PushP (int lev, Real dt,
             if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Boris){
                 amrex::ParallelFor( pti.numParticles(),
                     [=] AMREX_GPU_DEVICE (long i) {
-                        UpdateMomentumBoris( ux[i], uy[i], uz[i], gi[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) {
                 amrex::ParallelFor( pti.numParticles(),
                     [=] AMREX_GPU_DEVICE (long i) {
-                        UpdateMomentumVay( ux[i], uy[i], uz[i], gi[i],
+                        UpdateMomentumVay( ux[i], uy[i], uz[i],
                               Expp[i], Eypp[i], Ezpp[i], Bxpp[i], Bypp[i], Bzpp[i], q, m, dt);
                     }
                 );
@@ -1690,30 +1614,30 @@ PhysicalParticleContainer::PushP (int lev, Real dt,
     }
 }
 
-void PhysicalParticleContainer::copy_attribs(WarpXParIter& pti,const Real* xp,
-                                             const Real* yp, const Real* zp)
+void PhysicalParticleContainer::copy_attribs(WarpXParIter& pti,const ParticleReal* xp,
+                                             const ParticleReal* yp, const ParticleReal* zp)
 {
     auto& attribs = pti.GetAttribs();
-    Real* AMREX_RESTRICT uxp = attribs[PIdx::ux].dataPtr();
-    Real* AMREX_RESTRICT uyp = attribs[PIdx::uy].dataPtr();
-    Real* AMREX_RESTRICT uzp = attribs[PIdx::uz].dataPtr();
-    
+    ParticleReal* AMREX_RESTRICT uxp = attribs[PIdx::ux].dataPtr();
+    ParticleReal* AMREX_RESTRICT uyp = attribs[PIdx::uy].dataPtr();
+    ParticleReal* AMREX_RESTRICT uzp = attribs[PIdx::uz].dataPtr();
+
     const auto np = pti.numParticles();
     const auto lev = pti.GetLevel();
     const auto index = pti.GetPairIndex();
-    Real* AMREX_RESTRICT xpold  = tmp_particle_data[lev][index][TmpIdx::xold ].dataPtr();
-    Real* AMREX_RESTRICT ypold  = tmp_particle_data[lev][index][TmpIdx::yold ].dataPtr();
-    Real* AMREX_RESTRICT zpold  = tmp_particle_data[lev][index][TmpIdx::zold ].dataPtr();
-    Real* AMREX_RESTRICT uxpold = tmp_particle_data[lev][index][TmpIdx::uxold].dataPtr();
-    Real* AMREX_RESTRICT uypold = tmp_particle_data[lev][index][TmpIdx::uyold].dataPtr();
-    Real* AMREX_RESTRICT uzpold = tmp_particle_data[lev][index][TmpIdx::uzold].dataPtr();
+    ParticleReal* AMREX_RESTRICT xpold  = tmp_particle_data[lev][index][TmpIdx::xold ].dataPtr();
+    ParticleReal* AMREX_RESTRICT ypold  = tmp_particle_data[lev][index][TmpIdx::yold ].dataPtr();
+    ParticleReal* AMREX_RESTRICT zpold  = tmp_particle_data[lev][index][TmpIdx::zold ].dataPtr();
+    ParticleReal* AMREX_RESTRICT uxpold = tmp_particle_data[lev][index][TmpIdx::uxold].dataPtr();
+    ParticleReal* AMREX_RESTRICT uypold = tmp_particle_data[lev][index][TmpIdx::uyold].dataPtr();
+    ParticleReal* AMREX_RESTRICT uzpold = tmp_particle_data[lev][index][TmpIdx::uzold].dataPtr();
 
     ParallelFor( np,
                  [=] AMREX_GPU_DEVICE (long i) {
                      xpold[i]=xp[i];
                      ypold[i]=yp[i];
                      zpold[i]=zp[i];
-            
+
                      uxpold[i]=uxp[i];
                      uypold[i]=uyp[i];
                      uzpold[i]=uzp[i];
@@ -1753,7 +1677,7 @@ void PhysicalParticleContainer::GetParticleSlice(const int direction, const Real
     slice_box.setHi(direction, z_max);
 
     diagnostic_particles.resize(finestLevel()+1);
-    
+
     for (int lev = 0; lev < nlevs; ++lev) {
 
         const Real* dx  = Geom(lev).CellSize();
@@ -1800,7 +1724,7 @@ void PhysicalParticleContainer::GetParticleSlice(const int direction, const Real
                 auto& uzp_old = tmp_particle_data[lev][index][TmpIdx::uzold];
 
                 const long np = pti.numParticles();
-                
+
                 Real uzfrm = -WarpX::gamma_boost*WarpX::beta_boost*PhysConst::c;
                 Real inv_c2 = 1.0/PhysConst::c/PhysConst::c;
 
@@ -1835,7 +1759,7 @@ void PhysicalParticleContainer::GetParticleSlice(const int direction, const Real
                     Real uzp = uz_old_p  *weight_old + uz_new_p  *weight_new;
 
                     diagnostic_particles[lev][index].GetRealData(DiagIdx::w).push_back(wp[i]);
-                    
+
                     diagnostic_particles[lev][index].GetRealData(DiagIdx::x).push_back(xp);
                     diagnostic_particles[lev][index].GetRealData(DiagIdx::y).push_back(yp);
                     diagnostic_particles[lev][index].GetRealData(DiagIdx::z).push_back(zp);
@@ -1860,7 +1784,7 @@ PhysicalParticleContainer::ContinuousInjection(const RealBox& injection_box)
     AddPlasma(lev, injection_box);
 }
 
-/* \brief Gather fields from FArrayBox exfab, eyfab, ezfab, bxfab, byfab, 
+/* \brief Gather fields from FArrayBox exfab, eyfab, ezfab, bxfab, byfab,
  * bzfab into arrays of fields on particles Exp, Eyp, Ezp, Bxp, Byp, Bzp.
  * \param Exp-Bzp: fields on particles.
  * \param exfab-bzfab: FAB of electric and magnetic fields for particles in pti
@@ -1870,7 +1794,7 @@ PhysicalParticleContainer::ContinuousInjection(const RealBox& injection_box)
  * \param np_to_gather: number of particles onto which fields are gathered
  * \param thread_num: if using OpenMP, thread number
  * \param lev: level on which particles are located
- * \param gather_lev: level from which particles gather fields (lev-1) for 
+ * \param gather_lev: level from which particles gather fields (lev-1) for
           particles in buffers.
  */
 void
@@ -1897,14 +1821,14 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti,
     AMREX_ALWAYS_ASSERT_WITH_MESSAGE((gather_lev==(lev-1)) ||
                                      (gather_lev==(lev  )),
                                      "Gather buffers only work for lev-1");
-    
+
     // If no particles, do not do anything
     if (np_to_gather == 0) return;
     // Get cell size on gather_lev
     const std::array<Real,3>& dx = WarpX::CellSize(std::max(gather_lev,0));
     // Set staggering shift depending on e_is_nodal
     const Real stagger_shift = e_is_nodal ? 0.0 : 0.5;
-    
+
     // Get box from which field is gathered.
     // If not gathering from the finest level, the box is coarsened.
     Box box;
@@ -1914,26 +1838,26 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti,
         const IntVect& ref_ratio = WarpX::RefRatio(gather_lev);
         box = amrex::coarsen(pti.tilebox(),ref_ratio);
     }
-    
+
     // Add guard cells to the box.
     box.grow(ngE);
-    
+
     const Array4<const Real>& ex_arr = exfab->array();
     const Array4<const Real>& ey_arr = eyfab->array();
     const Array4<const Real>& ez_arr = ezfab->array();
     const Array4<const Real>& bx_arr = bxfab->array();
     const Array4<const Real>& by_arr = byfab->array();
     const Array4<const Real>& bz_arr = bzfab->array();
-    
-    const Real * const AMREX_RESTRICT xp = m_xp[thread_num].dataPtr() + offset;
-    const Real * const AMREX_RESTRICT zp = m_zp[thread_num].dataPtr() + offset;
-    const Real * const AMREX_RESTRICT yp = m_yp[thread_num].dataPtr() + offset;
-    
+
+    const ParticleReal * const AMREX_RESTRICT xp = m_xp[thread_num].dataPtr() + offset;
+    const ParticleReal * const AMREX_RESTRICT zp = m_zp[thread_num].dataPtr() + offset;
+    const ParticleReal * const AMREX_RESTRICT yp = m_yp[thread_num].dataPtr() + offset;
+
     // Lower corner of tile box physical domain
     const std::array<Real, 3>& xyzmin = WarpX::LowerCorner(box, gather_lev);
-    
+
     const Dim3 lo = lbound(box);
-    
+
     // Depending on l_lower_in_v and WarpX::nox, call
     // different versions of template function doGatherShapeN
     if (WarpX::l_lower_order_in_v){
@@ -2016,10 +1940,10 @@ void PhysicalParticleContainer::InitIonizationModule ()
     }
     // Compute ADK prefactors (See Chen, JCP 236 (2013), equation (2))
     // For now, we assume l=0 and m=0.
-    // The approximate expressions are used, 
+    // The approximate expressions are used,
     // without Gamma function
     Real wa = std::pow(PhysConst::alpha,3) * PhysConst::c / PhysConst::r_e;
-    Real Ea = PhysConst::m_e * PhysConst::c*PhysConst::c /PhysConst::q_e * 
+    Real Ea = PhysConst::m_e * PhysConst::c*PhysConst::c /PhysConst::q_e *
         std::pow(PhysConst::alpha,4)/PhysConst::r_e;
     Real UH = table_ionization_energies[0];
     Real l_eff = std::sqrt(UH/ionization_energies[0]) - 1.;
@@ -2034,18 +1958,18 @@ void PhysicalParticleContainer::InitIonizationModule ()
         Real C2 = std::pow(2,2*n_eff)/(n_eff*tgamma(n_eff+l_eff+1)*tgamma(n_eff-l_eff));
         adk_power[i] = -(2*n_eff - 1);
         Real Uion = ionization_energies[i];
-        adk_prefactor[i] = dt * wa * C2 * ( Uion/(2*UH) ) 
+        adk_prefactor[i] = dt * wa * C2 * ( Uion/(2*UH) )
             * std::pow(2*std::pow((Uion/UH),3./2)*Ea,2*n_eff - 1);
         adk_exp_prefactor[i] = -2./3 * std::pow( Uion/UH,3./2) * Ea;
     }
 }
 
 /* \brief create mask of ionized particles (1 if ionized, 0 otherwise)
- * 
+ *
  * \param mfi: tile or grid
  * \param lev: MR level
  * \param ionization_mask: Array with as many elements as particles in mfi.
- * This function initialized the array, and set each element to 1 or 0 
+ * This function initialized the array, and set each element to 1 or 0
  * depending on whether the particle is ionized or not.
  */
 void
@@ -2074,15 +1998,15 @@ PhysicalParticleContainer::buildIonizationMask (const amrex::MFIter& mfi, const
     // Otherwise, resize ionization_mask, and get poiters to attribs arrays.
     ionization_mask.resize(np);
     int * const AMREX_RESTRICT p_ionization_mask = ionization_mask.data();
-    const Real * const AMREX_RESTRICT ux = soa.GetRealData(PIdx::ux).data();
-    const Real * const AMREX_RESTRICT uy = soa.GetRealData(PIdx::uy).data();
-    const Real * const AMREX_RESTRICT uz = soa.GetRealData(PIdx::uz).data();
-    const Real * const AMREX_RESTRICT ex = soa.GetRealData(PIdx::Ex).data();
-    const Real * const AMREX_RESTRICT ey = soa.GetRealData(PIdx::Ey).data();
-    const Real * const AMREX_RESTRICT ez = soa.GetRealData(PIdx::Ez).data();
-    const Real * const AMREX_RESTRICT bx = soa.GetRealData(PIdx::Bx).data();
-    const Real * const AMREX_RESTRICT by = soa.GetRealData(PIdx::By).data();
-    const Real * const AMREX_RESTRICT bz = soa.GetRealData(PIdx::Bz).data();
+    const ParticleReal * const AMREX_RESTRICT ux = soa.GetRealData(PIdx::ux).data();
+    const ParticleReal * const AMREX_RESTRICT uy = soa.GetRealData(PIdx::uy).data();
+    const ParticleReal * const AMREX_RESTRICT uz = soa.GetRealData(PIdx::uz).data();
+    const ParticleReal * const AMREX_RESTRICT ex = soa.GetRealData(PIdx::Ex).data();
+    const ParticleReal * const AMREX_RESTRICT ey = soa.GetRealData(PIdx::Ey).data();
+    const ParticleReal * const AMREX_RESTRICT ez = soa.GetRealData(PIdx::Ez).data();
+    const ParticleReal * const AMREX_RESTRICT bx = soa.GetRealData(PIdx::Bx).data();
+    const ParticleReal * const AMREX_RESTRICT by = soa.GetRealData(PIdx::By).data();
+    const ParticleReal * const AMREX_RESTRICT bz = soa.GetRealData(PIdx::Bz).data();
     int* ion_lev = soa.GetIntData(particle_icomps["ionization_level"]).data();
 
     Real c = PhysConst::c;
@@ -2091,7 +2015,7 @@ PhysicalParticleContainer::buildIonizationMask (const amrex::MFIter& mfi, const
 
     // Loop over all particles in grid/tile. If ionized, set mask to 1
     // and increment ionization level.
-    ParallelFor( 
+    ParallelFor(
         np,
         [=] AMREX_GPU_DEVICE (long i) {
             // Get index of ionization_level