Fix uninitialized particle data (#2282)

1 files changed, 85 insertions, 6 deletions

diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp
index c704e824a..5fd1b9ced 100644
--- a/Source/Particles/PhysicalParticleContainer.cpp
+++ b/Source/Particles/PhysicalParticleContainer.cpp
@@ -107,6 +107,8 @@ using namespace amrex;
 
 namespace
 {
+    using ParticleType = WarpXParticleContainer::ParticleType;
+
     // Since the user provides the density distribution
     // at t_lab=0 and in the lab-frame coordinates,
     // we need to find the lab-frame position of this
@@ -155,6 +157,53 @@ namespace
 #endif
         return pos;
     }
+
+    /**
+     * \brief This function is called in AddPlasma when we want a particle to be removed at the
+     * next call to redistribute. It initializes all the particle properties to zero (to be safe
+     * and avoid any possible undefined behavior before the next call to redistribute) and sets
+     * the particle id to -1 so that it can be effectively deleted.
+     *
+     * \param p particle aos data
+     * \param pa particle soa data
+     * \param ip index for soa data
+     * \param do_field_ionization whether species has ionization
+     * \param pi ionization level data
+     * \param has_quantum_sync whether species has quantum synchrotron
+     * \param p_optical_depth_QSR quantum synchrotron optical depth data
+     * \param has_breit_wheeler whether species has Breit-Wheeler
+     * \param p_optical_depth_BW Breit-Wheeler optical depth data
+     */
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    void ZeroInitializeAndSetNegativeID (
+        ParticleType& p, const GpuArray<ParticleReal*,PIdx::nattribs>& pa, long& ip,
+        const bool& do_field_ionization, int* pi
+#ifdef WARPX_QED
+        ,const bool& has_quantum_sync, amrex::ParticleReal* p_optical_depth_QSR
+        ,const bool& has_breit_wheeler, amrex::ParticleReal* p_optical_depth_BW
+#endif
+        ) noexcept
+    {
+        p.pos(0) = 0._rt;
+        p.pos(1) = 0._rt;
+#if (AMREX_SPACEDIM == 3)
+        p.pos(2) = 0._rt;
+#endif
+        pa[PIdx::w ][ip] = 0._rt;
+        pa[PIdx::ux][ip] = 0._rt;
+        pa[PIdx::uy][ip] = 0._rt;
+        pa[PIdx::uz][ip] = 0._rt;
+#ifdef WARPX_DIM_RZ
+        pa[PIdx::theta][ip] = 0._rt;
+#endif
+        if (do_field_ionization) {pi[ip] = 0;}
+#ifdef WARPX_QED
+        if (has_quantum_sync) {p_optical_depth_QSR[ip] = 0._rt;}
+        if (has_breit_wheeler) {p_optical_depth_BW[ip] = 0._rt;}
+#endif
+
+        p.id() = -1;
+    }
 }
 
 PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int ispecies,
@@ -883,13 +932,23 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox)
 
 #if (AMREX_SPACEDIM == 3)
                 if (!tile_realbox.contains(XDim3{pos.x,pos.y,pos.z})) {
-                    p.id() = -1;
+                    ZeroInitializeAndSetNegativeID(p, pa, ip, loc_do_field_ionization, pi
+#ifdef WARPX_QED
+                                                   ,loc_has_quantum_sync, p_optical_depth_QSR
+                                                   ,loc_has_breit_wheeler, p_optical_depth_BW
+#endif
+                                                   );
                     continue;
                 }
 #else
                 amrex::ignore_unused(k);
                 if (!tile_realbox.contains(XDim3{pos.x,pos.z,0.0_rt})) {
-                    p.id() = -1;
+                    ZeroInitializeAndSetNegativeID(p, pa, ip, loc_do_field_ionization, pi
+#ifdef WARPX_QED
+                                                   ,loc_has_quantum_sync, p_optical_depth_QSR
+                                                   ,loc_has_breit_wheeler, p_optical_depth_BW
+#endif
+                                                   );
                     continue;
                 }
 #endif
@@ -926,7 +985,12 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox)
                     const Real z0 = applyBallisticCorrection(pos, inj_mom, gamma_boost,
                                                              beta_boost, t);
                     if (!inj_pos->insideBounds(xb, yb, z0)) {
-                        p.id() = -1;
+                        ZeroInitializeAndSetNegativeID(p, pa, ip, loc_do_field_ionization, pi
+#ifdef WARPX_QED
+                                                   ,loc_has_quantum_sync, p_optical_depth_QSR
+                                                   ,loc_has_breit_wheeler, p_optical_depth_BW
+#endif
+                                                   );
                         continue;
                     }
 
@@ -935,7 +999,12 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox)
 
                     // Remove particle if density below threshold
                     if ( dens < density_min ){
-                        p.id() = -1;
+                        ZeroInitializeAndSetNegativeID(p, pa, ip, loc_do_field_ionization, pi
+#ifdef WARPX_QED
+                                                   ,loc_has_quantum_sync, p_optical_depth_QSR
+                                                   ,loc_has_breit_wheeler, p_optical_depth_BW
+#endif
+                                                   );
                         continue;
                     }
                     // Cut density if above threshold
@@ -948,14 +1017,24 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox)
                     // If the particle is not within the lab-frame zmin, zmax, etc.
                     // go to the next generated particle.
                     if (!inj_pos->insideBounds(xb, yb, z0_lab)) {
-                        p.id() = -1;
+                        ZeroInitializeAndSetNegativeID(p, pa, ip, loc_do_field_ionization, pi
+#ifdef WARPX_QED
+                                                   ,loc_has_quantum_sync, p_optical_depth_QSR
+                                                   ,loc_has_breit_wheeler, p_optical_depth_BW
+#endif
+                                                   );
                         continue;
                     }
                     // call `getDensity` with lab-frame parameters
                     dens = inj_rho->getDensity(pos.x, pos.y, z0_lab);
                     // Remove particle if density below threshold
                     if ( dens < density_min ){
-                        p.id() = -1;
+                        ZeroInitializeAndSetNegativeID(p, pa, ip, loc_do_field_ionization, pi
+#ifdef WARPX_QED
+                                                   ,loc_has_quantum_sync, p_optical_depth_QSR
+                                                   ,loc_has_breit_wheeler, p_optical_depth_BW
+#endif
+                                                   );
                         continue;
                     }
                     // Cut density if above threshold