ORANGE: implement distance-to-boundary calculation

src/orange/Data.hh

@@ -206,7 +206,10 @@ struct OrangeStateData
     StateItems<Sense>     sense;
 
     // Scratch space
-    Items<Sense> temp_senses; // [track][max_faces]
+    Items<Sense>     temp_sense;    // [track][max_faces]
+    Items<FaceId>    temp_face;     // [track][max_intersections]
+    Items<real_type> temp_distance; // [track][max_intersections]
+    Items<size_type> temp_isect;    // [track][max_intersections]
 
     //// METHODS ////
 
@@ -219,7 +222,10 @@ struct OrangeStateData
             && vol.size() == pos.size()
             && surf.size() == pos.size()
             && sense.size() == pos.size()
-            && !temp_senses.empty();
+            && !temp_sense.empty()
+            && !temp_face.empty()
+            && temp_distance.size() == temp_face.size()
+            && temp_isect.size() == temp_face.size();
         // clang-format on
     }
 
@@ -238,7 +244,10 @@ struct OrangeStateData
         surf  = other.surf;
         sense = other.sense;
 
-        temp_senses = other.temp_senses;
+        temp_sense    = other.temp_sense;
+        temp_face     = other.temp_face;
+        temp_distance = other.temp_distance;
+        temp_isect    = other.temp_isect;
 
         CELER_ENSURE(*this);
         return *this;
@@ -264,7 +273,12 @@ resize(OrangeStateData<Ownership::value, M>* data,
     make_builder(&data->sense).resize(size);
 
     size_type face_states = params.scalars.max_faces * size;
-    make_builder(&data->temp_senses).resize(face_states);
+    make_builder(&data->temp_sense).resize(face_states);
+
+    size_type isect_states = params.scalars.max_intersections * size;
+    make_builder(&data->temp_face).resize(isect_states);
+    make_builder(&data->temp_distance).resize(isect_states);
+    make_builder(&data->temp_isect).resize(isect_states);
 
     CELER_ENSURE(*data);
 }

src/orange/universes/SimpleUnitTracker.hh

@@ -11,6 +11,7 @@
 #include "orange/surfaces/Surfaces.hh"
 #include "detail/LogicEvaluator.hh"
 #include "detail/SenseCalculator.hh"
+#include "detail/SurfaceFunctors.hh"
 #include "detail/Types.hh"
 #include "detail/Utils.hh"
 
@@ -33,6 +34,7 @@ class SimpleUnitTracker
     using ParamsRef
         = OrangeParamsData<Ownership::const_reference, MemSpace::native>;
     using Initialization = detail::Initialization;
+    using Intersection   = detail::Intersection;
     using LocalState     = detail::LocalState;
     //!@}
 
@@ -43,8 +45,18 @@ class SimpleUnitTracker
     // Find the local cell and possibly surface ID.
     inline CELER_FUNCTION Initialization initialize(LocalState state) const;
 
+    // Calculate the distance to an exiting face for the current volume.
+    inline CELER_FUNCTION Intersection intersect(LocalState state) const;
+
   private:
+    //// DATA ////
     const ParamsRef& params_;
+
+    //// METHODS ////
+    inline CELER_FUNCTION
+        Intersection simple_intersect(LocalState, VolumeView) const;
+    inline CELER_FUNCTION
+        Intersection complex_intersect(LocalState, VolumeView) const;
 };
 
 //---------------------------------------------------------------------------//
@@ -76,7 +88,7 @@ CELER_FUNCTION auto SimpleUnitTracker::initialize(LocalState state) const
     CELER_EXPECT(params_);
 
     detail::SenseCalculator calc_senses(
-        Surfaces{params_.surfaces}, state.pos, state.temp_senses);
+        Surfaces{params_.surfaces}, state.pos, state.temp_sense);
 
     for (VolumeId volid : range(VolumeId{params_.volumes.size()}))
     {
@@ -113,5 +125,182 @@ CELER_FUNCTION auto SimpleUnitTracker::initialize(LocalState state) const
     return {};
 }
 
+//---------------------------------------------------------------------------//
+/*!
+ * Calculate distance-to-intercept for the next surface.
+ */
+CELER_FUNCTION auto SimpleUnitTracker::intersect(LocalState state) const
+    -> Intersection
+{
+    CELER_EXPECT(state.volume && !state.temp_sense.empty());
+
+    // Resize temporaries based on volume properties
+    VolumeView vol{params_.volumes, state.volume};
+    CELER_ASSERT(state.temp_next.size >= vol.num_intersections());
+    state.temp_next.size = vol.num_intersections();
+    const bool is_simple = !(vol.flags() & VolumeView::internal_surfaces);
+
+    // Find all surface intersection distances inside this volume
+    {
+        auto calc_intersections = make_surface_action(
+            Surfaces{params_.surfaces},
+            detail::CalcIntersections{
+                state.pos,
+                state.dir,
+                state.surface ? vol.find_face(state.surface.id()) : FaceId{},
+                is_simple,
+                state.temp_next});
+        for (SurfaceId surface : vol.faces())
+        {
+            calc_intersections(surface);
+        }
+        CELER_ASSERT(calc_intersections.action().face_idx() == vol.num_faces());
+        CELER_ASSERT(calc_intersections.action().isect_idx()
+                     == vol.num_intersections());
+    }
+
+    if (is_simple)
+    {
+        // No interior surfaces: closest distance is next boundary
+        return this->simple_intersect(state, vol);
+    }
+    else
+    {
+        return this->complex_intersect(state, vol);
+    }
+}
+
+//---------------------------------------------------------------------------//
+/*!
+ * Calculate distance to the next boundary for nonreentrant cells.
+ */
+CELER_FUNCTION auto
+SimpleUnitTracker::simple_intersect(LocalState state, VolumeView vol) const
+    -> Intersection
+{
+    CELER_EXPECT(state.temp_next && vol.num_intersections() > 0);
+
+    // Crossing any surface will leave the cell; perform a linear search for
+    // the smallest (but positive) distance
+    size_type distance_idx;
+    {
+        const real_type* distance_ptr = celeritas::min_element(
+            state.temp_next.distance,
+            state.temp_next.distance + state.temp_next.size,
+            detail::CloserPositiveDistance{});
+        CELER_ASSERT(*distance_ptr > 0);
+        distance_idx = distance_ptr - state.temp_next.distance;
+    }
+
+    if (state.temp_next.distance[distance_idx] == no_intersection())
+        return {};
+
+    // Determine the crossing surface
+    SurfaceId surface;
+    {
+        FaceId face = state.temp_next.face[distance_idx];
+        CELER_ASSERT(face);
+        surface = vol.get_surface(face);
+        CELER_ASSERT(surface);
+    }
+
+    Sense cur_sense;
+    if (surface == state.surface.id())
+    {
+        // Crossing the same surface that we're currently on and inside (e.g.
+        // on the inside surface of a sphere, and the next intersection is the
+        // other side)
+        cur_sense = state.surface.sense();
+    }
+    else
+    {
+        auto calc_sense = make_surface_action(Surfaces{params_.surfaces},
+                                              detail::CalcSense{state.pos});
+
+        SignedSense ss = calc_sense(surface);
+        CELER_ASSERT(ss != SignedSense::on);
+        cur_sense = to_sense(ss);
+    }
+
+    // Post-surface sense will be on the other side of the surface
+    Intersection result;
+    result.surface  = {surface, flip_sense(cur_sense)};
+    result.distance = state.temp_next.distance[distance_idx];
+    return result;
+}
+
+//---------------------------------------------------------------------------//
+/*!
+ * Calculate boundary distance if internal surfaces are present.
+ *
+ * In "complex" cells, crossing a surface can still leave the particle in an
+ * "inside" state.
+ *
+ * We have to iteratively track through all surfaces, in order of minimum
+ * distance, to determine whether crossing them in sequence will cause us to
+ * exit the cell.
+ */
+CELER_FUNCTION auto
+SimpleUnitTracker::complex_intersect(LocalState state, VolumeView vol) const
+    -> Intersection
+{
+    // Partition intersections (enumerated from 0 as the `idx` array) into
+    // valid (finite positive) and invalid (infinite-or-negative) groups.
+    size_type num_isect = celeritas::partition(
+                              state.temp_next.isect,
+                              state.temp_next.isect + state.temp_next.size,
+                              detail::IntersectionPartitioner{state.temp_next})
+                          - state.temp_next.isect;
+    CELER_ASSERT(num_isect <= state.temp_next.size);
+
+    // Sort these finite distances in ascending order
+    celeritas::sort(state.temp_next.isect,
+                    state.temp_next.isect + num_isect,
+                    [&state](size_type a, size_type b) {
+                        return state.temp_next.distance[a]
+                               < state.temp_next.distance[b];
+                    });
+
+    // Calculate local senses, taking current face into account
+    auto logic_state = detail::SenseCalculator(
+        Surfaces{params_.surfaces}, state.pos, state.temp_sense)(
+        vol, detail::find_face(vol, state.surface));
+
+    // Current senses should put us inside the cell
+    detail::LogicEvaluator is_inside(vol.logic());
+    CELER_ASSERT(is_inside(logic_state.senses));
+
+    // Loop over distances and surface indices to cross by iterating over
+    // temp_next.isect[:num_isect].
+    // Evaluate the logic expression at each crossing to determine whether
+    // we're actually leaving the cell.
+    for (size_type isect_idx = 0; isect_idx != num_isect; ++isect_idx)
+    {
+        // Index into the distance/face arrays
+        const size_type isect = state.temp_next.isect[isect_idx];
+        // Face being crossed in this ordered intersection
+        FaceId face = state.temp_next.face[isect];
+        // Flip the sense of the face being crossed
+        Sense new_sense = flip_sense(logic_state.senses[face.get()]);
+        logic_state.senses[face.unchecked_get()] = new_sense;
+        if (!is_inside(logic_state.senses))
+        {
+            // Flipping this sense puts us outside the current volume: in
+            // other words, only after crossing all the internal surfaces along
+            // this direction do we hit a surface that actually puts us
+            // outside.
+            Intersection result;
+            result.surface  = {vol.get_surface(face), new_sense};
+            result.distance = state.temp_next.distance[isect];
+            CELER_ENSURE(result.distance > 0 && !std::isinf(result.distance));
+            return result;
+        }
+    }
+
+    // No intersection: perhaps leaving an exterior cell? Perhaps geometry
+    // error.
+    return {};
+}
+
 //---------------------------------------------------------------------------//
 } // namespace celeritas

src/orange/universes/detail/SurfaceFunctors.hh

@@ -102,26 +102,22 @@ struct CalcSafetyDistance
  */
 class CalcIntersections
 {
-  public:
-    //!@{
-    //! Types
-    using face_int = FaceId::size_type;
-    //!@}
-
   public:
     //! Construct from the particle point, direction, face ID, and temp storage
-    CELER_FUNCTION CalcIntersections(const Real3& pos,
-                                     const Real3& dir,
-                                     FaceId       on_face,
-                                     TempNextFace next_face)
+    CELER_FUNCTION CalcIntersections(const Real3&        pos,
+                                     const Real3&        dir,
+                                     FaceId              on_face,
+                                     bool                is_simple,
+                                     const TempNextFace& next_face)
         : pos_(pos)
         , dir_(dir)
         , on_face_idx_(on_face.unchecked_get())
-        , face_idx_(0)
+        , fill_isect_(!is_simple)
         , face_(next_face.face)
-        , dist_(next_face.distance)
+        , distance_(next_face.distance)
+        , isect_(next_face.isect)
     {
-        CELER_EXPECT(face_ && dist_);
+        CELER_EXPECT(face_ && distance_);
     }
 
     //! Operate on a surface
@@ -138,24 +134,33 @@ class CalcIntersections
         for (real_type dist : all_dist)
         {
             CELER_ASSERT(dist >= 0);
-            *face_++ = FaceId{face_idx_};
-            *dist_++ = dist;
+            face_[isect_idx_]     = FaceId{face_idx_};
+            distance_[isect_idx_] = dist;
+            if (fill_isect_)
+            {
+                isect_[isect_idx_] = isect_idx_;
+            }
+            ++isect_idx_;
         }
         // Increment to next face
         ++face_idx_;
     }
 
-    CELER_FUNCTION face_int face_idx() const { return face_idx_; }
+    CELER_FUNCTION size_type face_idx() const { return face_idx_; }
+    CELER_FUNCTION size_type isect_idx() const { return isect_idx_; }
 
   private:
     //// DATA ////
 
-    const Real3&   pos_;
-    const Real3&   dir_;
-    const face_int on_face_idx_;
-    face_int       face_idx_;
-    FaceId*        face_;
-    real_type*     dist_;
+    const Real3&     pos_;
+    const Real3&     dir_;
+    const size_type  on_face_idx_;
+    const bool       fill_isect_;
+    FaceId* const    face_;
+    real_type* const distance_;
+    size_type* const isect_;
+    size_type        face_idx_{0};
+    size_type        isect_idx_{0};
 };
 
 //---------------------------------------------------------------------------//