Fix multithreading in Geant4 and HitManager (#694)

* Make StepCollector thread-safe
* Create thread-local hit processors
* Add debug log messages to sensitive detectors
* Extend InitializedValue to do finalization
* Use InitializedValue to correctly finalize local data

app/CMakeLists.txt

@@ -431,9 +431,8 @@ if(CELERITAS_BUILD_DEMOS AND CELERITAS_USE_Geant4)
       "CELER_DISABLE_PARALLEL=1"
     )
     if(Geant4_multithreaded_FOUND)
-      # TODO: use 2 threads when we support multithreaded hit manager
       list(APPEND _env
-        "G4FORCENUMBEROFTHREADS=1"
+        "G4FORCENUMBEROFTHREADS=2"
       )
     endif()
     if(CELERITAS_USE_OpenMP)

src/accel/LocalTransporter.cc

@@ -23,6 +23,7 @@
 
 #include "SetupOptions.hh"
 #include "SharedParams.hh"
+#include "detail/HitManager.hh"
 
 namespace celeritas
 {
@@ -50,7 +51,9 @@ inline Real3 convert_from_geant(G4ThreeVector const& vec, double units)
  */
 LocalTransporter::LocalTransporter(SetupOptions const& options,
                                    SharedParams const& params)
-    : auto_flush_(options.max_num_tracks), max_steps_(options.max_steps)
+    : auto_flush_(options.max_num_tracks)
+    , max_steps_(options.max_steps)
+    , hit_manager_{params.hit_manager()}
 {
     CELER_EXPECT(params);
     particles_ = params.Params()->particle();
@@ -202,5 +205,17 @@ void LocalTransporter::Finalize()
     CELER_ENSURE(!*this);
 }
 
+//---------------------------------------------------------------------------//
+/*!
+ * Clear thread-local hit manager on destruction.
+ */
+void LocalTransporter::HMFinalizer::operator()(SPHitManger& hm) const
+{
+    if (hm)
+    {
+        hm->finalize();
+    }
+}
+
 //---------------------------------------------------------------------------//
 }  // namespace celeritas

src/accel/LocalTransporter.hh

@@ -12,6 +12,7 @@
 #include <G4Track.hh>
 
 #include "corecel/Types.hh"
+#include "corecel/cont/InitializedValue.hh"
 #include "corecel/io/Logger.hh"
 #include "celeritas/Types.hh"
 #include "celeritas/global/CoreParams.hh"
@@ -20,6 +21,11 @@
 
 namespace celeritas
 {
+namespace detail
+{
+class HitManager;
+}
+
 struct SetupOptions;
 class SharedParams;
 
@@ -69,6 +75,13 @@ class LocalTransporter
     explicit operator bool() const { return static_cast<bool>(step_); }
 
   private:
+    using SPHitManger = std::shared_ptr<detail::HitManager>;
+
+    struct HMFinalizer
+    {
+        void operator()(SPHitManger& hm) const;
+    };
+
     std::shared_ptr<ParticleParams const> particles_;
     std::shared_ptr<StepperInterface> step_;
     std::vector<Primary> buffer_;
@@ -78,6 +91,9 @@ class LocalTransporter
 
     size_type auto_flush_{};
     size_type max_steps_{};
+
+    // Shared pointer across threads, "finalize" called when clearing
+    InitializedValue<SPHitManger, HMFinalizer> hit_manager_;
 };
 
 //---------------------------------------------------------------------------//

src/accel/SharedParams.cc

@@ -375,8 +375,6 @@ void SharedParams::initialize_core(SetupOptions const& options)
     // Construct sensitive detector callback
     if (options.sd)
     {
-        CELER_VALIDATE(params.max_streams == 1,
-                       << "HitManager currently supports only serial mode");
         hit_manager_ = std::make_shared<detail::HitManager>(*params.geometry,
                                                             options.sd);
         step_collector_ = std::make_shared<StepCollector>(

src/accel/SharedParams.hh

@@ -41,6 +41,7 @@ class SharedParams
     //!@{
     //! \name Type aliases
     using SPConstParams = std::shared_ptr<CoreParams const>;
+    using SPHitManager = std::shared_ptr<detail::HitManager>;
     //!@}
 
   public:
@@ -70,6 +71,9 @@ class SharedParams
     //! Whether this instance is initialized
     explicit operator bool() const { return static_cast<bool>(params_); }
 
+    //! Hit manager, to be used only by LocalTransporter
+    SPHitManager const& hit_manager() const { return hit_manager_; }
+
   private:
     //// DATA ////
 

src/accel/detail/HitManager.cc

@@ -10,13 +10,17 @@
 #include <utility>
 #include <G4LogicalVolume.hh>
 #include <G4LogicalVolumeStore.hh>
+#include <G4RunManager.hh>
+#include <G4Threading.hh>
+#include <G4VSensitiveDetector.hh>
 
 #include "celeritas_cmake_strings.h"
 #include "corecel/cont/EnumArray.hh"
 #include "corecel/cont/Label.hh"
 #include "corecel/cont/Range.hh"
 #include "corecel/io/Logger.hh"
 #include "celeritas/Types.hh"
+#include "celeritas/ext/GeantSetup.hh"
 #include "celeritas/ext/detail/GeantVolumeVisitor.hh"
 #include "celeritas/geo/GeoParams.hh"  // IWYU pragma: keep
 #include "celeritas/io/ImportVolume.hh"
@@ -47,21 +51,22 @@ void update_selection(StepPointSelection* selection,
  */
 HitManager::HitManager(GeoParams const& geo, SDSetupOptions const& setup)
     : nonzero_energy_deposition_(setup.ignore_zero_deposition)
+    , locate_touchable_(setup.locate_touchable)
 {
     CELER_EXPECT(setup.enabled);
 
     // Convert setup options to step data
     selection_.energy_deposition = setup.energy_deposition;
     update_selection(&selection_.points[StepPoint::pre], setup.pre);
     update_selection(&selection_.points[StepPoint::post], setup.post);
-    if (setup.locate_touchable)
+    if (locate_touchable_)
     {
         selection_.points[StepPoint::pre].pos = true;
         selection_.points[StepPoint::pre].dir = true;
     }
 
     // Logical volumes to pass to hit processor
-    std::vector<G4LogicalVolume*> lv_with_sd;
+    std::vector<G4LogicalVolume*> geant_vols;
 
     // Helper class to extract GDML names+labels from Geant4 volume
     GeantVolumeVisitor visitor(true);
@@ -115,16 +120,29 @@ HitManager::HitManager(GeoParams const& geo, SDSetupOptions const& setup)
                        << "failed to find " << celeritas_geometry
                        << " volume corresponding to Geant4 volume '"
                        << lv->GetName() << "'");
+        CELER_LOG(debug) << "Mapped sensitive detector '" << sd->GetName()
+                         << "' (logical volume '" << lv->GetName() << "') to "
+                         << celeritas_geometry << " volume '"
+                         << geo.id_to_label(id) << "'";
 
         // Add Geant4 volume and corresponding volume ID to list
-        lv_with_sd.push_back(lv);
+        geant_vols.push_back(lv);
         vecgeom_vols_.push_back(id);
     }
     CELER_VALIDATE(!vecgeom_vols_.empty(),
                    << "no sensitive detectors were found");
 
-    process_hits_ = std::make_unique<HitProcessor>(
-        std::move(lv_with_sd), selection_, setup.locate_touchable);
+    // Hit processors *must* be allocated on the thread they're used because of
+    // geant4 thread-local SDs. There must be one per thread.
+    auto* run_man = G4RunManager::GetRunManager();
+    CELER_VALIDATE(run_man,
+                   << "G4RunManager was not created before setting up "
+                      "HitManager");
+    processors_.resize(celeritas::get_num_threads(*run_man));
+
+    geant_vols_ = std::make_shared<std::vector<G4LogicalVolume*>>(
+        std::move(geant_vols));
+    CELER_ENSURE(geant_vols_ && geant_vols_->size() == vecgeom_vols_.size());
 }
 
 //---------------------------------------------------------------------------//
@@ -155,11 +173,8 @@ auto HitManager::filters() const -> Filters
  */
 void HitManager::execute(StateHostRef const& data)
 {
-    copy_steps(&steps_, data);
-    if (steps_)
-    {
-        (*process_hits_)(steps_);
-    }
+    auto&& process_hits = this->get_local_hit_processor();
+    process_hits(data);
 }
 
 //---------------------------------------------------------------------------//
@@ -168,11 +183,43 @@ void HitManager::execute(StateHostRef const& data)
  */
 void HitManager::execute(StateDeviceRef const& data)
 {
-    copy_steps(&steps_, data);
-    if (steps_)
+    auto&& process_hits = this->get_local_hit_processor();
+    process_hits(data);
+}
+
+//---------------------------------------------------------------------------//
+/*!
+ * Destroy local data to avoid Geant4 crashes.
+ *
+ * This deallocates the local hit processor on the Geant4 thread that was
+ * active when allocating it. This is necessary because Geant4 has thread-local
+ * allocators that crash if trying to deallocate data allocated on another
+ * thread.
+ */
+void HitManager::finalize()
+{
+    CELER_LOG_LOCAL(debug) << "Deallocating hit processor";
+    int local_thread = G4Threading::G4GetThreadId();
+    CELER_ASSERT(static_cast<std::size_t>(local_thread) < processors_.size());
+    processors_[local_thread].reset();
+}
+
+//---------------------------------------------------------------------------//
+/*!
+ * Ensure the local hit processor exists, and return it.
+ */
+HitProcessor& HitManager::get_local_hit_processor()
+{
+    int local_thread = G4Threading::G4GetThreadId();
+    CELER_ASSERT(static_cast<std::size_t>(local_thread) < processors_.size());
+    if (CELER_UNLIKELY(!processors_[local_thread]))
     {
-        (*process_hits_)(steps_);
+        CELER_LOG_LOCAL(debug) << "Allocating hit processor";
+        // Allocate the hit processor locally
+        processors_[local_thread] = std::make_unique<HitProcessor>(
+            geant_vols_, selection_, locate_touchable_);
     }
+    return *processors_[local_thread];
 }
 
 //---------------------------------------------------------------------------//

src/accel/detail/HitManager.hh

@@ -12,9 +12,10 @@
 
 #include "orange/Types.hh"
 #include "celeritas/geo/GeoParamsFwd.hh"
-#include "celeritas/user/DetectorSteps.hh"
 #include "celeritas/user/StepInterface.hh"
 
+class G4LogicalVolume;
+
 namespace celeritas
 {
 struct SDSetupOptions;
@@ -31,14 +32,9 @@ class HitProcessor;
  * - Created during SharedParams::Initialize alongside the step collector
  * - Is shared across threads
  * - Finds *all* logical volumes that have SDs attached (TODO: add list of
- *   exclusions?)
+ *   exclusions for SDs that are implemented natively on GPU)
  * - Maps those volumes to VecGeom geometry
  * - Creates a HitProcessor for each Geant4 thread
- *
- * Execute:
- * - Can share DetectorStepOutput across threads for now since StepGatherAction
- *   is mutexed across all threads
- * - Calls a single HitProcessor (thread safe because of caller's mutex)
  */
 class HitManager final : public StepInterface
 {
@@ -61,14 +57,19 @@ class HitManager final : public StepInterface
     // Process device-generated hits
     void execute(StateDeviceRef const&) final;
 
+    // Destroy local data to avoid Geant4 crashes
+    void finalize();
+
   private:
     bool nonzero_energy_deposition_{};
+    bool locate_touchable_{};
     StepSelection selection_;
-    DetectorStepOutput steps_;
+    std::shared_ptr<const std::vector<G4LogicalVolume*>> geant_vols_;
     std::vector<VolumeId> vecgeom_vols_;
-    std::unique_ptr<HitProcessor> process_hits_;
+    std::vector<std::unique_ptr<HitProcessor>> processors_;
 
-    void call_local_processor() const;
+    // Ensure thread-local hit processor exists and return it
+    HitProcessor& get_local_hit_processor();
 };
 
 //---------------------------------------------------------------------------//