Bugfix to particle defragmentation

CHANGELOG.md

@@ -3,6 +3,7 @@
 ## Current develop
 
 ### Added (new features/APIs/variables/...)
+- [[PR 1185]](https://github.com/parthenon-hpc-lab/parthenon/pull/1185/files) Bugfix to particle defragmentation
 - [[PR 1183]](https://github.com/parthenon-hpc-lab/parthenon/pull/1183) Fix particle leapfrog example initialization data
 - [[PR 1179]](https://github.com/parthenon-hpc-lab/parthenon/pull/1179) Make a global variable for whether simulation is a restart
 - [[PR 1171]](https://github.com/parthenon-hpc-lab/parthenon/pull/1171) Add PARTHENON_USE_SYSTEM_PACKAGES build option

example/particles/parthinput.particles

@@ -59,3 +59,4 @@ rng_seed = 23487
 const_dt = 0.5
 deposition_method = per_cell
 destroy_particles_frac = 0.1
+defrag_threshold = 0.9

example/particles/particles.cpp

@@ -62,6 +62,11 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
       destroy_particles_frac >= 0. && destroy_particles_frac <= 1.,
       "Fraction of particles to destroy each timestep must be between 0 and 1");
 
+  Real defrag_threshold = pin->GetOrAddReal("Particles", "defrag_threshold", 0.9);
+  pkg->AddParam<>("defrag_threshold", defrag_threshold);
+  PARTHENON_REQUIRE(defrag_threshold >= 0. && defrag_threshold <= 1.,
+                    "Defragmentation threshold must be between 0 and 1 inclusive.");
+
   std::string deposition_method =
       pin->GetOrAddString("Particles", "deposition_method", "per_particle");
   if (deposition_method == "per_particle") {
@@ -609,6 +614,9 @@ TaskCollection ParticleDriver::MakeFinalizationTaskCollection() const {
     auto &sc1 = pmb->meshblock_data.Get();
     auto &tl = async_region1[i];
 
+    auto pkg = pmb->packages.Get("particles_package");
+    const auto defrag_threshold = pkg->Param<Real>("defrag_threshold");
+
     auto destroy_some_particles = tl.AddTask(none, DestroySomeParticles, pmb.get());
 
     auto sort_particles = tl.AddTask(destroy_some_particles,
@@ -617,7 +625,8 @@ TaskCollection ParticleDriver::MakeFinalizationTaskCollection() const {
     auto deposit_particles = tl.AddTask(sort_particles, DepositParticles, pmb.get());
 
     // Defragment if swarm memory pool occupancy is 90%
-    auto defrag = tl.AddTask(deposit_particles, &SwarmContainer::Defrag, sc.get(), 0.9);
+    auto defrag = tl.AddTask(deposit_particles, &SwarmContainer::Defrag, sc.get(),
+                             defrag_threshold);
 
     // estimate next time step
     auto new_dt = tl.AddTask(

src/interface/swarm.cpp

@@ -354,15 +354,18 @@ void Swarm::Defrag() {
   auto &mask = mask_;
 
   const int &num_active = num_active_;
+  auto empty_indices = empty_indices_;
   parthenon::par_scan(
-      "Set empty indices prefix sum", 0, nmax_pool_ - num_active_ - 1,
-      KOKKOS_LAMBDA(const int nn, int &update, const bool &final) {
-        const int n = nn + num_active;
+      "Set empty indices prefix sum", num_active, nmax_pool_ - 1,
+      KOKKOS_LAMBDA(const int n, int &update, const bool &final) {
         const int val = mask(n);
         if (val) {
           update += 1;
         }
-        if (final) scan_scratch_toread(n) = update;
+        if (final) {
+          scan_scratch_toread(n) = update;
+          empty_indices(n - num_active) = n;
+        }
       });
 
   parthenon::par_for(