Lazily evaluate senses

src/orange/OrangeData.hh

@@ -473,6 +473,7 @@ struct OrangeStateData
 
     // Scratch space with dimensions {track}{max_faces}
     Items<Sense> temp_sense;
+    Items<SenseModFlags> temp_sense_mod;
 
     // Scratch space with dimensions {track}{max_intersections}
     Items<FaceId> temp_face;
@@ -500,6 +501,7 @@ struct OrangeStateData
             && vol.size() == max_depth  * this->size()
             && universe.size() == max_depth  * this->size()
             && !temp_sense.empty()
+            && !temp_sense_mod.empty()
             && !temp_face.empty()
             && temp_distance.size() == temp_face.size()
             && temp_isect.size() == temp_face.size();
@@ -533,6 +535,7 @@ struct OrangeStateData
         universe = other.universe;
 
         temp_sense = other.temp_sense;
+        temp_sense_mod = other.temp_sense_mod;
 
         temp_face = other.temp_face;
         temp_distance = other.temp_distance;
@@ -576,6 +579,7 @@ inline void resize(OrangeStateData<Ownership::value, M>* data,
 
     size_type face_states = params.scalars.max_faces * num_tracks;
     resize(&data->temp_sense, face_states);
+    resize(&data->temp_sense_mod, face_states);
 
     size_type isect_states = params.scalars.max_intersections * num_tracks;
     resize(&data->temp_face, isect_states);

src/orange/OrangeTrackView.hh

@@ -198,6 +198,7 @@ class OrangeTrackView
 
     // Create local sense reference
     inline CELER_FUNCTION Span<Sense> make_temp_sense() const;
+    inline CELER_FUNCTION Span<SenseModFlags> make_temp_sense_mod() const;
 
     // Create local distance
     inline CELER_FUNCTION detail::TempNextFace make_temp_next() const;
@@ -276,6 +277,7 @@ OrangeTrackView::operator=(Initializer_t const& init)
     local.volume = {};
     local.surface = {};
     local.temp_sense = this->make_temp_sense();
+    local.temp_sense_mod = this->make_temp_sense_mod();
 
     // Helpers for applying parent-to-daughter transformations
     TransformVisitor apply_transform{params_};
@@ -689,6 +691,7 @@ CELER_FUNCTION void OrangeTrackView::cross_boundary()
         local.volume = lsa.vol();
         local.surface = {this->surf(), this->sense()};
         local.temp_sense = this->make_temp_sense();
+        local.temp_sense_mod = this->make_temp_sense_mod();
     }
 
     TrackerVisitor visit_tracker{params_};
@@ -1047,6 +1050,18 @@ CELER_FUNCTION Span<Sense> OrangeTrackView::make_temp_sense() const
         offset, max_faces);
 }
 
+//---------------------------------------------------------------------------//
+/*!
+ * Get a reference to the current volume, or to world volume if outside.
+ */
+CELER_FUNCTION Span<SenseModFlags> OrangeTrackView::make_temp_sense_mod() const
+{
+    auto const max_faces = params_.scalars.max_faces;
+    auto offset = track_slot_.get() * max_faces;
+    return states_.temp_sense_mod[AllItems<SenseModFlags, MemSpace::native>{}]
+        .subspan(offset, max_faces);
+}
+
 //---------------------------------------------------------------------------//
 /*!
  * Set up intersection scratch space.
@@ -1089,6 +1104,7 @@ OrangeTrackView::make_local_state(LevelId level) const
         local.surface = {};
     }
     local.temp_sense = this->make_temp_sense();
+    local.temp_sense_mod = this->make_temp_sense_mod();
     local.temp_next = this->make_temp_next();
     return local;
 }

src/orange/OrangeTypes.hh

@@ -101,6 +101,18 @@ enum class Sense : bool
     outside,  //!< Expression is greater than zero
 };
 
+//---------------------------------------------------------------------------//
+/*!
+ * Transformations to apply to senses when using lazy sense evaluation.
+ */
+enum class SenseMod : unsigned char
+{
+    normal = 0,
+    flipped = 1 << 0,
+    cached = 1 << 1,
+};
+using SenseModFlags = std::underlying_type_t<SenseMod>;
+
 //---------------------------------------------------------------------------//
 /*!
  * Enumeration for mapping surface classes to integers.
@@ -378,6 +390,46 @@ CELER_CONSTEXPR_FUNCTION Sense to_sense(bool s)
     return static_cast<SignedSense>(-static_cast<IntT>(orig));
 }
 
+//---------------------------------------------------------------------------//
+/*!
+ * Check if a sense modifier is set.
+ */
+[[nodiscard]] CELER_CONSTEXPR_FUNCTION bool
+is_sense_mod_set(SenseMod mod, SenseModFlags flags)
+{
+    return (flags & static_cast<SenseModFlags>(mod)) != 0;
+}
+
+//---------------------------------------------------------------------------//
+/*!
+ * Set a sense modifier.
+ */
+[[nodiscard]] CELER_CONSTEXPR_FUNCTION SenseModFlags
+set_sense_mod(SenseMod mod, SenseModFlags flags)
+{
+    return flags | static_cast<SenseModFlags>(mod);
+}
+
+//---------------------------------------------------------------------------//
+/*!
+ * Unset a sense modifier.
+ */
+[[nodiscard]] CELER_CONSTEXPR_FUNCTION SenseModFlags
+unset_sense_mod(SenseMod mod, SenseModFlags flags)
+{
+    return flags & ~static_cast<SenseModFlags>(mod);
+}
+
+//---------------------------------------------------------------------------//
+/*!
+ * Flip a sense modifier.
+ */
+[[nodiscard]] CELER_CONSTEXPR_FUNCTION SenseModFlags
+flip_sense_mod(SenseMod mod, SenseModFlags flags)
+{
+    return flags ^ static_cast<SenseModFlags>(mod);
+}
+
 //---------------------------------------------------------------------------//
 /*!
  * Change whether a boundary crossing is reentrant or exiting.

src/orange/univ/SimpleUnitTracker.hh

@@ -10,10 +10,12 @@
 #include "corecel/Assert.hh"
 #include "corecel/math/Algorithms.hh"
 #include "orange/OrangeData.hh"
+#include "orange/OrangeTypes.hh"
 #include "orange/detail/BIHEnclosingVolFinder.hh"
 #include "orange/surf/LocalSurfaceVisitor.hh"
 
 #include "detail/InfixEvaluator.hh"
+#include "detail/LazySenseCalculator.hh"
 #include "detail/LogicEvaluator.hh"
 #include "detail/SenseCalculator.hh"
 #include "detail/SurfaceFunctors.hh"
@@ -165,22 +167,24 @@ SimpleUnitTracker::initialize(LocalState const& state) const -> Initialization
     CELER_EXPECT(params_);
     CELER_EXPECT(!state.surface && !state.volume);
 
-    detail::SenseCalculator calc_senses(
-        this->make_surface_visitor(), state.pos, state.temp_sense);
+    detail::LazySenseCalculator calc_senses(this->make_surface_visitor(),
+                                            state.pos,
+                                            state.temp_sense,
+                                            state.temp_sense_mod);
 
     // Use the BIH to locate a position that's inside, and save whether it's on
     // a surface in the found volume
-    bool on_surface{false};
-    auto is_inside
-        = [this, &calc_senses, &on_surface](LocalVolumeId id) -> bool {
+    auto is_inside = [this, &calc_senses](LocalVolumeId id) -> bool {
         VolumeView vol = this->make_local_volume(id);
-        auto logic_state = calc_senses(vol);
-        on_surface = static_cast<bool>(logic_state.face);
-        return detail::LogicEvaluator(vol.logic())(logic_state.senses);
+        auto bind_calc_sense
+            = [&](FaceId face_id) { return calc_senses(vol, face_id); };
+        auto inside = detail::LogicEvaluator(vol.logic())(bind_calc_sense);
+        calc_senses.invalidate_cache();
+        return inside;
     };
     LocalVolumeId id = this->find_volume_where(state.pos, is_inside);
 
-    if (on_surface)
+    if (static_cast<bool>(calc_senses.on_face()))
     {
         // Prohibit initialization on a surface
         id = {};
@@ -202,26 +206,29 @@ CELER_FUNCTION auto
 SimpleUnitTracker::cross_boundary(LocalState const& state) const -> Initialization
 {
     CELER_EXPECT(state.surface && state.volume);
-    detail::SenseCalculator calc_senses(
-        this->make_surface_visitor(), state.pos, state.temp_sense);
 
     detail::OnLocalSurface on_surface;
-    auto is_inside = [this, &state, &calc_senses, &on_surface](
-                         LocalVolumeId const& id) -> bool {
+    auto is_inside
+        = [this, &state, &on_surface](LocalVolumeId const& id) -> bool {
+        detail::LazySenseCalculator calc_senses(this->make_surface_visitor(),
+                                                state.pos,
+                                                state.temp_sense,
+                                                state.temp_sense_mod);
         if (id == state.volume)
         {
             // Cannot cross surface into the same volume
             return false;
         }
 
         VolumeView vol = this->make_local_volume(id);
-        auto logic_state
-            = calc_senses(vol, detail::find_face(vol, state.surface));
-
-        if (detail::LogicEvaluator(vol.logic())(logic_state.senses))
+        auto on_face = detail::find_face(vol, state.surface);
+        auto bind_calc_sense = [&](FaceId face_id) {
+            return calc_senses(vol, face_id, on_face);
+        };
+        if (detail::LogicEvaluator(vol.logic())(bind_calc_sense))
         {
             // Inside: find and save the local surface ID, and end the search
-            on_surface = get_surface(vol, logic_state.face);
+            on_surface = get_surface(vol, calc_senses.on_face());
             return true;
         }
         return false;
@@ -528,13 +535,16 @@ SimpleUnitTracker::complex_intersect(LocalState const& state,
     CELER_ASSERT(num_isect > 0);
 
     // Calculate local senses, taking current face into account
-    auto logic_state = detail::SenseCalculator(
-        this->make_surface_visitor(), state.pos, state.temp_sense)(
-        vol, detail::find_face(vol, state.surface));
-
+    auto calc_senses = detail::LazySenseCalculator(this->make_surface_visitor(),
+                                                   state.pos,
+                                                   state.temp_sense,
+                                                   state.temp_sense_mod);
+    auto bind_calc_sense = [&](FaceId face_id) {
+        return calc_senses(vol, face_id, detail::find_face(vol, state.surface));
+    };
     // Current senses should put us inside the volume
     detail::LogicEvaluator is_inside(vol.logic());
-    CELER_ASSERT(is_inside(logic_state.senses));
+    CELER_ASSERT(is_inside(bind_calc_sense));
 
     // Loop over distances and surface indices to cross by iterating over
     // temp_next.isect[:num_isect].
@@ -547,16 +557,17 @@ SimpleUnitTracker::complex_intersect(LocalState const& state,
         // 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))
+        calc_senses.flip_sense(face);
+        if (!is_inside(bind_calc_sense))
         {
             // 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), flip_sense(new_sense)};
+
+            result.surface
+                = {vol.get_surface(face), flip_sense(bind_calc_sense(face))};
             result.distance = state.temp_next.distance[isect];
             CELER_ENSURE(result.distance > 0 && !std::isinf(result.distance));
             return result;
@@ -620,10 +631,15 @@ CELER_FUNCTION auto SimpleUnitTracker::background_intersect(
         {
             CELER_ASSERT(vid != state.volume);
             VolumeView vol = this->make_local_volume(vid);
-            auto logic_state = detail::SenseCalculator{
-                this->make_surface_visitor(), pos, state.temp_sense}(vol);
-
-            if (detail::LogicEvaluator{vol.logic()}(logic_state.senses))
+            auto calc_senses
+                = detail::LazySenseCalculator{this->make_surface_visitor(),
+                                              pos,
+                                              state.temp_sense,
+                                              state.temp_sense_mod};
+            auto bind_calc_sense
+                = [&](FaceId face_id) { return calc_senses(vol, face_id); };
+
+            if (detail::LogicEvaluator{vol.logic()}(bind_calc_sense))
             {
                 // We are in this new volume by crossing the tested surface.
                 // Get the sense corresponding to this "crossed" surface.
@@ -633,8 +649,7 @@ CELER_FUNCTION auto SimpleUnitTracker::background_intersect(
                 Intersection result;
                 result.distance = state.temp_next.distance[isect];
                 result.surface = detail::OnLocalSurface{
-                    surface,
-                    flip_sense(logic_state.senses[face.unchecked_get()])};
+                    surface, flip_sense(bind_calc_sense(face))};
                 return result;
             }
         }