Remove bumps and relocation-on-find from vecgeom tracker

src/celeritas/alongstep/detail/PropagationApplier.hh

@@ -114,6 +114,11 @@ PropagationApplierBaseImpl<MP>::operator()(CoreTrackView& track)
 #endif
         auto propagate = make_propagator(track);
         p = propagate(sim.step_length());
+        if (CELER_UNLIKELY(p.failed))
+        {
+            track.apply_errored();
+            return;
+        }
         tracks_can_loop = propagate.tracks_can_loop();
         CELER_ASSERT(p.distance > 0);
 #if CELERITAS_DEBUG

src/celeritas/field/FieldPropagator.hh

@@ -197,7 +197,13 @@ FieldPropagator<SubstepperT, GTV>::operator()(real_type step) -> result_type
         real_type const update_length = substep.length * linear_step.distance
                                         / chord.length;
 
-        if (!linear_step.boundary)
+        if (CELER_UNLIKELY(geo_.failed()))
+        {
+            // Time to die
+            result.failed = true;
+            return result;
+        }
+        else if (!linear_step.boundary)
         {
             // No boundary intersection along the chord: accept substep
             // movement inside the current volume and reset the remaining

src/celeritas/field/LinearPropagator.hh

@@ -6,6 +6,7 @@
 //---------------------------------------------------------------------------//
 #pragma once
 
+#include "corecel/Macros.hh"
 #include "corecel/math/Algorithms.hh"
 #include "geocel/Types.hh"
 
@@ -78,7 +79,11 @@ CELER_FUNCTION auto LinearPropagator<GTV>::operator()(real_type dist)
 
     result_type result = geo_.find_next_step(dist);
 
-    if (result.boundary)
+    if (CELER_UNLIKELY(geo_.failed()))
+    {
+        result.failed = true;
+    }
+    else if (result.boundary)
     {
         geo_.move_to_boundary();
     }

src/corecel/math/Algorithms.hh

@@ -537,6 +537,8 @@ inline CELER_FUNCTION T fastpow(T a, T b)
 /*!
  * Use fused multiply-add for generic calculations.
  *
+ * \f[ x \gets a \times b + y \f]
+ *
  * This provides a floating point specialization so that \c fma can be used in
  * code that is accelerated for floating point calculations but still works
  * correctly with integer arithmetic.

src/geocel/Types.hh

@@ -109,12 +109,17 @@ struct GeoTrackInitializer
  *
  * The boundary flag means that the geometry is step limiting, but the surface
  * crossing must be called externally.
+ *
+ * \todo Change to a status type, or have separate "propagation" and "next
+ * step" result types. Currently this is copied from the geometry state to the
+ * propagation result because the track view failure flags are ephemeral.
  */
 struct Propagation
 {
     real_type distance{0};  //!< Distance traveled
     bool boundary{false};  //!< True if hit a boundary before given distance
     bool looping{false};  //!< True if track is looping in the field propagator
+    bool failed{false};  //!< True if a failure state occurred
 };
 
 //---------------------------------------------------------------------------//

src/geocel/detail/LengthUnits.hh

@@ -23,14 +23,17 @@ namespace lengthunits
 CELER_ICC meter{100};
 CELER_ICC centimeter{1};
 CELER_ICC millimeter{0.1};
+inline constexpr char const label[] = "cm";
 #elif CELERITAS_UNITS == CELERITAS_UNITS_SI
 CELER_ICC meter{1};
 CELER_ICC centimeter{0.01};
 CELER_ICC millimeter{0.001};
+inline constexpr char const label[] = "m";
 #elif CELERITAS_UNITS == CELERITAS_UNITS_CLHEP
 CELER_ICC meter{1000};
 CELER_ICC centimeter{10};
 CELER_ICC millimeter{1};
+inline constexpr char const label[] = "mm";
 #else
 #    error "CELERITAS_UNITS is undefined"
 #endif

src/geocel/vg/VecgeomTrackView.hh

@@ -41,6 +41,12 @@
 #    include "detail/VgNavStateWrapper.hh"
 #endif
 
+#if !CELER_DEVICE_COMPILE
+#    include "corecel/io/Logger.hh"
+#    include "corecel/io/Repr.hh"
+#    include "geocel/detail/LengthUnits.hh"
+#endif
+
 namespace celeritas
 {
 //---------------------------------------------------------------------------//
@@ -54,8 +60,9 @@ namespace celeritas
     VecgeomTrackView geom(vg_params_ref, vg_state_ref, trackslot_id);
    \endcode
  *
- * The "next distance" is cached as part of `find_next_step`, but it is only
- * used when the immediate next call is `move_to_boundary`.
+ * \note VecGeom's solid model is imprcise and requires bumping. The bump
+ * values are arbitrary and may lead to geometry errors. The values of the bump
+ * constants do not account for the problem's length scale.
  */
 class VecgeomTrackView
 {
@@ -86,11 +93,6 @@ class VecgeomTrackView
     inline CELER_FUNCTION VecgeomTrackView&
     operator=(Initializer_t const& init);
 
-    //// STATIC ACCESSORS ////
-
-    //! A tiny push to make sure tracks do not get stuck at boundaries
-    static CELER_CONSTEXPR_FUNCTION real_type extra_push() { return 1e-13; }
-
     //// ACCESSORS ////
 
     //!@{
@@ -121,7 +123,7 @@ class VecgeomTrackView
     // Whether the track is exactly on a surface
     CELER_FORCEINLINE_FUNCTION bool is_on_boundary() const;
     //! Whether the last operation resulted in an error
-    CELER_FORCEINLINE_FUNCTION bool failed() const { return false; }
+    CELER_FORCEINLINE_FUNCTION bool failed() const { return failed_; }
     // Get the normal vector of the current surface
     inline CELER_FUNCTION Real3 normal() const;
 
@@ -184,6 +186,19 @@ class VecgeomTrackView
 
     // Temporary data
     real_type next_step_{0};
+    bool failed_{false};
+
+    // Static data
+
+    //! Absolute tolerance in vecgeom::kTolerance
+    //! 1e-9 for double, 1e-3 for single
+    static constexpr vg_real_type abs_tolerance_
+        = CELERITAS_REAL_TYPE == CELERITAS_REAL_TYPE_DOUBLE ? 1e-9 : 1e-3;
+
+    //! Tolerance used for solid model relocation bump
+    //! V1 solids by default;
+    static constexpr vg_real_type relocate_bump_
+        = CELERITAS_VECGEOM_SURFACE ? 0 : 10 * abs_tolerance_;
 
     //// HELPER FUNCTIONS ////
 
@@ -193,9 +208,6 @@ class VecgeomTrackView
     // Whether the next distance-to-boundary is to a surface
     inline CELER_FUNCTION bool is_next_boundary() const;
 
-    // Get a reference to the world volume
-    inline CELER_FUNCTION VgPlacedVol const& world() const;
-
     // Get a reference to the current volume instance
     inline CELER_FUNCTION VgPlacedVol const& physical_volume() const;
 
@@ -254,6 +266,7 @@ CELER_FUNCTION VecgeomTrackView&
 VecgeomTrackView::operator=(Initializer_t const& init)
 {
     CELER_EXPECT(is_soft_unit_vector(init.dir));
+    failed_ = false;
 
     // Initialize direction
     dir_ = init.dir;
@@ -285,14 +298,27 @@ VecgeomTrackView::operator=(Initializer_t const& init)
     // Set up current state and locate daughter volume
     vgstate_.Clear();
 #if CELERITAS_VECGEOM_SURFACE
-    auto world = vecgeom::NavigationState::WorldId();
+    // VecGeom's BVHSurfNav takes `int pvol_id` but vecgeom's navtuple/index
+    // return NavIndex_t via `NavInd` :(
+    VgPlacedVolumeInt world = vecgeom::NavigationState::WorldId();
 #else
-    auto const* world = &this->world();
+    auto const* world = params_.scalars.world<MemSpace::native>();
 #endif
     // LocatePointIn sets `vgstate_`
     constexpr bool contains_point = true;
     Navigator::LocatePointIn(
         world, detail::to_vector(pos_), vgstate_, contains_point);
+
+    if (CELER_UNLIKELY(vgstate_.IsOutside()))
+    {
+#if !CELER_DEVICE_COMPILE
+        auto msg = CELER_LOG_LOCAL(error);
+        msg << "Failed to initialize geometry state at " << repr(pos_)
+            << lengthunits::label;
+#endif
+        failed_ = true;
+    }
+
     return *this;
 }
 
@@ -425,32 +451,12 @@ CELER_FUNCTION Real3 VecgeomTrackView::normal() const
 /*!
  * Find the distance to the next geometric boundary.
  *
- * This function is allowed to be called from the exterior for ray tracing.
+ * \todo To support ray tracing from unknown distances from the world, we
+ * could add another special function \c find_first_step .
  */
 CELER_FUNCTION Propagation VecgeomTrackView::find_next_step()
 {
-    if (this->is_outside())
-    {
-        // SPECIAL CASE: find distance to interior from outside world volume
-        auto const& world_pv = this->world();
-        next_step_ = world_pv.DistanceToIn(detail::to_vector(pos_),
-                                           detail::to_vector(dir_),
-                                           vecgeom::kInfLength);
-
-        next_step_ = max(next_step_, this->extra_push());
-        vgnext_.Clear();
-        Propagation result;
-        result.distance = next_step_;
-        result.boundary = next_step_ < vecgeom::kInfLength;
-
-        if (result.boundary)
-        {
-            vgnext_.Push(&world_pv);
-            vgnext_.SetBoundaryState(true);
-        }
-
-        return result;
-    }
+    CELER_EXPECT(!this->is_outside());
 
     return this->find_next_step(vecgeom::kInfLength);
 }
@@ -468,6 +474,7 @@ CELER_FUNCTION Propagation VecgeomTrackView::find_next_step(real_type max_step)
     {
         *next_surf_ = null_surface();
     }
+
     // TODO: vgnext is simply copied and the boundary flag optionally set
     next_step_ = Navigator::ComputeStepAndNextVolume(detail::to_vector(pos_),
                                                      detail::to_vector(dir_),
@@ -479,15 +486,30 @@ CELER_FUNCTION Propagation VecgeomTrackView::find_next_step(real_type max_step)
                                                      *next_surf_
 #endif
     );
+
+    if (CELER_UNLIKELY(!(next_step_ > 0)))
+    {
+#if !CELER_DEVICE_COMPILE
+        auto msg = CELER_LOG_LOCAL(error);
+        msg << "Failed to find next step at " << repr(pos_) << ' '
+            << lengthunits::label << " along " << repr(dir_)
+            << ": computed step is " << repr(next_step_) << ' '
+            << lengthunits::label;
+#endif
+        failed_ = true;
+        Propagation result;
+        result.distance = next_step_;
+        result.boundary = false;
+        return result;
+    }
+
     if constexpr (CELERITAS_VECGEOM_SURFACE)
     {
         // Our accessor uses the next_surf_ state, but the temporary used for
         // vgnext_ should reflect the same result
         CELER_ASSERT((*next_surf_ != null_surface()) == vgnext_.IsOnBoundary());
     }
 
-    next_step_ = max(next_step_, this->extra_push());
-
     if (!this->is_next_boundary())
     {
         // Soft equivalence between distance and max step is because the
@@ -503,9 +525,8 @@ CELER_FUNCTION Propagation VecgeomTrackView::find_next_step(real_type max_step)
 
     CELER_ENSURE(this->has_next_step());
     CELER_ENSURE(result.distance > 0);
-    CELER_ENSURE(result.distance <= max(max_step, this->extra_push()));
-    CELER_ENSURE(result.boundary || result.distance == max_step
-                 || max_step < this->extra_push());
+    CELER_ENSURE(result.distance <= max_step);
+    CELER_ENSURE(result.boundary || result.distance == max_step);
     return result;
 }
 
@@ -531,13 +552,8 @@ CELER_FUNCTION real_type VecgeomTrackView::find_safety(real_type max_radius)
     CELER_EXPECT(!this->is_on_boundary());
     CELER_EXPECT(max_radius > 0);
 
-    real_type safety = Navigator::ComputeSafety(detail::to_vector(this->pos()),
-                                                vgstate_
-#if VECGEOM_VERSION >= 0x200000
-                                                ,
-                                                max_radius
-#endif
-    );
+    real_type safety = Navigator::ComputeSafety(
+        detail::to_vector(this->pos()), vgstate_, max_radius);
     safety = min<real_type>(safety, max_radius);
 
     // Since the reported "safety" is negative if we've moved slightly beyond
@@ -570,22 +586,47 @@ CELER_FUNCTION void VecgeomTrackView::move_to_boundary()
  */
 CELER_FUNCTION void VecgeomTrackView::cross_boundary()
 {
+    CELER_EXPECT(!this->is_outside());
     CELER_EXPECT(this->is_on_boundary());
     CELER_EXPECT(this->is_next_boundary());
 
     // Relocate to next tracking volume (maybe across multiple boundaries)
     if (vgnext_.Top() != nullptr)
     {
-        if (!CELERITAS_VECGEOM_SURFACE || !vgstate_.IsOutside())
+        VgReal3 bumped_pos;
+        if constexpr (CELERITAS_VECGEOM_SURFACE)
         {
-            // In surf model, relocation does not work from [OUTSIDE]
-            Navigator::RelocateToNextVolume(detail::to_vector(this->pos_),
-                                            detail::to_vector(this->dir_),
+            // Surface relocation is exact, assuming a well constructed
+            // geometry
+            bumped_pos = detail::to_vector(pos_);
+        }
+        else
+        {
+            for (auto i : range(3))
+            {
+                bumped_pos[i] = fma(relocate_bump_, dir_[i], pos_[i]);
+            }
+        }
+
+        Navigator::RelocateToNextVolume(bumped_pos,
+                                        detail::to_vector(this->dir_),
 #if CELERITAS_VECGEOM_SURFACE
-                                            *next_surf_,
+                                        *next_surf_,
 #endif
-                                            vgnext_);
-        }
+                                        vgnext_);
+    }
+
+    // Relocation should have changed volume
+    if (CELER_UNLIKELY(vgstate_.HasSamePathAsOther(vgnext_)))
+    {
+#if !CELER_DEVICE_COMPILE
+        auto msg = CELER_LOG_LOCAL(error);
+        msg << "Failed to cross boundary: unique volume instance is the same "
+               "before and after at "
+            << repr(pos_) << ' ' << lengthunits::label;
+#endif
+        failed_ = true;
+        return;
     }
 
     vgstate_ = vgnext_;
@@ -636,6 +677,9 @@ CELER_FUNCTION void VecgeomTrackView::move_internal(Real3 const& pos)
  *
  * This happens after a scattering event or movement inside a magnetic field.
  * It resets the calculated distance-to-boundary.
+ *
+ * \todo If on a boundary, determine as with ORANGE whether we should cancel
+ * the surface crossing.
  */
 CELER_FUNCTION void VecgeomTrackView::set_dir(Real3 const& newdir)
 {
@@ -672,17 +716,6 @@ CELER_FUNCTION bool VecgeomTrackView::is_next_boundary() const
     }
 }
 
-//---------------------------------------------------------------------------//
-/*!
- * Get a reference to the world volume instance.
- */
-CELER_FUNCTION auto VecgeomTrackView::world() const -> VgPlacedVol const&
-{
-    auto* pv = params_.scalars.world<MemSpace::native>();
-    CELER_ENSURE(pv);
-    return *pv;
-}
-
 //---------------------------------------------------------------------------//
 /*!
  * Get a reference to the current volume.