Fix geometry initialization when new primaries are inserted in the mi…

…ddle of a simulation (#1369)
celeritas-project · Aug 16, 2024 · 9ed50a6 · 9ed50a6
1 parent 1e4aa73
commit 9ed50a6
Show file tree

Hide file tree

Showing 5 changed files with 22 additions and 5 deletions.
diff --git a/src/celeritas/track/ExtendFromPrimariesAction.cc b/src/celeritas/track/ExtendFromPrimariesAction.cc
@@ -9,6 +9,7 @@
 
 #include "corecel/Assert.hh"
 #include "corecel/Macros.hh"
+#include "corecel/data/CollectionAlgorithms.hh"
 #include "corecel/sys/MultiExceptionHandler.hh"
 #include "celeritas/global/ActionLauncher.hh"
 #include "celeritas/global/CoreParams.hh"
@@ -58,6 +59,12 @@ void ExtendFromPrimariesAction::execute_impl(CoreParams const& params,
 
     // Mark that the primaries have been processed
     state.clear_primaries();
+
+    // Clear the track slot IDs of the track initializers' parent tracks. This
+    // is necessary when new primaries are inserted in the middle of a
+    // simulation and the parent IDs of secondaries produced in the previous
+    // step have been stored.
+    fill(TrackSlotId{}, &state.ref().init.parents);
 }
 
 //---------------------------------------------------------------------------//

diff --git a/src/celeritas/track/detail/InitTracksExecutor.hh b/src/celeritas/track/detail/InitTracksExecutor.hh
@@ -94,13 +94,21 @@ CELER_FUNCTION void InitTracksExecutor::operator()(ThreadId tid) const
     // Initialize the geometry
     {
         auto geo = vacancy.make_geo_view();
-        if (tid < counters.num_secondaries)
+        auto parent_id = [&] {
+            if (!(tid < counters.num_secondaries))
+            {
+                return TrackSlotId{};
+            }
+            return data
+                .parents[TrackSlotId(index_before(data.parents.size(), tid))];
+        }();
+
+        if (parent_id)
         {
             // Copy the geometry state from the parent for improved performance
-            TrackSlotId parent_id = data.parents[TrackSlotId{
-                index_before(data.parents.size(), tid)}];
             GeoTrackView const parent_geo(
                 params->geometry, state->geometry, parent_id);
+            CELER_ASSERT(parent_geo.pos() == init.geo.pos);
             geo = GeoTrackView::DetailedInitializer{parent_geo, init.geo.dir};
             CELER_ASSERT(!geo.is_outside());
         }

diff --git a/src/corecel/data/detail/Filler.cu b/src/corecel/data/detail/Filler.cu
@@ -15,6 +15,7 @@ namespace detail
 template struct Filler<real_type, MemSpace::device>;
 template struct Filler<size_type, MemSpace::device>;
 template struct Filler<int, MemSpace::device>;
+template struct Filler<TrackSlotId, MemSpace::device>;
 //---------------------------------------------------------------------------//
 }  // namespace detail
 }  // namespace celeritas
diff --git a/src/corecel/data/detail/Filler.hh b/src/corecel/data/detail/Filler.hh
@@ -59,6 +59,7 @@ void Filler<T, MemSpace::device>::operator()(Span<T>) const
 extern template struct Filler<real_type, MemSpace::device>;
 extern template struct Filler<size_type, MemSpace::device>;
 extern template struct Filler<int, MemSpace::device>;
+extern template struct Filler<TrackSlotId, MemSpace::device>;
 #endif
 
 //---------------------------------------------------------------------------//

diff --git a/test/celeritas/global/StepperGeant.test.cc b/test/celeritas/global/StepperGeant.test.cc
@@ -391,8 +391,8 @@ TEST_F(TestEm3NoMsc, host_multi)
     counts = step();
     if (this->is_default_build())
     {
-        EXPECT_EQ(44, counts.active);
-        EXPECT_EQ(43, counts.alive);
+        EXPECT_EQ(36, counts.active);
+        EXPECT_EQ(35, counts.alive);
     }
 }