Use 63 bits for Beam ID

docs/source/run/parameters.rst

@@ -514,9 +514,6 @@ which are valid only for certain beam types, are introduced further below under
     Whether the beam particles are pushed along the z-axis. The momentum is still fully updated.
     Note: using ``do_z_push = 0`` results in unphysical behavior.
 
-* ``<beam name> or beams.do_reset_id_init`` (`bool`) optional (default `0`)
-    Whether to reset the ID incrementor to 1 before initializing beam particles.
-
 * ``<beam name> or beams.reorder_period`` (`int`) optional (default `0`)
     Reorder particles periodically to speed-up current deposition and particle push on GPU.
     A good starting point is a period of 1 to reorder beam particles on every timestep.

src/diagnostics/OpenPMDWriter.cpp

@@ -249,7 +249,9 @@ OpenPMDWriter::WriteBeamParticleData (MultiBeam& beams, openPMD::Iteration itera
             uint64_t * const uint64_data = m_uint64_beam_data[ibeam][idx].get();
 
             for (uint64_t i=0; i<np_total; ++i) {
-                amrex::Long id = amrex::ConstParticleIDWrapper(uint64_data[i]);
+                uint64_t id = uint64_data[i];
+                // in the amrex format valid idcpus start with 1 and invalid with 0
+                amrex::ParticleIDWrapper{id}.make_invalid();
                 uint64_data[i] = id;
             }
 

src/particles/beam/BeamParticleContainer.H

@@ -35,6 +35,7 @@ struct BeamIdx
         // nsubcycles: by how many subcycles was this particle pushed already
         // nsubcycles is not stored or communicated in MultiBuffer
         nsubcycles=int_nattribs_in_buffer,
+        mr_level,
         int_nattribs
     };
 };
@@ -202,7 +203,7 @@ public:
         return rcomp < numRealComponents();
     }
     bool communicateIntComponent (int icomp) const {
-        // don't communicate nsubcycles
+        // don't communicate nsubcycles or mr_level
         return icomp < BeamIdx::int_nattribs_in_buffer;
     }
 
@@ -222,8 +223,6 @@ public:
     bool m_do_radiation_reaction {false}; /**< whether to calculate radiation losses */
     /** Number of particles on upstream rank (required for IO) */
     bool m_do_salame = false; /**< Whether this beam uses salame */
-    /** Whether to reset the global ID incrementor to 1 before initializing this beam */
-    bool m_do_reset_id_init {false};
     bool m_initialize_on_cpu {false};
     /** How often the insitu beam diagnostics should be computed and written
      * Default is 0, meaning no output */

src/particles/beam/BeamParticleContainer.cpp

@@ -63,7 +63,6 @@ BeamParticleContainer::ReadParameters ()
     queryWithParserAlt(pp, "insitu_radius", m_insitu_radius, pp_alt);
     queryWithParser(pp, "n_subcycles", m_n_subcycles);
     AMREX_ALWAYS_ASSERT_WITH_MESSAGE( m_n_subcycles >= 1, "n_subcycles must be >= 1");
-    queryWithParserAlt(pp, "do_reset_id_init", m_do_reset_id_init, pp_alt);
     queryWithParser(pp, "do_salame", m_do_salame);
     queryWithParserAlt(pp, "reorder_period", m_reorder_period, pp_alt);
     amrex::Array<int, 2> idx_array
@@ -121,7 +120,6 @@ BeamParticleContainer::InitData (const amrex::Geometry& geom)
     amrex::ParmParse pp(m_name);
     amrex::ParmParse pp_alt("beams");
     amrex::Real ptime {0.};
-    if (m_do_reset_id_init) BeamTileInit::ParticleType::NextID(1);
     if (m_injection_type == "fixed_ppc") {
 
         queryWithParser(pp, "ppc", m_ppc);
@@ -315,7 +313,7 @@ void BeamParticleContainer::TagByLevel (const int current_N_level,
     auto& soa = getBeamSlice(which_slice).GetStructOfArrays();
     const amrex::Real * const pos_x = soa.GetRealData(BeamIdx::x).data();
     const amrex::Real * const pos_y = soa.GetRealData(BeamIdx::y).data();
-    auto * const idcpup = soa.GetIdCPUData().data();
+    int * const p_mr_level = soa.GetIntData(BeamIdx::mr_level).data();
 
     const int lev1_idx = std::min(1, current_N_level-1);
     const int lev2_idx = std::min(2, current_N_level-1);
@@ -341,15 +339,15 @@ void BeamParticleContainer::TagByLevel (const int current_N_level,
                 lo_x_lev2 < xp && xp < hi_x_lev2 &&
                 lo_y_lev2 < yp && yp < hi_y_lev2) {
                 // level 2
-                amrex::ParticleCPUWrapper{idcpup[ip]} = 2;
+                p_mr_level[ip] = 2;
             } else if (current_N_level > 1 &&
                 lo_x_lev1 < xp && xp < hi_x_lev1 &&
                 lo_y_lev1 < yp && yp < hi_y_lev1) {
                 // level 1
-                amrex::ParticleCPUWrapper{idcpup[ip]} = 1;
+                p_mr_level[ip] = 1;
             } else {
                 // level 0
-                amrex::ParticleCPUWrapper{idcpup[ip]} = 0;
+                p_mr_level[ip] = 0;
             }
         }
     );

src/particles/beam/BeamParticleContainerInit.cpp

@@ -41,8 +41,8 @@ namespace
     void AddOneBeamParticle (
         const BeamTileInit::ParticleTileDataType& ptd, const amrex::Real& x,
         const amrex::Real& y, const amrex::Real& z, const amrex::Real& ux, const amrex::Real& uy,
-        const amrex::Real& uz, const amrex::Real& weight, const int& pid,
-        const int& ip, const amrex::Real& speed_of_light) noexcept
+        const amrex::Real& uz, const amrex::Real& weight, const amrex::Long pid,
+        const amrex::Long ip, const amrex::Real& speed_of_light) noexcept
     {
         ptd.rdata(BeamIdx::x  )[ip] = x;
         ptd.rdata(BeamIdx::y  )[ip] = y;
@@ -52,8 +52,8 @@ namespace
         ptd.rdata(BeamIdx::uz )[ip] = uz * speed_of_light;
         ptd.rdata(BeamIdx::w  )[ip] = std::abs(weight);
 
-        ptd.id(ip) = pid > 0 ? pid + ip : pid;
-        ptd.cpu(ip) = 0;
+        ptd.idcpu(ip) = pid + ip;
+        ptd.id(ip).make_valid();
     }
 
     /** \brief Adds a single beam particle into the per-slice BeamTile
@@ -86,15 +86,15 @@ namespace
         ptd.rdata(BeamIdx::uz )[ip] = uz * speed_of_light;
         ptd.rdata(BeamIdx::w  )[ip] = std::abs(weight);
 
-        // ensure id is always valid
-        // it will repeat if there are more than 2^31-1 particles
-        int id = int((pid + ip) & ((1ull << 31) - 1));
-        if (id == 0) id = 1;
-        if (!is_valid) id = -1;
+        ptd.idcpu(ip) = pid + ip;
+        if (is_valid) {
+            ptd.id(ip).make_valid(); // ensure id is valid
+        } else {
+            ptd.id(ip).make_invalid();
+        }
 
-        ptd.id(ip) = id;
-        ptd.cpu(ip) = 0;
         ptd.idata(BeamIdx::nsubcycles)[ip] = 0;
+        ptd.idata(BeamIdx::mr_level)[ip] = 0;
     }
 }
 
@@ -389,7 +389,7 @@ InitBeamFixedWeightSlice (int slice, int which_slice)
     amrex::Real * const pos_z = m_z_array.dataPtr();
 
     const uint64_t pid = m_id64;
-    m_id64 += num_to_add;
+    m_id64 += m_do_symmetrize ? 4*num_to_add : num_to_add;
 
     const bool can = m_can_profile;
     const bool do_symmetrize = m_do_symmetrize;
@@ -587,7 +587,7 @@ InitBeamFixedWeightPDFSlice (int slice, int which_slice)
         const auto ptd = particle_tile.getParticleTileData();
 
         const uint64_t pid = m_id64;
-        m_id64 += num_to_add;
+        m_id64 += m_do_symmetrize ? 4*num_to_add : num_to_add;
 
         const amrex::Real clight = get_phys_const().c;
         const bool do_symmetrize = m_do_symmetrize;
@@ -947,7 +947,7 @@ InitBeamFromFile (const std::string input_file,
     const auto num_to_add = electrons[name_r][name_rx].getExtent()[0];
 
     if(num_to_add >= 2147483647) {
-        amrex::Abort("Beam can't have more than 2'147'483'646 Particles\n");
+        amrex::Abort("Beam from file can't have more than 2'147'483'646 Particles\n");
     }
 
     auto del = [](input_type *p){ amrex::The_Pinned_Arena()->free(reinterpret_cast<void*>(p)); };
@@ -1031,8 +1031,9 @@ InitBeamFromFile (const std::string input_file,
     auto new_size = old_size + num_to_add;
     particle_tile.resize(new_size);
     const auto ptd = particle_tile.getParticleTileData();
-    const int pid = BeamTileInit::ParticleType::NextID();
-    BeamTileInit::ParticleType::NextID(pid + num_to_add);
+
+    const uint64_t pid = m_id64;
+    m_id64 += num_to_add;
 
     const input_type * const r_x_ptr = r_x_data.get();
     const input_type * const r_y_ptr = r_y_data.get();
@@ -1042,8 +1043,8 @@ InitBeamFromFile (const std::string input_file,
     const input_type * const u_z_ptr = u_z_data.get();
     const input_type * const w_w_ptr = w_w_data.get();
 
-    amrex::ParallelFor(int(num_to_add),
-        [=] AMREX_GPU_DEVICE (const int i) {
+    amrex::ParallelFor(amrex::Long(num_to_add),
+        [=] AMREX_GPU_DEVICE (const amrex::Long i) {
             AddOneBeamParticle(ptd,
                 static_cast<amrex::Real>(r_x_ptr[i] * unit_rx),
                 static_cast<amrex::Real>(r_y_ptr[i] * unit_ry),

src/particles/deposition/BeamDepositCurrent.cpp

@@ -99,7 +99,11 @@ DepositCurrentSlice (BeamParticleContainer& beam, Fields& fields,
 
             // Skip invalid particles and ghost particles not in the last slice
             // beam deposits only up to its finest level
-            if (!ptd.id(ip).is_valid() || (only_highest ? (ptd.cpu(ip)!=lev) : (ptd.cpu(ip)<lev))) return;
+            if (!ptd.id(ip).is_valid() ||
+                (only_highest ?
+                    (ptd.idata(BeamIdx::mr_level)[ip] != lev) :
+                    (ptd.idata(BeamIdx::mr_level)[ip] < lev))
+            ) return;
 
             // --- Get particle quantities
             const amrex::Real ux = ptd.rdata(BeamIdx::ux)[ip];

src/particles/plasma/PlasmaParticleContainer.cpp

@@ -396,9 +396,6 @@ IonizationModule (const int lev,
         ptile_elec.resize(new_size);
 
         // Load electron soa and aos after resize
-        const long pid_start = ParticleType::NextID();
-        ParticleType::NextID(pid_start + num_new_electrons.dataValue());
-
         auto arrdata_ion = ptile_ion.GetStructOfArrays().realarray();
         auto arrdata_elec = ptile_elec.GetStructOfArrays().realarray();
         auto int_arrdata_elec = ptile_elec.GetStructOfArrays().intarray();
@@ -417,7 +414,7 @@ IonizationModule (const int lev,
                 const long pidx = pid + old_size;
 
                 // Copy ion data to new electron
-                amrex::ParticleIDWrapper{idcpu_elec[pidx]} = pid_start + pid;
+                amrex::ParticleIDWrapper{idcpu_elec[pidx]} = 2; // only for valid/invalid
                 amrex::ParticleCPUWrapper{idcpu_elec[pidx]} = lev; // current level
                 arrdata_elec[PlasmaIdx::x      ][pidx] = arrdata_ion[PlasmaIdx::x     ][ip];
                 arrdata_elec[PlasmaIdx::y      ][pidx] = arrdata_ion[PlasmaIdx::y     ][ip];

src/particles/plasma/PlasmaParticleContainerInit.cpp

@@ -59,9 +59,6 @@ InitParticles (const amrex::RealVect& a_u_std,
         y_offset = -0.5_rt;
     }
 
-    // Reset first ID to 1 as we create a fresh plasma each time we call this function
-    ParticleType::NextID(1);
-
     for(amrex::MFIter mfi = MakeMFIter(lev, DfltMfi); mfi.isValid(); ++mfi)
     {
         amrex::Box tile_box  = mfi.tilebox(box_nodal, box_grow);
@@ -182,9 +179,6 @@ InitParticles (const amrex::RealVect& a_u_std,
         const auto new_size = old_size + total_num_particles;
         particle_tile.resize(new_size);
 
-        const int pid = ParticleType::NextID();
-        ParticleType::NextID(pid + total_num_particles);
-
         auto ptd = particle_tile.getParticleTileData();
         const int init_ion_lev = m_init_ion_lev;
 
@@ -301,7 +295,7 @@ InitParticles (const amrex::RealVect& a_u_std,
                 amrex::Real u[3] = {0.,0.,0.};
                 ParticleUtil::get_gaussian_random_momentum(u, a_u_mean, a_u_std, engine);
 
-                ptd.id(pidx) = pid + int(pidx);
+                ptd.id(pidx) = 1; // plasma id is only used to distinguish between valid/invalid
                 ptd.cpu(pidx) = 0; // level 0
                 ptd.rdata(PlasmaIdx::x)[pidx] = x;
                 ptd.rdata(PlasmaIdx::y)[pidx] = y;
@@ -358,7 +352,7 @@ InitParticles (const amrex::RealVect& a_u_std,
                 for (int imirror=0; imirror<3; ++imirror) {
                     const amrex::Long midx = (imirror+1)*mirror_offset +pidx;
 
-                    ptd.id(midx) = pid + int(midx);
+                    ptd.id(midx) = 1; // plasma id is only used to distinguish between valid/invalid
                     ptd.cpu(midx) = 0; // level 0
                     ptd.rdata(PlasmaIdx::x)[midx] = x_arr[imirror];
                     ptd.rdata(PlasmaIdx::y)[midx] = y_arr[imirror];

tests/checksum/benchmarks_json/production.SI.2Rank_pwfa.json

@@ -19,7 +19,7 @@
   },
   "driver": {
     "charge": 1.601913877032e-13,
-    "id": 137641610358,
+    "id": 499836513366,
     "mass": 9.107889762573e-25,
     "x": 3.1337620903426,
     "y": 3.1349173868172,
@@ -31,7 +31,7 @@
   },
   "witness": {
     "charge": 1.602176634e-13,
-    "id": 199678708984,
+    "id": 500000500000,
     "mass": 9.1093837015e-25,
     "x": 2.457586712044,
     "y": 2.4574680947762,