The surface model uses BVH-accelerated Locate and Safety. No solids o…

…n GPU (#304)

This PR adapts the AdePT surface navigator to the index-based
interfaces. All navigation calls use the BVH surface navigator. The
solid model is no longer copied to GPU if `ADEPT_USE_SURF ` is ON.

The new interfaces seem to work correctly with non-overlapping
geometries, but we will need an example generating hits to validate
that.

CMakeLists.txt

@@ -92,16 +92,18 @@ endif()
 if(ADEPT_USE_SURF)
   if(VecGeom_SURF_FOUND)
     if (ADEPT_USE_SURF_SINGLE)
-      message(STATUS "${Green}Surface model in VecGeom uses single precision${ColorReset}")
+      message(STATUS "${Green}Using the surface model in VecGeom in single precision${ColorReset}")
       add_compile_definitions(ADEPT_USE_SURF_SINGLE)
     else()
-      message(STATUS "${Green}Surface model in VecGeom uses double precision${ColorReset}")
+      message(STATUS "${Green}Using the surface model in VecGeom in double precision${ColorReset}")
     endif()
     add_compile_definitions(ADEPT_USE_SURF)
   else()
     message(STATUS "${Magenta}No VecGeom surface support. Forcing ADEPT_USE_SURF to OFF${ColorReset}")
     set(ADEPT_USE_SURF OFF CACHE BOOL "Disable using the surface model" FORCE)
   endif()
+else()
+  message(STATUS "${Green}Using the solid model in VecGeom${ColorReset}")
 endif()
 
 # Find Geant4, optional for now

include/AdePT/core/AdePTTransport.cuh

@@ -111,8 +111,8 @@ __global__ void InitTracks(adeptint::TrackData *trackinfo, int ntracks, int star
     };
     assert(trackmgr != nullptr && "Unsupported pdg type");
 
-    Track &track    = trackmgr->NextTrack();
-    track.parentID  = trackinfo[i].parentID;
+    Track &track   = trackmgr->NextTrack();
+    track.parentID = trackinfo[i].parentID;
 
     track.rngState.SetSeed(1234567 * event + startTrack + i);
     track.eKin         = trackinfo[i].eKin;
@@ -134,12 +134,20 @@ __global__ void InitTracks(adeptint::TrackData *trackinfo, int ntracks, int star
     track.navState.Clear();
     // We locate the pushed point because we run the risk that the
     // point is not located in the GPU region
+#ifndef ADEPT_USE_SURF
     AdePTNavigator::LocatePointIn(world, track.pos + tolerance * track.dir, track.navState, true);
+#else
+    AdePTNavigator::LocatePointIn(vecgeom::NavigationState::WorldId(), track.pos + tolerance * track.dir,
+                                  track.navState, true);
+#endif
     // The track must be on boundary at this point
     track.navState.SetBoundaryState(true);
     // nextState is initialized as needed.
-    auto volume = track.navState.Top();
-    int lvolID  = volume->GetLogicalVolume()->id();
+#ifndef ADEPT_USE_SURF
+    int lvolID  = track.navState.Top()->GetLogicalVolume()->id();
+#else
+    int lvolID  = track.navState.GetLogicalId();
+#endif
     assert(auxDataArray[lvolID].fGPUregion);
   }
 }

include/AdePT/core/AdePTTransport.icc

@@ -66,6 +66,8 @@ void AdePTTransport<IntegrationLayer>::AddTrack(int pdg, int parent_id, double e
 template <typename IntegrationLayer>
 bool AdePTTransport<IntegrationLayer>::InitializeGeometry(const vecgeom::cxx::VPlacedVolume *world)
 {
+  auto &cudaManager = vecgeom::cxx::CudaManager::Instance();
+  bool success = true;
 #ifdef ADEPT_USE_SURF
 #ifdef ADEPT_USE_SURF_SINGLE
   using BrepHelper = vgbrep::BrepHelper<float>;
@@ -77,14 +79,17 @@ bool AdePTTransport<IntegrationLayer>::InitializeGeometry(const vecgeom::cxx::VP
   if (!BrepHelper::Instance().Convert()) return 1;
   BrepHelper::Instance().PrintSurfData();
   std::cout << "== Conversion to surface model done in " << timer.Stop() << " [s]\n";
-#endif
-  // Upload geometry to GPU.
-  auto &cudaManager = vecgeom::cxx::CudaManager::Instance();
+  // Upload only navigation table to the GPU
+  cudaManager.SynchronizeNavigationTable();
+#else
+  // Upload solid geometry to GPU.
   cudaManager.LoadGeometry(world);
   auto world_dev = cudaManager.Synchronize();
+  success = world_dev != nullptr;
+#endif
   // Initialize BVH
   InitBVH();
-  return (world_dev != nullptr);
+  return success;
 }
 
 template <typename IntegrationLayer>
@@ -168,7 +173,9 @@ template <typename IntegrationLayer>
 void AdePTTransport<IntegrationLayer>::InitBVH()
 {
   vecgeom::BVHManager::Init();
+#ifndef ADEPT_USE_SURF
   vecgeom::BVHManager::DeviceInit();
+#endif
 }
 
 template <typename IntegrationLayer>

include/AdePT/kernels/electrons.cuh

@@ -67,10 +67,14 @@ static __device__ __forceinline__ void TransportElectrons(adept::TrackManager<Tr
     double localTime  = currentTrack.localTime;
     double properTime = currentTrack.properTime;
     auto navState     = currentTrack.navState;
-    const auto volume = navState.Top();
     adeptint::TrackData trackdata;
     // the MCC vector is indexed by the logical volume id
-    const int lvolID          = volume->GetLogicalVolume()->id();
+#ifndef ADEPT_USE_SURF
+    const int lvolID = navState.Top()->GetLogicalVolume()->id();
+#else
+    const int lvolID = navState.GetLogicalId();
+#endif
+
     VolAuxData const &auxData = auxDataArray[lvolID];
 
     auto survive = [&](bool leak = false) {
@@ -301,14 +305,17 @@ static __device__ __forceinline__ void TransportElectrons(adept::TrackManager<Tr
       // For now, just count that we hit something.
 
       // Kill the particle if it left the world.
-      if (nextState.Top() != nullptr) {
+      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
-        const auto nextvolume         = navState.Top();
-        const int nextlvolID          = nextvolume->GetLogicalVolume()->id();
+#ifndef ADEPT_USE_SURF
+        const int nextlvolID          = navState.Top()->GetLogicalVolume()->id();
+#else
+        const int nextlvolID          = navState.GetLogicalId();
+#endif
         VolAuxData const &nextauxData = auxDataArray[nextlvolID];
         if (nextauxData.fGPUregion > 0)
           survive();

include/AdePT/kernels/gammas.cuh

@@ -45,10 +45,13 @@ __global__ void TransportGammas(adept::TrackManager<Track> *gammas, Secondaries
     double localTime  = currentTrack.localTime;
     double properTime = currentTrack.properTime;
     auto navState     = currentTrack.navState;
-    const auto volume = navState.Top();
     adeptint::TrackData trackdata;
     // the MCC vector is indexed by the logical volume id
-    int lvolID                = volume->GetLogicalVolume()->id();
+#ifndef ADEPT_USE_SURF
+    int lvolID = navState.Top()->GetLogicalVolume()->id();
+#else
+    int lvolID = navState.GetLogicalId();
+#endif
     VolAuxData const &auxData = auxDataArray[lvolID];
 
     auto survive = [&](bool leak = false) {
@@ -117,14 +120,17 @@ __global__ void TransportGammas(adept::TrackManager<Track> *gammas, Secondaries
       // For now, just count that we hit something.
 
       // Kill the particle if it left the world.
-      if (nextState.Top() != nullptr) {
+      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
-        const auto nextvolume         = navState.Top();
-        const int nextlvolID          = nextvolume->GetLogicalVolume()->id();
+#ifndef ADEPT_USE_SURF
+        const int nextlvolID          = navState.Top()->GetLogicalVolume()->id();
+#else
+        const int nextlvolID          = navState.GetLogicalId();
+#endif
         VolAuxData const &nextauxData = auxDataArray[nextlvolID];
         if (nextauxData.fGPUregion > 0)
           survive();

include/AdePT/navigation/SurfNavigator.h

@@ -26,27 +26,34 @@ template <typename Real_t>
 class SurfNavigator {
 
 public:
-  using Precision          = vecgeom::Precision;
-  using Vector3D           = vecgeom::Vector3D<vecgeom::Precision>;
-  using VPlacedVolumePtr_t = vecgeom::VPlacedVolume const *;
-  using SurfData           = vgbrep::SurfData<Real_t>;
-  using Real_b             = typename SurfData::Real_b;
+  using Precision = vecgeom::Precision;
+  using Vector3D  = vecgeom::Vector3D<vecgeom::Precision>;
+  using SurfData  = vgbrep::SurfData<Real_t>;
+  using Real_b    = typename SurfData::Real_b;
 
   static constexpr Precision kBoundaryPush = 10 * vecgeom::kTolerance;
 
-  // Locates the point in the geometry volume tree
-  __host__ __device__ static VPlacedVolumePtr_t LocatePointIn(VPlacedVolumePtr_t vol, Vector3D const &point,
-                                                              vecgeom::NavigationState &path, bool top,
-                                                              VPlacedVolumePtr_t exclude = nullptr)
+  /// @brief Locates the point in the geometry volume tree
+  /// @param pvol_id Placed volume id to be checked first
+  /// @param point Point to be checked, in the local frame of pvol
+  /// @param path Path to a parent of pvol that must contain the point
+  /// @param top Check first if pvol contains the point
+  /// @param exclude Placed volume id to exclude from the search
+  /// @return Index of the placed volume that contains the point
+  __host__ __device__ static int LocatePointIn(int pvol_id, Vector3D const &point, vecgeom::NavigationState &path,
+                                               bool top, int *exclude = nullptr)
   {
-    return vgbrep::protonav::LocatePointIn(vol, point, path, top);
+    return vgbrep::protonav::BVHSurfNavigator<Real_t>::LocatePointIn(pvol_id, point, path, top, exclude);
   }
 
-  // Computes the isotropic safety from the globalpoint.
+  /// @brief Computes the isotropic safety from the globalpoint.
+  /// @param globalpoint Point in global coordinates
+  /// @param state Path where to compute safety
+  /// @return Isotropic safe distance
   __host__ __device__ static Precision ComputeSafety(Vector3D const &globalpoint, vecgeom::NavigationState const &state)
   {
-    int closest_surf = 0;
-    return vgbrep::protonav::ComputeSafety(globalpoint, state, closest_surf);
+    auto safety = vgbrep::protonav::BVHSurfNavigator<Real_t>::ComputeSafety(globalpoint, state);
+    return safety;
   }
 
   // Computes a step from the globalpoint (which must be in the current volume)
@@ -67,8 +74,8 @@ class SurfNavigator {
     }
 
     vgbrep::CrossedSurface crossed_surf;
-    auto step =
-        vgbrep::protonav::BVHSurfNavigator<Real_t>::ComputeStepAndHit(globalpoint, globaldir, in_state, out_state, crossed_surf, step_limit);
+    auto step = vgbrep::protonav::BVHSurfNavigator<Real_t>::ComputeStepAndHit(globalpoint, globaldir, in_state,
+                                                                              out_state, crossed_surf, step_limit);
     return step;
   }
 

test/testField/electrons.cu

@@ -7,8 +7,6 @@
 
 #include <AdePT/copcore/PhysicalConstants.h>
 
-#define NOFLUCTUATION
-
 // Classes for Runge-Kutta integration
 #include <AdePT/magneticfield/MagneticFieldEquation.h>
 #include <AdePT/magneticfield/DormandPrinceRK45.h>
@@ -90,10 +88,14 @@ static __device__ __forceinline__ void TransportElectrons(Track *electrons, cons
     auto pos            = currentTrack.pos;
     auto dir            = currentTrack.dir;
     auto navState       = currentTrack.navState;
-    const auto volume   = navState.Top();
-    const int volumeID  = volume->id();
+#ifndef ADEPT_USE_SURF
+    const int volumeID  = navState.Top()->GetLogicalVolume()->id();
+#else
+    const int volumeID  = navState.GetLogicalId();
+#endif
+
     // the MCC vector is indexed by the logical volume id
-    const int theMCIndex = MCIndex[volume->GetLogicalVolume()->id()];
+    const int theMCIndex = MCIndex[volumeID];
 
     auto survive = [&] {
       currentTrack.energy   = energy;
@@ -213,8 +215,8 @@ static __device__ __forceinline__ void TransportElectrons(Track *electrons, cons
       constexpr Precision thresholdDiff = 3.0e-3;
       bool diffLength = false, badPosition = false, badDirection = false;
       vecgeom::NavigationState &currNavState = navState;
-      bool sameLevel                       = nextState.GetLevel() == nextStateHx.GetLevel();
-      bool sameIndex                       = nextState.GetNavIndex() == nextStateHx.GetNavIndex();
+      bool sameLevel                         = nextState.GetLevel() == nextStateHx.GetLevel();
+      bool sameIndex                         = nextState.GetNavIndex() == nextStateHx.GetNavIndex();
 
       if (std::fabs(helixStepLength - geometryStepLength) > 1.0e-4 * helixStepLength) {
         bool sameNextVol = (nextState.GetLevel() == nextStateHx.GetLevel()) &&
@@ -247,8 +249,8 @@ static __device__ __forceinline__ void TransportElectrons(Track *electrons, cons
       }
 #endif
     } else {
-      geometryStepLength = AdePTNavigator::ComputeStepAndNextVolume(pos, dir, geometricalStepLengthFromPhysics, navState,
-                                                                  nextState, kPush);
+      geometryStepLength = AdePTNavigator::ComputeStepAndNextVolume(pos, dir, geometricalStepLengthFromPhysics,
+                                                                    navState, nextState, kPush);
       pos += geometryStepLength * dir;
     }
 
@@ -353,7 +355,7 @@ static __device__ __forceinline__ void TransportElectrons(Track *electrons, cons
       atomicAdd(&globalScoring->hits, 1);
 
       // Kill the particle if it left the world.
-      if (nextState.Top() != nullptr) {
+      if (!nextState.IsOutside()) {
         AdePTNavigator::RelocateToNextVolume(pos, dir, nextState);
 
         // Move to the next boundary.

test/testField/gammas.cu

@@ -36,10 +36,13 @@ __global__ void TransportGammas(Track *gammas, const adept::MParray *active, Sec
     auto pos            = currentTrack.pos;
     auto dir            = currentTrack.dir;
     auto navState       = currentTrack.navState;
-    const auto volume   = navState.Top();
-    const int volumeID  = volume->id();
+#ifndef ADEPT_USE_SURF
+    const int volumeID  = navState.Top()->GetLogicalVolume()->id();
+#else
+    const int volumeID  = navState.GetLogicalId();
+#endif
     // the MCC vector is indexed by the logical volume id
-    const int theMCIndex = MCIndex[volume->GetLogicalVolume()->id()];
+    const int theMCIndex = MCIndex[volumeID];
 
     auto survive = [&] {
       currentTrack.energy   = energy;
@@ -102,7 +105,7 @@ __global__ void TransportGammas(Track *gammas, const adept::MParray *active, Sec
       atomicAdd(&globalScoring->hits, 1);
 
       // Kill the particle if it left the world.
-      if (nextState.Top() != nullptr) {
+      if (!nextState.IsOutside()) {
         AdePTNavigator::RelocateToNextVolume(pos, dir, nextState);
 
         // Move to the next boundary.

test/testField/testField.cpp

@@ -28,6 +28,9 @@
 #include <VecGeom/management/BVHManager.h>
 #include <VecGeom/gdml/Frontend.h>
 #include <VecGeom/management/CudaManager.h>
+#ifdef ADEPT_USE_SURF
+#include <VecGeom/surfaces/BrepHelper.h>
+#endif
 
 #include <AdePT/copcore/SystemOfUnits.h>
 #include <AdePT/base/ArgParser.h>
@@ -171,7 +174,9 @@ void FreeG4HepEm(G4HepEmState *state)
 void InitBVH()
 {
   vecgeom::cxx::BVHManager::Init();
+#ifndef ADEPT_USE_SURF
   vecgeom::cxx::BVHManager::DeviceInit();
+#endif
 }
 
 void PrintScoringPerVolume(const vecgeom::VPlacedVolume *placed, const ScoringPerVolume *scoring, int level = 0)
@@ -243,8 +248,22 @@ int main(int argc, char *argv[])
   // Load and synchronize the geometry on the GPU
   std::cout << "synchronizing VecGeom geometry to GPU ...\n";
   auto &cudaManager = vecgeom::cxx::CudaManager::Instance();
+#ifdef ADEPT_USE_SURF
+#ifdef ADEPT_USE_SURF_SINGLE
+  using BrepHelper = vgbrep::BrepHelper<float>;
+#else
+  using BrepHelper = vgbrep::BrepHelper<double>;
+#endif
+#endif
+
+#ifndef ADEPT_USE_SURF
   cudaManager.LoadGeometry(world);
   cudaManager.Synchronize();
+#else
+  if (!BrepHelper::Instance().Convert()) return 1;
+  BrepHelper::Instance().PrintSurfData();
+  cudaManager.SynchronizeNavigationTable();
+#endif
   InitBVH();
 
   auto time_cpu = timer.Stop();

test/testField/testField.cu

@@ -132,7 +132,11 @@ __global__ void InitPrimaries(ParticleGenerator generator, int startEvent, int n
       track.dir = {1.0, 0, 0};
     }
     track.navState.Clear();
+#ifndef ADEPT_USE_SURF
     AdePTNavigator::LocatePointIn(world, track.pos, track.navState, true);
+#else
+    AdePTNavigator::LocatePointIn(vecgeom::NavigationState::WorldId(), track.pos, track.navState, true);
+#endif
 
     atomicAdd(&globalScoring->numElectrons, 1);
   }