Fix bug in stalling of the GPUs when running out of hitslots (#460)

This fixes a bug in the stalling of the GPU in case the hitBuffer is
overflowing. Since the calculation of `nextStepMightFail` is based on
unsigned int parameters, the subtraction was causing false negatives in
case of wraparounds if `hitBufferSafetyFactor * maxInFlight ` was larger
than ` gpuState.fHitScoring->HitCapacity() / numThread`. By changing the
condition to an addition
```
        bool nextStepMightFail = gpuState.stats->hitBufferOccupancy + hitBufferSafetyFactor * maxInFlight >=
                                 gpuState.fHitScoring->HitCapacity() / numThreads;
```
this is safely resolved.

Furthermore, as a track can cause multiple (3) hits per track, (one for
each secondary, one for the energy it deposits itself), the previous
hardcoded hitBufferSafetyFactor of 1.5 was dangerous, so now it is an UI
parameter. It doesn't make sense to hardcode a higher (safe) value, as
it would slow down the settings where a lower factor is safe.

Furthermore, two printouts that are per track are increased to a higher
verbosity, as those render larger simulations useless, but other
printouts (e.g., whether the hitbuffer is copied) were at the same
level.

Author: SeverinDiederichs

include/AdePT/core/AdePTConfiguration.hh

@@ -42,6 +42,7 @@ public:
   void SetMillionsOfHitSlots(double millionSlots) { fMillionsOfHitSlots = millionSlots; }
   void SetHitBufferFlushThreshold(float threshold) { fHitBufferFlushThreshold = threshold; }
   void SetCPUCapacityFactor(float CPUCapacityFactor) { fCPUCapacityFactor = CPUCapacityFactor; }
+  void SetHitBufferSafetyFactor(double HitBufferSafetyFactor) { fHitBufferSafetyFactor = HitBufferSafetyFactor; }
   void SetAdePTSeed(int seed) { fAdePTSeed = static_cast<uint64_t>(seed); }
 
   void SetCUDAStackLimit(int limit) { fCUDAStackLimit = limit; }
@@ -81,6 +82,7 @@ public:
   double GetWDTKineticEnergyLimit() { return fWDTKineticEnergyLimit; }
   float GetHitBufferFlushThreshold() { return fHitBufferFlushThreshold; }
   float GetCPUCapacityFactor() { return fCPUCapacityFactor; }
+  double GetHitBufferSafetyFactor() { return fHitBufferSafetyFactor; }
   double GetMillionsOfTrackSlots() { return fMillionsOfTrackSlots; }
   double GetMillionsOfLeakSlots() { return fMillionsOfLeakSlots; }
   double GetMillionsOfHitSlots() { return fMillionsOfHitSlots; }
@@ -108,6 +110,7 @@ private:
   uint64_t fAdePTSeed{1234567};
   float fHitBufferFlushThreshold{0.8};
   float fCPUCapacityFactor{2.5};
+  double fHitBufferSafetyFactor{1.5};
   double fMillionsOfTrackSlots{1};
   double fMillionsOfLeakSlots{1};
   double fMillionsOfHitSlots{1};

include/AdePT/core/AsyncAdePTTransport.cuh

@@ -912,7 +912,8 @@ void HitProcessingLoop(HitProcessingContext *const context, GPUstate &gpuState,
 void TransportLoop(int trackCapacity, int leakCapacity, int scoringCapacity, int numThreads, TrackBuffer &trackBuffer,
                    GPUstate &gpuState, std::vector<std::atomic<EventState>> &eventStates,
                    std::condition_variable &cvG4Workers, std::vector<AdePTScoring> &scoring, int adeptSeed,
-                   int debugLevel, bool returnAllSteps, bool returnLastStep, unsigned short lastNParticlesOnCPU)
+                   int debugLevel, bool returnAllSteps, bool returnLastStep, unsigned short lastNParticlesOnCPU,
+                   const double hitBufferSafetyFactor)
 {
   // NVTXTracer tracer{"TransportLoop"};
 
@@ -1624,9 +1625,11 @@ void TransportLoop(int trackCapacity, int leakCapacity, int scoringCapacity, int
                        "/adept/setMillionsOfHitSlots or reduce /run/numberOfThreads!"
                     << std::endl;
 
-        // next step could fail because there are more tracks in flight than hit slots left..
-        bool nextStepMightFail =
-            gpuState.stats->hitBufferOccupancy >= gpuState.fHitScoring->HitCapacity() / numThreads - 1.5 * maxInFlight;
+        // next step could fail because there are too tracks in flight. A track can cause multiple hits (one for each
+        // secondary and one for itself). The safety factor should be as low as possible to prevent stalling, but must
+        // be as high as needed to avoid crashes.
+        bool nextStepMightFail = gpuState.stats->hitBufferOccupancy + hitBufferSafetyFactor * maxInFlight >=
+                                 gpuState.fHitScoring->HitCapacity() / numThreads;
 
         if (!gpuState.fHitScoring->ReadyToSwapBuffers() && !nextStepMightFail) {
           hitProcessing->cv.notify_one();
@@ -1777,7 +1780,7 @@ std::thread LaunchGPUWorker(int trackCapacity, int leakCapacity, int scoringCapa
                             TrackBuffer &trackBuffer, GPUstate &gpuState,
                             std::vector<std::atomic<EventState>> &eventStates, std::condition_variable &cvG4Workers,
                             std::vector<AdePTScoring> &scoring, int adeptSeed, int debugLevel, bool returnAllSteps,
-                            bool returnLastStep, unsigned short lastNParticlesOnCPU)
+                            bool returnLastStep, unsigned short lastNParticlesOnCPU, const double hitBufferSafetyFactor)
 {
   return std::thread{&TransportLoop,
                      trackCapacity,
@@ -1793,7 +1796,8 @@ std::thread LaunchGPUWorker(int trackCapacity, int leakCapacity, int scoringCapa
                      debugLevel,
                      returnAllSteps,
                      returnLastStep,
-                     lastNParticlesOnCPU};
+                     lastNParticlesOnCPU,
+                     hitBufferSafetyFactor};
 }
 
 void FreeGPU(std::unique_ptr<AsyncAdePT::GPUstate, AsyncAdePT::GPUstateDeleter> &gpuState, G4HepEmState &g4hepem_state,

include/AdePT/core/AsyncAdePTTransport.hh

@@ -78,6 +78,8 @@ private:
   ///< Filling fraction of the ScoringCapacity on host when the hits are copied out and not taken directly by the
   ///< G4workers
   double fCPUCopyFraction{0.5};
+  ///< Needed to stall the GPU, in case the nPartInFlight * fHitBufferSafetyFactor > available HitSlots
+  double fHitBufferSafetyFactor{1.5};
 
   void Initialize(G4HepEmTrackingManagerSpecialized *hepEmTM);
   void InitBVH();

include/AdePT/core/AsyncAdePTTransport.icc

@@ -56,7 +56,7 @@ std::pair<GPUHit *, GPUHit *> GetGPUHitsFromBuffer(unsigned int, unsigned int, A
 void CloseGPUBuffer(unsigned int, AsyncAdePT::GPUstate &, GPUHit *, const bool);
 std::thread LaunchGPUWorker(int, int, int, int, AsyncAdePT::TrackBuffer &, AsyncAdePT::GPUstate &,
                             std::vector<std::atomic<AsyncAdePT::EventState>> &, std::condition_variable &,
-                            std::vector<AdePTScoring> &, int, int, bool, bool, unsigned short);
+                            std::vector<AdePTScoring> &, int, int, bool, bool, unsigned short, const double);
 std::unique_ptr<AsyncAdePT::GPUstate, AsyncAdePT::GPUstateDeleter> InitializeGPU(
     int trackCapacity, int leakCapacity, int scoringCapacity, int numThreads, AsyncAdePT::TrackBuffer &trackBuffer,
     std::vector<AdePTScoring> &scoring, double CPUCapacityFactor, double CPUCopyFraction);
@@ -100,7 +100,8 @@ AsyncAdePTTransport<IntegrationLayer>::AsyncAdePTTransport(AdePTConfiguration &c
       fReturnAllSteps{configuration.GetCallUserSteppingAction()},
       fReturnFirstAndLastStep{configuration.GetCallUserTrackingAction() || configuration.GetCallUserSteppingAction()},
       fBfieldFile{configuration.GetCovfieBfieldFile()}, fCPUCapacityFactor{configuration.GetCPUCapacityFactor()},
-      fCPUCopyFraction{configuration.GetHitBufferFlushThreshold()}
+      fCPUCopyFraction{configuration.GetHitBufferFlushThreshold()},
+      fHitBufferSafetyFactor{configuration.GetHitBufferSafetyFactor()}
 {
   if (fNThread > kMaxThreads)
     throw std::invalid_argument("AsyncAdePTTransport limited to " + std::to_string(kMaxThreads) + " threads");
@@ -157,7 +158,7 @@ void AsyncAdePTTransport<IntegrationLayer>::AddTrack(
                          static_cast<short>(threadId)};
   if (fDebugLevel >= 2) {
     fGPUNetEnergy[threadId] += energy;
-    if (fDebugLevel >= 6) {
+    if (fDebugLevel >= 8) {
       G4cout << "\n[_in," << eventId << "," << trackId << "]: " << track << "\tGPU net energy " << std::setprecision(6)
              << fGPUNetEnergy[threadId] << G4endl;
     }
@@ -317,11 +318,12 @@ void AsyncAdePTTransport<IntegrationLayer>::Initialize(G4HepEmTrackingManagerSpe
 
   assert(fBuffer != nullptr);
 
-  fGPUstate  = async_adept_impl::InitializeGPU(fTrackCapacity, fLeakCapacity, fScoringCapacity, fNThread, *fBuffer,
-                                               fScoring, fCPUCapacityFactor, fCPUCopyFraction);
-  fGPUWorker = async_adept_impl::LaunchGPUWorker(
-      fTrackCapacity, fLeakCapacity, fScoringCapacity, fNThread, *fBuffer, *fGPUstate, fEventStates, fCV_G4Workers,
-      fScoring, fAdePTSeed, fDebugLevel, fReturnAllSteps, fReturnFirstAndLastStep, fLastNParticlesOnCPU);
+  fGPUstate = async_adept_impl::InitializeGPU(fTrackCapacity, fLeakCapacity, fScoringCapacity, fNThread, *fBuffer,
+                                              fScoring, fCPUCapacityFactor, fCPUCopyFraction);
+  fGPUWorker =
+      async_adept_impl::LaunchGPUWorker(fTrackCapacity, fLeakCapacity, fScoringCapacity, fNThread, *fBuffer, *fGPUstate,
+                                        fEventStates, fCV_G4Workers, fScoring, fAdePTSeed, fDebugLevel, fReturnAllSteps,
+                                        fReturnFirstAndLastStep, fLastNParticlesOnCPU, fHitBufferSafetyFactor);
 }
 
 template <typename IntegrationLayer>
@@ -412,7 +414,7 @@ void AsyncAdePTTransport<IntegrationLayer>::Flush(G4int threadId, G4int eventId)
     for (const auto &track : tracks) {
 
       fGPUNetEnergy[threadId] -= track.eKin;
-      if (fDebugLevel >= 5) {
+      if (fDebugLevel >= 8) {
         G4cout << "\n[out," << track.eventId << "," << trackId++ << "]: " << track << "\tGPU net energy "
                << std::setprecision(6) << fGPUNetEnergy[threadId] << G4endl;
       }

include/AdePT/integration/AdePTConfigurationMessenger.hh

@@ -56,6 +56,7 @@ private:
   G4UIcmdWithADouble *fSetMillionsOfHitSlotsCmd;
   G4UIcmdWithADouble *fSetHitBufferFlushThresholdCmd;
   G4UIcmdWithADouble *fSetCPUCapacityFactorCmd;
+  G4UIcmdWithADouble *fSetHitBufferSafetyFactorCmd;
 
   // Temporary method for setting the VecGeom geometry.
   // In the future the geometry will be converted from Geant4 rather than loaded from GDML.

src/AdePTConfigurationMessenger.cc

@@ -92,6 +92,11 @@ AdePTConfigurationMessenger::AdePTConfigurationMessenger(AdePTConfiguration *ade
   fSetCPUCapacityFactorCmd->SetParameterName("CPUCapacityFactor", false);
   fSetCPUCapacityFactorCmd->SetRange("CPUCapacityFactor>=2.5");
 
+  fSetHitBufferSafetyFactorCmd = new G4UIcmdWithADouble("/adept/setHitBufferSafetyFactor", this);
+  fSetHitBufferSafetyFactorCmd->SetGuidance(
+      "Sets the HitBuffer safety factor for stalling the GPU. If nParticlesInFlight * HitBufferSafetyFactor > "
+      "NumHitSlotsLeft, the GPU will stall. Default: 1.5 ");
+
   fSetGDMLCmd = new G4UIcmdWithAString("/adept/setVecGeomGDML", this);
   fSetGDMLCmd->SetGuidance("Temporary method for setting the geometry to use with VecGeom");
 
@@ -155,6 +160,7 @@ AdePTConfigurationMessenger::~AdePTConfigurationMessenger()
   delete fSetFinishOnCpuCmd;
   delete fSetSpeedOfLightCmd;
   delete fSetCPUCapacityFactorCmd;
+  delete fSetHitBufferSafetyFactorCmd;
   delete fSetMultipleStepsInMSCWithTransportationCmd;
   delete fSetEnergyLossFluctuationCmd;
   delete fSetMaxWDTIterCmd;
@@ -202,6 +208,8 @@ void AdePTConfigurationMessenger::SetNewValue(G4UIcommand *command, G4String new
     fAdePTConfiguration->SetHitBufferFlushThreshold(fSetHitBufferFlushThresholdCmd->GetNewDoubleValue(newValue));
   } else if (command == fSetCPUCapacityFactorCmd) {
     fAdePTConfiguration->SetCPUCapacityFactor(fSetCPUCapacityFactorCmd->GetNewDoubleValue(newValue));
+  } else if (command == fSetHitBufferSafetyFactorCmd) {
+    fAdePTConfiguration->SetHitBufferSafetyFactor(fSetHitBufferSafetyFactorCmd->GetNewDoubleValue(newValue));
   } else if (command == fSetGDMLCmd) {
     fAdePTConfiguration->SetVecGeomGDML(newValue);
   } else if (command == fSetCovfieFileCmd) {