Better use of safety.

Implementing two safety features:

1. Safety caching: the tracks hold the point where safety is computed, and the safety value in that point. When safety is needed in a different point, the method track.GetSafety(new_point) will subtract from the cached value and return the positive remainder or zero otherwise. The GetSafety method can take an accurate limit value, and when the remainder is smaller than this limit the returned safety is zero. It is the user responsibility to call track.SetSafety(pos, safety) whenever the safety is fully recomputed in a new position.

2. Safety calculation usses a new surface model feature to only compute safety accurately when close to boundaries, and only use the aligned bounding box safety when far away. The distance limit is passed as third parameter to AdePTNavigator::ComputeSafety, and must be typically equal to the discrete interaction step.

Author: agheata

include/AdePT/core/Track.cuh

@@ -21,19 +21,21 @@ struct Track {
   using Precision = vecgeom::Precision;
 
   RanluxppDouble rngState;
-  double eKin;
-  double numIALeft[4];
-  double initialRange;
-  double dynamicRangeFactor;
-  double tlimitMin;
-
-  double globalTime{0};
-  double localTime{0};
-  double properTime{0};
-
-  vecgeom::Vector3D<Precision> pos;
-  vecgeom::Vector3D<Precision> dir;
-  vecgeom::NavigationState navState;
+  double eKin{0.};
+  double numIALeft[4]{0., 0., 0., 0.};
+  double initialRange{0.};
+  double dynamicRangeFactor{0.};
+  double tlimitMin{0.};
+
+  double globalTime{0.};
+  double localTime{0.};
+  double properTime{0.};
+
+  vecgeom::Vector3D<Precision> pos;   ///< track position
+  vecgeom::Vector3D<Precision> dir;   ///< track direction
+  vecgeom::Vector3D<float> safetyPos; ///< last position where the safety was computed
+  float safety{0.f};                  ///< last computed safety value
+  vecgeom::NavigationState navState;  ///< current navigation state
 
 #ifdef USE_SPLIT_KERNELS
   // Variables used to store track info needed for scoring
@@ -46,7 +48,6 @@ struct Track {
 
   // Variables used to store navigation results
   double geometryStepLength{0};
-  double safety{0};
   long hitsurfID{0};
 #endif
 
@@ -60,6 +61,29 @@ struct Track {
   bool stopped{false};
 #endif
 
+  /// @brief Get recomputed cached safety ay a given track position
+  /// @param new_pos Track position
+  /// @param accurate_limit Only return non-zero if the recomputed safety if larger than the accurate_limit
+  /// @return Recomputed safety.
+  __host__ __device__ VECGEOM_FORCE_INLINE float GetSafety(vecgeom::Vector3D<Precision> const &new_pos,
+                                                           float accurate_limit = 0.f) const
+  {
+    float dsafe = safety - accurate_limit;
+    if (dsafe <= 0.f) return 0.f;
+    float distSq = (vecgeom::Vector3D<float>(new_pos) - safetyPos).Mag2();
+    if (dsafe * dsafe < distSq) return 0.f;
+    return (safety - vecCore::math::Sqrt(distSq));
+  }
+
+  /// @brief Set Safety value computed in a new point
+  /// @param new_pos Position where the safety is computed
+  /// @param safe Safety value
+  __host__ __device__ VECGEOM_FORCE_INLINE void SetSafety(vecgeom::Vector3D<Precision> const &new_pos, float safe)
+  {
+    safetyPos.Set(static_cast<float>(new_pos[0]), static_cast<float>(new_pos[1]), static_cast<float>(new_pos[2]));
+    safety = vecCore::math::Max(safe, 0.f);
+  }
+
   __host__ __device__ double Uniform() { return rngState.Rndm(); }
 
   __host__ __device__ void InitAsSecondary(const vecgeom::Vector3D<Precision> &parentPos,
@@ -77,7 +101,9 @@ struct Track {
 
     // A secondary inherits the position of its parent; the caller is responsible
     // to update the directions.
-    this->pos      = parentPos;
+    this->pos = parentPos;
+    this->safetyPos.Set(0.f, 0.f, 0.f);
+    this->safety   = 0.0f;
     this->navState = parentNavState;
 
     // The global time is inherited from the parent

include/AdePT/kernels/electrons.cuh

@@ -110,14 +110,6 @@ static __device__ __forceinline__ void TransportElectrons(adept::TrackManager<Tr
     // used, this provides a fresh round of random numbers and reduces thread
     // divergence because the RNG state doesn't need to be advanced later.
     RanluxppDouble newRNG(currentTrack.rngState.BranchNoAdvance());
-
-    // Compute safety, needed for MSC step limit.
-    double safety = 0;
-    if (!navState.IsOnBoundary()) {
-      safety = AdePTNavigator::ComputeSafety(pos, navState);
-    }
-    theTrack->SetSafety(safety);
-
     G4HepEmRandomEngine rnge(&currentTrack.rngState);
 
     // Sample the `number-of-interaction-left` and put it into the track.
@@ -132,6 +124,24 @@ static __device__ __forceinline__ void TransportElectrons(adept::TrackManager<Tr
 
     G4HepEmElectronManager::HowFarToDiscreteInteraction(&g4HepEmData, &g4HepEmPars, &elTrack);
 
+    auto physicalStepLength = elTrack.GetPStepLength();
+    // Compute safety, needed for MSC step limit. The accuracy range is physicalStepLength
+    double safety = 0.;
+    if (!navState.IsOnBoundary()) {
+      // Get the remaining safety only if larger than physicalStepLength
+      safety = currentTrack.GetSafety(pos);
+      if (safety < physicalStepLength) {
+        // Recompute safety and update it in the track.
+#ifdef ADEPT_USE_SURF
+        // Use maximum accuracy only if safety is samller than physicalStepLength
+        safety = AdePTNavigator::ComputeSafety(pos, navState, physicalStepLength);
+#else
+        safety = AdePTNavigator::ComputeSafety(pos, navState);
+#endif
+        currentTrack.SetSafety(pos, safety);
+      }
+    }
+    theTrack->SetSafety(safety);
     bool restrictedPhysicalStepLength = false;
     if (BzFieldValue != 0) {
       const double momentumMag = sqrt(eKin * (eKin + 2.0 * restMass));
@@ -142,11 +152,11 @@ static __device__ __forceinline__ void TransportElectrons(adept::TrackManager<Tr
       double limit                = MaxSafeLength * safeLength;
       limit                       = safety > limit ? safety : limit;
 
-      double physicalStepLength = elTrack.GetPStepLength();
       if (physicalStepLength > limit) {
         physicalStepLength           = limit;
         restrictedPhysicalStepLength = true;
         elTrack.SetPStepLength(physicalStepLength);
+
         // Note: We are limiting the true step length, which is converted to
         // a shorter geometry step length in HowFarToMSC. In that sense, the
         // limit is an over-approximation, but that is fine for our purpose.
@@ -201,6 +211,7 @@ static __device__ __forceinline__ void TransportElectrons(adept::TrackManager<Tr
     // correct information (navState = nextState only if relocated
     // in case of a boundary; see below)
     navState.SetBoundaryState(nextState.IsOnBoundary());
+    if (nextState.IsOnBoundary()) currentTrack.SetSafety(pos, 0.);
 
     // Propagate information from geometrical step to MSC.
     theTrack->SetDirection(dir.x(), dir.y(), dir.z());
@@ -222,7 +233,7 @@ static __device__ __forceinline__ void TransportElectrons(adept::TrackManager<Tr
       if (dLength2 > kGeomMinLength2) {
         const double dispR = std::sqrt(dLength2);
         // Estimate safety by subtracting the geometrical step length.
-        safety -= geometryStepLength;
+        safety                 = currentTrack.GetSafety(pos);
         constexpr double sFact = 0.99;
         double reducedSafety   = sFact * safety;
 
@@ -232,7 +243,13 @@ static __device__ __forceinline__ void TransportElectrons(adept::TrackManager<Tr
           pos += displacement;
         } else {
           // Recompute safety.
-          safety        = AdePTNavigator::ComputeSafety(pos, navState);
+#ifdef ADEPT_USE_SURF
+          // Use maximum accuracy only if safety is samller than physicalStepLength
+          safety = AdePTNavigator::ComputeSafety(pos, navState, dispR);
+#else
+          safety = AdePTNavigator::ComputeSafety(pos, navState);
+#endif
+          currentTrack.SetSafety(pos, safety);
           reducedSafety = sFact * safety;
 
           // 1b. Far away from geometry boundary:

include/AdePT/magneticfield/fieldPropagatorConstBz.h

@@ -142,7 +142,12 @@ __host__ __device__ Precision fieldPropagatorConstBz::ComputeStepAndNextVolume(
     } else {
       Precision newSafety = 0;
       if (stepDone > 0) {
+#ifdef ADEPT_USE_SURF
+        // Use maximum accuracy only if safety is samller than physicalStepLength
+        newSafety = Navigator::ComputeSafety(position, current_state, remains);
+#else
         newSafety = Navigator::ComputeSafety(position, current_state);
+#endif
       }
       if (newSafety > chordLen) {
         move         = chordLen;

include/AdePT/navigation/SurfNavigator.h

@@ -49,10 +49,12 @@ class SurfNavigator {
   /// @brief Computes the isotropic safety from the globalpoint.
   /// @param globalpoint Point in global coordinates
   /// @param state Path where to compute safety
+  /// @param limit Limit to which safety should be computed accurately
   /// @return Isotropic safe distance
-  __host__ __device__ static Precision ComputeSafety(Vector3D const &globalpoint, vecgeom::NavigationState const &state)
+  __host__ __device__ static Precision ComputeSafety(Vector3D const &globalpoint, vecgeom::NavigationState const &state,
+                                                     Precision limit = vecgeom::InfinityLength<Real_t>())
   {
-    auto safety = vgbrep::protonav::BVHSurfNavigator<Real_t>::ComputeSafety(globalpoint, state);
+    auto safety = vgbrep::protonav::BVHSurfNavigator<Real_t>::ComputeSafety(globalpoint, state, limit);
     return safety;
   }