diff --git a/CMakeLists.txt b/CMakeLists.txt index 648f7722..1a4dbc28 100644 --- a/CMakeLists.txt +++ b/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 diff --git a/include/AdePT/core/AdePTTransport.cuh b/include/AdePT/core/AdePTTransport.cuh index d4469c01..15a665d0 100644 --- a/include/AdePT/core/AdePTTransport.cuh +++ b/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); } } diff --git a/include/AdePT/core/AdePTTransport.icc b/include/AdePT/core/AdePTTransport.icc index 6b207401..4eba5cbb 100644 --- a/include/AdePT/core/AdePTTransport.icc +++ b/include/AdePT/core/AdePTTransport.icc @@ -66,6 +66,8 @@ void AdePTTransport::AddTrack(int pdg, int parent_id, double e template bool AdePTTransport::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; @@ -77,14 +79,17 @@ bool AdePTTransport::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 @@ -168,7 +173,9 @@ template void AdePTTransport::InitBVH() { vecgeom::BVHManager::Init(); +#ifndef ADEPT_USE_SURF vecgeom::BVHManager::DeviceInit(); +#endif } template diff --git a/include/AdePT/kernels/electrons.cuh b/include/AdePT/kernels/electrons.cuh index 00369993..38bd20f9 100644 --- a/include/AdePT/kernels/electrons.cuh +++ b/include/AdePT/kernels/electrons.cuh @@ -67,10 +67,14 @@ static __device__ __forceinline__ void TransportElectrons(adept::TrackManagerGetLogicalVolume()->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::TrackManagerGetLogicalVolume()->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(); diff --git a/include/AdePT/kernels/gammas.cuh b/include/AdePT/kernels/gammas.cuh index 5983e72c..5250267b 100644 --- a/include/AdePT/kernels/gammas.cuh +++ b/include/AdePT/kernels/gammas.cuh @@ -45,10 +45,13 @@ __global__ void TransportGammas(adept::TrackManager *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 *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(); diff --git a/include/AdePT/navigation/SurfNavigator.h b/include/AdePT/navigation/SurfNavigator.h index 35257c24..90b60219 100644 --- a/include/AdePT/navigation/SurfNavigator.h +++ b/include/AdePT/navigation/SurfNavigator.h @@ -26,27 +26,34 @@ template class SurfNavigator { public: - using Precision = vecgeom::Precision; - using Vector3D = vecgeom::Vector3D; - using VPlacedVolumePtr_t = vecgeom::VPlacedVolume const *; - using SurfData = vgbrep::SurfData; - using Real_b = typename SurfData::Real_b; + using Precision = vecgeom::Precision; + using Vector3D = vecgeom::Vector3D; + using SurfData = vgbrep::SurfData; + 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::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::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::ComputeStepAndHit(globalpoint, globaldir, in_state, out_state, crossed_surf, step_limit); + auto step = vgbrep::protonav::BVHSurfNavigator::ComputeStepAndHit(globalpoint, globaldir, in_state, + out_state, crossed_surf, step_limit); return step; } diff --git a/test/testField/electrons.cu b/test/testField/electrons.cu index 5a4be50f..c0d3fa5e 100644 --- a/test/testField/electrons.cu +++ b/test/testField/electrons.cu @@ -7,8 +7,6 @@ #include -#define NOFLUCTUATION - // Classes for Runge-Kutta integration #include #include @@ -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. diff --git a/test/testField/gammas.cu b/test/testField/gammas.cu index 522d7793..c5d43235 100644 --- a/test/testField/gammas.cu +++ b/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. diff --git a/test/testField/testField.cpp b/test/testField/testField.cpp index 2deb7edd..c36fd1bc 100644 --- a/test/testField/testField.cpp +++ b/test/testField/testField.cpp @@ -28,6 +28,9 @@ #include #include #include +#ifdef ADEPT_USE_SURF +#include +#endif #include #include @@ -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; +#else + using BrepHelper = vgbrep::BrepHelper; +#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(); diff --git a/test/testField/testField.cu b/test/testField/testField.cu index 3fa7303e..0d4dbe82 100644 --- a/test/testField/testField.cu +++ b/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); }