Draft: Copy Lost Particles

to do:
- remove in beam species
- output

src/ImpactX.H

@@ -134,6 +134,9 @@ namespace impactx
         /** these are the physical/beam particles of the simulation */
         std::unique_ptr<ImpactXParticleContainer> m_particle_container;
 
+        /** former beam particles that got lost in apertures, the wall, etc. */
+        std::unique_ptr<ImpactXParticleContainer> m_particles_lost;
+
         /** charge density per level */
         std::unordered_map<int, amrex::MultiFab> m_rho;
         /** scalar potential per level */

src/ImpactX.cpp

@@ -9,6 +9,7 @@
  */
 #include "ImpactX.H"
 #include "initialization/InitAmrCore.H"
+#include "particles/CollectLost.H"
 #include "particles/ImpactXParticleContainer.H"
 #include "particles/Push.H"
 #include "particles/diagnostics/DiagnosticOutput.H"
@@ -186,6 +187,9 @@ namespace impactx
                     // push all particles with external maps
                     Push(*m_particle_container, element_variant, global_step);
 
+                    // move "lost" particles to another particle container
+                    collect_lost_particles(*m_particle_container, *m_particles_lost);
+
                     // just prints an empty newline at the end of the slice_step
                     amrex::Print() << "\n";
 

src/particles/CMakeLists.txt

@@ -1,6 +1,7 @@
 target_sources(ImpactX
   PRIVATE
     ChargeDeposition.cpp
+    CollectLost.cpp
     ImpactXParticleContainer.cpp
     Push.cpp
 )

src/particles/CollectLost.H

@@ -0,0 +1,29 @@
+/* Copyright 2022-2023 The Regents of the University of California, through Lawrence
+ *           Berkeley National Laboratory (subject to receipt of any required
+ *           approvals from the U.S. Dept. of Energy). All rights reserved.
+ *
+ * This file is part of ImpactX.
+ *
+ * Authors: Axel Huebl, Chad Mitchell
+ * License: BSD-3-Clause-LBNL
+ */
+#ifndef IMPACTX_COLLECT_LOST_H
+#define IMPACTX_COLLECT_LOST_H
+
+#include "particles/ImpactXParticleContainer.H"
+
+
+namespace impactx
+{
+    /** ...
+     *
+     * ...
+     *
+     * @param source
+     * @param dest
+     */
+    void collect_lost_particles (ImpactXParticleContainer& source, ImpactXParticleContainer& dest);
+
+} // namespace impactx
+
+#endif // IMPACTX_COLLECT_LOST_H

src/particles/CollectLost.cpp

@@ -0,0 +1,58 @@
+/* Copyright 2022-2023 The Regents of the University of California, through Lawrence
+ *           Berkeley National Laboratory (subject to receipt of any required
+ *           approvals from the U.S. Dept. of Energy). All rights reserved.
+ *
+ * This file is part of ImpactX.
+ *
+ * Authors: Axel Huebl, Chad Mitchell
+ * License: BSD-3-Clause-LBNL
+ */
+#include "CollectLost.H"
+
+#include <AMReX_GpuLaunch.H>
+#include <AMReX_GpuQualifiers.H>
+
+
+namespace impactx
+{
+    void collect_lost_particles (ImpactXParticleContainer& source, ImpactXParticleContainer& dest)
+    {
+        using SrcData = ImpactXParticleContainer::ParticleTileType::ConstParticleTileDataType;
+
+        // copy all particles marked with a negative ID from source to destination
+        bool const local = true;
+        dest.copyParticles(
+                source,
+                [=] AMREX_GPU_HOST_DEVICE(const SrcData &src, int ip) {
+                    return src.id(ip) < 0;
+                },
+                local
+        );
+
+        // flip IDs back to positive in destination
+        int const nLevel = dest.finestLevel();
+        for (int lev = 0; lev <= nLevel; ++lev) {
+            // loop over all particle boxes
+            using ParIt = ImpactXParticleContainer::iterator;
+
+#ifdef AMREX_USE_OMP
+#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
+#endif
+            for (ParIt pti(dest, lev); pti.isValid(); ++pti) {
+                const int np = pti.numParticles();
+
+                // preparing access to particle data: AoS
+                using PType = ImpactXParticleContainer::ParticleType;
+                auto &aos = pti.GetArrayOfStructs();
+                PType *AMREX_RESTRICT aos_ptr = aos().dataPtr();
+
+                amrex::ParallelFor(np, [=] AMREX_GPU_DEVICE(long i) {
+                    PType &p = aos_ptr[i];
+                    p.id() = -p.id();
+                });
+            }
+        }
+
+        // TODO remove particles with negative ids in source
+    }
+} // namespace impactx

src/particles/elements/Aperture.H

@@ -88,16 +88,12 @@ namespace impactx
             {
                 case Shape::rectangular :  // default
                   if (pow(u,2)>1 || pow(v,2) > 1_prt) {
-                     p.pos(RealAoS::x) = 0.0_prt; // replace with id change of sign
-                     p.pos(RealAoS::y) = 0.0_prt; // replace with id change of sign
                      p.id() = -id;
                   }
                   break;
 
                case Shape::elliptical :
                   if (pow(u,2)+pow(v,2) > 1_prt) {
-                     p.pos(RealAoS::x) = 0.0_prt; // replace with id change of sign
-                     p.pos(RealAoS::y) = 0.0_prt; // replace with id change of sign
                      p.id() = -id;
                   }
                   break;