Surface navigation now split in finding the frame and relocation (#308)

The AdePT interfaces for navigation normally compute first the distance
to the next hit volume, then relocates to find out the next path. The
first interface for surface navigation did systematically
step+relocation in one go, which is very inefficient in case of non-zero
magnetic field, where the field propagator has to call the distance
computation in several points along the trajectory.

The PR introduces the split: distance computation does not change the
output state, but sets its boundary flag according to a surface being
hit within the physics step, returning also the index of this surface.
The relocation method used this index as input. Due to the addition of
this index in the interfaces, currently, the calling interface of the
surface navigator is for the moment different than the solid navigators.

---------

Co-authored-by: Severin Diederichs <65728274+SeverinDiederichs@users.noreply.github.com>
Co-authored-by: Severin Diederichs <severin.diederichs@cern.ch>
Co-authored-by: Juan González <juangonzalezcaminero@protonmail.com>

examples/Example1/include/SensitiveDetector.hh

@@ -33,6 +33,7 @@
 #include "SimpleHit.hh"
 
 #include "G4VSensitiveDetector.hh"
+#include <set>
 
 class G4HCofThisEvent;
 class G4TouchableHistory;

examples/IntegrationBenchmark/include/SensitiveDetector.hh

@@ -34,6 +34,7 @@
 
 #include "G4VFastSimSensitiveDetector.hh"
 #include "G4VSensitiveDetector.hh"
+#include <set>
 
 class G4HCofThisEvent;
 class G4TouchableHistory;

include/AdePT/core/AdePTTransport.icc

@@ -67,7 +67,7 @@ template <typename IntegrationLayer>
 bool AdePTTransport<IntegrationLayer>::InitializeGeometry(const vecgeom::cxx::VPlacedVolume *world)
 {
   auto &cudaManager = vecgeom::cxx::CudaManager::Instance();
-  bool success = true;
+  bool success      = true;
 #ifdef ADEPT_USE_SURF
 #ifdef ADEPT_USE_SURF_SINGLE
   using BrepHelper = vgbrep::BrepHelper<float>;
@@ -83,9 +83,10 @@ bool AdePTTransport<IntegrationLayer>::InitializeGeometry(const vecgeom::cxx::VP
   cudaManager.SynchronizeNavigationTable();
 #else
   // Upload solid geometry to GPU.
+  cudaDeviceSetLimit(cudaLimitStackSize, 1024 * 8);
   cudaManager.LoadGeometry(world);
   auto world_dev = cudaManager.Synchronize();
-  success = world_dev != nullptr;
+  success        = world_dev != nullptr;
 #endif
   // Initialize BVH
   InitBVH();

include/AdePT/core/HostScoringImpl.cuh

@@ -169,14 +169,14 @@ namespace adept_scoring
     aGPUHit->fStepLength         = aStepLength;
     aGPUHit->fTotalEnergyDeposit = aTotalEnergyDeposit;
     // Pre step point
-    aGPUHit->fPreStepPoint.fNavigationStateIndex = aPreState->GetNavIndex();
+    aGPUHit->fPreStepPoint.fNavigationState = *aPreState;
     Copy3DVector(aPrePosition, &(aGPUHit->fPreStepPoint.fPosition));
     Copy3DVector(aPreMomentumDirection, &(aGPUHit->fPreStepPoint.fMomentumDirection));
     // Copy3DVector(aPrePolarization, aGPUHit.fPreStepPoint.fPolarization);
     aGPUHit->fPreStepPoint.fEKin   = aPreEKin;
     aGPUHit->fPreStepPoint.fCharge = aPreCharge;
     // Post step point
-    aGPUHit->fPostStepPoint.fNavigationStateIndex = aPostState->GetNavIndex();
+    aGPUHit->fPostStepPoint.fNavigationState = *aPostState;
     Copy3DVector(aPostPosition, &(aGPUHit->fPostStepPoint.fPosition));
     Copy3DVector(aPostMomentumDirection, &(aGPUHit->fPostStepPoint.fMomentumDirection));
     // Copy3DVector(aPostPolarization, aGPUHit.fPostStepPoint.fPolarization);

include/AdePT/core/HostScoringStruct.cuh

@@ -15,7 +15,7 @@ struct GPUStepPoint {
   double fEKin;
   double fCharge;
   // Data needed to reconstruct G4 Touchable history
-  NavIndex_t fNavigationStateIndex{0}; // VecGeom navigation state index, used to identify the touchable
+  vecgeom::NavigationState fNavigationState{0}; // VecGeom navigation state, used to identify the touchable
 };
 
 // Stores the necessary data to reconstruct GPU hits on the host , and

include/AdePT/integration/AdePTGeant4Integration.hh

@@ -63,7 +63,7 @@ public:
 
 private:
   /// @brief Reconstruct G4TouchableHistory from a VecGeom Navigation index
-  void FillG4NavigationHistory(unsigned int aNavIndex, G4NavigationHistory *aG4NavigationHistory);
+  void FillG4NavigationHistory(vecgeom::NavigationState aNavState, G4NavigationHistory *aG4NavigationHistory);
 
   void FillG4Step(GPUHit *aGPUHit, G4Step *aG4Step, G4TouchableHandle &aPreG4TouchableHandle,
                   G4TouchableHandle &aPostG4TouchableHandle);

include/AdePT/kernels/electrons.cuh

@@ -170,15 +170,22 @@ static __device__ __forceinline__ void TransportElectrons(adept::TrackManager<Tr
     // also need to carry them over!
 
     // Check if there's a volume boundary in between.
-    bool propagated = true;
+    bool propagated    = true;
+    long hitsurf_index = -1;
     double geometryStepLength;
     vecgeom::NavigationState nextState;
     if (BzFieldValue != 0) {
       geometryStepLength = fieldPropagatorBz.ComputeStepAndNextVolume<AdePTNavigator>(
-          eKin, restMass, Charge, geometricalStepLengthFromPhysics, pos, dir, navState, nextState, propagated, safety);
+          eKin, restMass, Charge, geometricalStepLengthFromPhysics, pos, dir, navState, nextState, hitsurf_index,
+          propagated, safety);
     } else {
+#ifdef ADEPT_USE_SURF
+      geometryStepLength = AdePTNavigator::ComputeStepAndNextVolume(pos, dir, geometricalStepLengthFromPhysics,
+                                                                    navState, nextState, hitsurf_index, kPush);
+#else
       geometryStepLength = AdePTNavigator::ComputeStepAndNextVolume(pos, dir, geometricalStepLengthFromPhysics,
                                                                     navState, nextState, kPush);
+#endif
       pos += geometryStepLength * dir;
     }
 
@@ -244,6 +251,18 @@ static __device__ __forceinline__ void TransportElectrons(adept::TrackManager<Tr
     localTime += deltaTime;
     properTime += deltaTime * (restMass / eKin);
 
+    if (nextState.IsOnBoundary()) {
+      // if the particle hit a boundary, and is neither stopped or outside, relocate to have the correct next state
+      // before RecordHit is called
+      if (!stopped && !nextState.IsOutside()) {
+#ifdef ADEPT_USE_SURF
+        AdePTNavigator::RelocateToNextVolume(pos, dir, hitsurf_index, nextState);
+#else
+        AdePTNavigator::RelocateToNextVolume(pos, dir, nextState);
+#endif
+      }
+    }
+
     if (auxData.fSensIndex >= 0)
       adept_scoring::RecordHit(userScoring, currentTrack.parentID,
                                IsElectron ? 0 : 1,       // Particle type
@@ -306,15 +325,14 @@ static __device__ __forceinline__ void TransportElectrons(adept::TrackManager<Tr
 
       // Kill the particle if it left the world.
       if (!nextState.IsOutside()) {
-        AdePTNavigator::RelocateToNextVolume(pos, dir, nextState);
 
         // Move to the next boundary.
         navState = nextState;
         // Check if the next volume belongs to the GPU region and push it to the appropriate queue
 #ifndef ADEPT_USE_SURF
-        const int nextlvolID          = navState.Top()->GetLogicalVolume()->id();
+        const int nextlvolID = navState.Top()->GetLogicalVolume()->id();
 #else
-        const int nextlvolID          = navState.GetLogicalId();
+        const int nextlvolID = navState.GetLogicalId();
 #endif
         VolAuxData const &nextauxData = auxDataArray[nextlvolID];
         if (nextauxData.fGPUregion > 0)

include/AdePT/kernels/gammas.cuh

@@ -95,8 +95,18 @@ __global__ void TransportGammas(adept::TrackManager<Track> *gammas, Secondaries
 
     // Check if there's a volume boundary in between.
     vecgeom::NavigationState nextState;
-    double geometryStepLength = AdePTNavigator::ComputeStepAndNextVolume(pos, dir, geometricalStepLengthFromPhysics,
-                                                                         navState, nextState, kPush);
+    double geometryStepLength;
+#ifdef ADEPT_USE_SURF
+    long hitsurf_index = -1;
+    geometryStepLength = AdePTNavigator::ComputeStepAndNextVolume(pos, dir, geometricalStepLengthFromPhysics, navState,
+                                                                  nextState, hitsurf_index, kPush);
+#else
+    geometryStepLength = AdePTNavigator::ComputeStepAndNextVolume(pos, dir, geometricalStepLengthFromPhysics, navState,
+                                                                  nextState, kPush);
+#endif
+    //  printf("pvol=%d  step=%g  onboundary=%d  pos={%g, %g, %g}  dir={%g, %g, %g}\n", navState.TopId(),
+    //  geometryStepLength,
+    //         nextState.IsOnBoundary(), pos[0], pos[1], pos[2], dir[0], dir[1], dir[2]);
     pos += geometryStepLength * dir;
 
     // Set boundary state in navState so the next step and secondaries get the
@@ -121,15 +131,20 @@ __global__ void TransportGammas(adept::TrackManager<Track> *gammas, Secondaries
 
       // Kill the particle if it left the world.
       if (!nextState.IsOutside()) {
+#ifdef ADEPT_USE_SURF
+        AdePTNavigator::RelocateToNextVolume(pos, dir, hitsurf_index, nextState);
+        if (nextState.IsOutside()) continue;
+#else
         AdePTNavigator::RelocateToNextVolume(pos, dir, nextState);
-
+#endif
         // Move to the next boundary.
         navState = nextState;
-        // Check if the next volume belongs to the GPU region and push it to the appropriate queue
+        // printf("  -> pvol=%d pos={%g, %g, %g} \n", navState.TopId(), pos[0], pos[1], pos[2]);
+        //  Check if the next volume belongs to the GPU region and push it to the appropriate queue
 #ifndef ADEPT_USE_SURF
-        const int nextlvolID          = navState.Top()->GetLogicalVolume()->id();
+        const int nextlvolID = navState.Top()->GetLogicalVolume()->id();
 #else
-        const int nextlvolID          = navState.GetLogicalId();
+        const int nextlvolID = navState.GetLogicalId();
 #endif
         VolAuxData const &nextauxData = auxDataArray[nextlvolID];
         if (nextauxData.fGPUregion > 0)