Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-01-18 10:05:59

0001 //----------------------------------*-C++-*----------------------------------//
0002 // Copyright 2021-2024 UT-Battelle, LLC, and other Celeritas developers.
0003 // See the top-level COPYRIGHT file for details.
0004 // SPDX-License-Identifier: (Apache-2.0 OR MIT)
0005 //---------------------------------------------------------------------------//
0006 //! \file orange/OrangeTrackView.hh
0007 //---------------------------------------------------------------------------//
0008 #pragma once
0009 
0010 #include "corecel/Macros.hh"
0011 #include "corecel/Types.hh"
0012 #include "corecel/cont/Array.hh"
0013 #include "corecel/sys/ThreadId.hh"
0014 
0015 #include "OrangeData.hh"
0016 #include "OrangeTypes.hh"
0017 #include "transform/TransformVisitor.hh"
0018 #include "univ/SimpleUnitTracker.hh"
0019 #include "univ/TrackerVisitor.hh"
0020 #include "univ/UniverseTypeTraits.hh"
0021 #include "univ/detail/Types.hh"
0022 
0023 #include "detail/LevelStateAccessor.hh"
0024 #include "detail/UniverseIndexer.hh"
0025 
0026 #if !CELER_DEVICE_COMPILE
0027 #    include "corecel/io/Logger.hh"
0028 #    include "corecel/io/Repr.hh"
0029 #endif
0030 
0031 namespace celeritas
0032 {
0033 //---------------------------------------------------------------------------//
0034 /*!
0035  * Navigate through an ORANGE geometry on a single thread.
0036  *
0037  * Ordering requirements:
0038  * - initialize (through assignment) must come first
0039  * - access (pos, dir, volume/surface/is_outside/is_on_boundary) good at any
0040  * time
0041  * - \c find_next_step
0042  * - \c find_safety or \c move_internal or \c move_to_boundary
0043  * - if on boundary, \c cross_boundary
0044  * - at any time, \c set_dir , but then must do \c find_next_step before any
0045  *   \c move or \c cross action above
0046  *
0047  * The main point is that \c find_next_step depends on the current
0048  * straight-line direction, \c move_to_boundary and \c move_internal (with
0049  * a step length) depends on that distance, and
0050  * \c cross_boundary depends on being on the boundary with a knowledge of the
0051  * post-boundary state.
0052  *
0053  * \c move_internal with a position \em should depend on the safety distance
0054  * but that's not yet implemented.
0055  */
0056 class OrangeTrackView
0057 {
0058   public:
0059     //!@{
0060     //! \name Type aliases
0061     using ParamsRef = NativeCRef<OrangeParamsData>;
0062     using StateRef = NativeRef<OrangeStateData>;
0063     using Initializer_t = GeoTrackInitializer;
0064     //!@}
0065 
0066     //! Helper struct for initializing from an existing geometry state
0067     struct DetailedInitializer
0068     {
0069         OrangeTrackView const& other;  //!< Existing geometry
0070         Real3 const& dir;  //!< New direction
0071     };
0072 
0073   public:
0074     // Construct from params and state
0075     inline CELER_FUNCTION OrangeTrackView(ParamsRef const& params,
0076                                           StateRef const& states,
0077                                           TrackSlotId tid);
0078 
0079     // Initialize the state
0080     inline CELER_FUNCTION OrangeTrackView& operator=(Initializer_t const& init);
0081     // Initialize the state from a parent state and new direction
0082     inline CELER_FUNCTION OrangeTrackView&
0083     operator=(DetailedInitializer const& init);
0084 
0085     //// ACCESSORS ////
0086 
0087     // The current position
0088     inline CELER_FUNCTION Real3 const& pos() const;
0089     // The current direction
0090     inline CELER_FUNCTION Real3 const& dir() const;
0091     // The current volume ID (null if outside)
0092     inline CELER_FUNCTION VolumeId volume_id() const;
0093     // The current surface ID
0094     inline CELER_FUNCTION SurfaceId surface_id() const;
0095     // After 'find_next_step', the next straight-line surface
0096     inline CELER_FUNCTION SurfaceId next_surface_id() const;
0097     // Whether the track is outside the valid geometry region
0098     inline CELER_FUNCTION bool is_outside() const;
0099     // Whether the track is exactly on a surface
0100     inline CELER_FUNCTION bool is_on_boundary() const;
0101     //! Whether the last operation resulted in an error
0102     CELER_FORCEINLINE_FUNCTION bool failed() const { return failed_; }
0103 
0104     //// OPERATIONS ////
0105 
0106     // Find the distance to the next boundary
0107     inline CELER_FUNCTION Propagation find_next_step();
0108 
0109     // Find the distance to the next boundary, up to and including a step
0110     inline CELER_FUNCTION Propagation find_next_step(real_type max_step);
0111 
0112     // Find the distance to the nearest boundary in any direction
0113     inline CELER_FUNCTION real_type find_safety();
0114 
0115     // Find the distance to the nearest nearby boundary in any direction
0116     inline CELER_FUNCTION real_type find_safety(real_type max_step);
0117 
0118     // Move to the boundary in preparation for crossing it
0119     inline CELER_FUNCTION void move_to_boundary();
0120 
0121     // Move within the volume
0122     inline CELER_FUNCTION void move_internal(real_type step);
0123 
0124     // Move within the volume to a specific point
0125     inline CELER_FUNCTION void move_internal(Real3 const& pos);
0126 
0127     // Cross from one side of the current surface to the other
0128     inline CELER_FUNCTION void cross_boundary();
0129 
0130     // Change direction
0131     inline CELER_FUNCTION void set_dir(Real3 const& newdir);
0132 
0133   private:
0134     //// TYPES ////
0135 
0136     using LSA = detail::LevelStateAccessor;
0137 
0138     //// DATA ////
0139 
0140     ParamsRef const& params_;
0141     StateRef const& states_;
0142     TrackSlotId track_slot_;
0143     bool failed_{false};
0144 
0145     //// PRIVATE STATE MUTATORS ////
0146 
0147     // The current level
0148     inline CELER_FUNCTION void level(LevelId);
0149 
0150     // The boundary on the current surface level
0151     inline CELER_FUNCTION void boundary(BoundaryResult);
0152 
0153     // The next step distance, as stored on the state
0154     inline CELER_FUNCTION void next_step(real_type dist);
0155 
0156     // The next surface to be encounted
0157     inline CELER_FUNCTION void next_surf(detail::OnLocalSurface const&);
0158 
0159     // The level of the next surface to be encounted
0160     inline CELER_FUNCTION void next_surface_level(LevelId);
0161 
0162     //// PRIVATE STATE ACCESSORS ////
0163 
0164     // The current level
0165     inline CELER_FUNCTION LevelId const& level() const;
0166 
0167     // The current surface level
0168     inline CELER_FUNCTION LevelId const& surface_level() const;
0169 
0170     // The local surface on the current surface level
0171     inline CELER_FUNCTION LocalSurfaceId const& surf() const;
0172 
0173     // The sense on the current surface level
0174     inline CELER_FUNCTION Sense const& sense() const;
0175 
0176     // The boundary on the current surface level
0177     inline CELER_FUNCTION BoundaryResult const& boundary() const;
0178 
0179     // The next step distance, as stored on the state
0180     inline CELER_FUNCTION real_type const& next_step() const;
0181 
0182     // The next surface to be encounted
0183     inline CELER_FUNCTION detail::OnLocalSurface next_surf() const;
0184 
0185     // The level of the next surface to be encounted
0186     inline CELER_FUNCTION LevelId const& next_surface_level() const;
0187 
0188     //// HELPER FUNCTIONS ////
0189 
0190     // Iterate over lower levels to find the next step
0191     inline CELER_FUNCTION Propagation
0192     find_next_step_impl(detail::Intersection isect);
0193 
0194     // Create local sense reference
0195     inline CELER_FUNCTION Span<Sense> make_temp_sense() const;
0196 
0197     // Create local distance
0198     inline CELER_FUNCTION detail::TempNextFace make_temp_next() const;
0199 
0200     inline CELER_FUNCTION detail::LocalState
0201     make_local_state(LevelId level) const;
0202 
0203     // Whether the next distance-to-boundary has been found
0204     inline CELER_FUNCTION bool has_next_step() const;
0205 
0206     // Whether the next surface has been found
0207     inline CELER_FUNCTION bool has_next_surface() const;
0208 
0209     // Invalidate the next distance-to-boundary and surface
0210     inline CELER_FUNCTION void clear_next();
0211 
0212     // Assign the surface on the current level
0213     inline CELER_FUNCTION void
0214     surface(LevelId level, detail::OnLocalSurface surf);
0215 
0216     // Clear the surface on the current level
0217     inline CELER_FUNCTION void clear_surface();
0218 
0219     // Make a LevelStateAccessor for the current thread and level
0220     inline CELER_FUNCTION LSA make_lsa() const;
0221 
0222     // Make a LevelStateAccessor for the current thread and a given level
0223     inline CELER_FUNCTION LSA make_lsa(LevelId level) const;
0224 
0225     // Get the daughter ID for the volume in the universe (or null)
0226     inline CELER_FUNCTION DaughterId get_daughter(LSA const& lsa);
0227 
0228     // Get the transform ID for the given daughter.
0229     inline CELER_FUNCTION TransformId get_transform(DaughterId daughter_id);
0230 
0231     // Get the transform ID to move from this level to the one below.
0232     inline CELER_FUNCTION TransformId get_transform(LevelId lev);
0233 };
0234 
0235 //---------------------------------------------------------------------------//
0236 // MEMBER FUNCTIONS
0237 //---------------------------------------------------------------------------//
0238 /*!
0239  * Construct from persistent and state data.
0240  */
0241 CELER_FUNCTION
0242 OrangeTrackView::OrangeTrackView(ParamsRef const& params,
0243                                  StateRef const& states,
0244                                  TrackSlotId tid)
0245     : params_(params), states_(states), track_slot_(tid)
0246 {
0247     CELER_EXPECT(params_);
0248     CELER_EXPECT(states_);
0249     CELER_EXPECT(track_slot_ < states.size());
0250 }
0251 
0252 //---------------------------------------------------------------------------//
0253 /*!
0254  * Construct the state.
0255  *
0256  * Expensive. This function should only be called to initialize an event from a
0257  * starting location and direction. Secondaries will initialize their states
0258  * from a copy of the parent.
0259  */
0260 CELER_FUNCTION OrangeTrackView&
0261 OrangeTrackView::operator=(Initializer_t const& init)
0262 {
0263     CELER_EXPECT(is_soft_unit_vector(init.dir));
0264 
0265     failed_ = false;
0266 
0267     // Create local state
0268     detail::LocalState local;
0269     local.pos = init.pos;
0270     local.dir = init.dir;
0271     local.volume = {};
0272     local.surface = {};
0273     local.temp_sense = this->make_temp_sense();
0274 
0275     // Helpers for applying parent-to-daughter transformations
0276     TransformVisitor apply_transform{params_};
0277     auto transform_down_local = [&local](auto&& t) {
0278         local.pos = t.transform_down(local.pos);
0279         local.dir = t.rotate_down(local.dir);
0280     };
0281 
0282     // Recurse into daughter universes starting with the outermost universe
0283     UniverseId uid = top_universe_id();
0284     DaughterId daughter_id;
0285     LevelId level{0};
0286     do
0287     {
0288         TrackerVisitor visit_tracker{params_};
0289         auto tinit = visit_tracker(
0290             [&local](auto&& t) { return t.initialize(local); }, uid);
0291 
0292         if (!tinit.volume || tinit.surface)
0293         {
0294 #if !CELER_DEVICE_COMPILE
0295             auto msg = CELER_LOG_LOCAL(error);
0296             msg << "Failed to initialize geometry state: ";
0297             if (!tinit.volume)
0298             {
0299                 msg << "could not find associated volume";
0300             }
0301             else
0302             {
0303                 msg << "started on a surface ("
0304                     << tinit.surface.id().unchecked_get() << ")";
0305             }
0306             msg << " in universe " << uid.unchecked_get()
0307                 << " at local position " << repr(local.pos);
0308 #endif
0309             // Mark as failed and place in local "exterior" to end the search
0310             // but preserve the current level information
0311             failed_ = true;
0312             tinit.volume = orange_exterior_volume;
0313         }
0314 
0315         auto lsa = this->make_lsa(level);
0316         lsa.vol() = tinit.volume;
0317         lsa.pos() = local.pos;
0318         lsa.dir() = local.dir;
0319         lsa.universe() = uid;
0320 
0321         daughter_id = visit_tracker(
0322             [&tinit](auto&& t) { return t.daughter(tinit.volume); }, uid);
0323 
0324         if (daughter_id)
0325         {
0326             auto const& daughter = params_.daughters[daughter_id];
0327             // Apply "transform down" based on stored transform
0328             apply_transform(transform_down_local, daughter.transform_id);
0329             // Update universe and increase level depth
0330             uid = daughter.universe_id;
0331             ++level;
0332         }
0333 
0334     } while (daughter_id);
0335 
0336     // Save found level
0337     this->level(level);
0338 
0339     // Reset surface/boundary information
0340     this->boundary(BoundaryResult::exiting);
0341     this->clear_surface();
0342     this->clear_next();
0343 
0344     CELER_ENSURE(!this->has_next_step());
0345     return *this;
0346 }
0347 
0348 //---------------------------------------------------------------------------//
0349 /*!
0350  * Construct the state from a direction and a copy of the parent state.
0351  */
0352 CELER_FUNCTION
0353 OrangeTrackView& OrangeTrackView::operator=(DetailedInitializer const& init)
0354 {
0355     CELER_EXPECT(is_soft_unit_vector(init.dir));
0356 
0357     failed_ = false;
0358 
0359     if (this != &init.other)
0360     {
0361         // Copy init track's position and logical state
0362         this->level(states_.level[init.other.track_slot_]);
0363         this->surface(init.other.surface_level(),
0364                       {init.other.surf(), init.other.sense()});
0365         this->boundary(init.other.boundary());
0366 
0367         for (auto lev : range(LevelId{this->level() + 1}))
0368         {
0369             // Copy all data accessed via LSA
0370             auto lsa = this->make_lsa(lev);
0371             lsa = init.other.make_lsa(lev);
0372         }
0373     }
0374 
0375     // Clear the next step information since we're changing direction or
0376     // initializing a new state
0377     this->clear_next();
0378 
0379     // Transform direction from global to local
0380     Real3 localdir = init.dir;
0381     auto apply_transform = TransformVisitor{params_};
0382     auto rotate_down
0383         = [&localdir](auto&& t) { localdir = t.rotate_down(localdir); };
0384     for (auto lev : range(LevelId{this->level()}))
0385     {
0386         auto lsa = this->make_lsa(lev);
0387         lsa.dir() = localdir;
0388         apply_transform(rotate_down,
0389                         this->get_transform(this->get_daughter(lsa)));
0390     }
0391 
0392     // Save final level
0393     this->make_lsa().dir() = localdir;
0394 
0395     CELER_ENSURE(!this->has_next_step());
0396     return *this;
0397 }
0398 
0399 //---------------------------------------------------------------------------//
0400 /*!
0401  * The current position.
0402  */
0403 CELER_FUNCTION Real3 const& OrangeTrackView::pos() const
0404 {
0405     return this->make_lsa(LevelId{0}).pos();
0406 }
0407 
0408 //---------------------------------------------------------------------------//
0409 /*!
0410  * The current direction.
0411  */
0412 CELER_FUNCTION Real3 const& OrangeTrackView::dir() const
0413 {
0414     return this->make_lsa(LevelId{0}).dir();
0415 }
0416 
0417 //---------------------------------------------------------------------------//
0418 /*!
0419  * The current volume ID.
0420  *
0421  * \note It is allowable to call this function when "outside", because the
0422  * outside in ORANGE is just a special volume. Other geometries may not have
0423  * that behavior.
0424  */
0425 CELER_FUNCTION VolumeId OrangeTrackView::volume_id() const
0426 {
0427     auto lsa = this->make_lsa();
0428     detail::UniverseIndexer ui(params_.universe_indexer_data);
0429     return ui.global_volume(lsa.universe(), lsa.vol());
0430 }
0431 
0432 //---------------------------------------------------------------------------//
0433 /*!
0434  * The current surface ID.
0435  */
0436 CELER_FUNCTION SurfaceId OrangeTrackView::surface_id() const
0437 {
0438     if (this->is_on_boundary())
0439     {
0440         auto lsa = this->make_lsa(this->surface_level());
0441         detail::UniverseIndexer ui{params_.universe_indexer_data};
0442         return ui.global_surface(lsa.universe(), this->surf());
0443     }
0444     else
0445     {
0446         return SurfaceId{};
0447     }
0448 }
0449 
0450 //---------------------------------------------------------------------------//
0451 /*!
0452  * After 'find_next_step', the next straight-line surface.
0453  */
0454 CELER_FUNCTION SurfaceId OrangeTrackView::next_surface_id() const
0455 {
0456     CELER_EXPECT(this->has_next_surface());
0457     auto lsa = this->make_lsa(this->next_surface_level());
0458     detail::UniverseIndexer ui{params_.universe_indexer_data};
0459     return ui.global_surface(lsa.universe(), this->next_surf().id());
0460 }
0461 
0462 //---------------------------------------------------------------------------//
0463 /*!
0464  * Whether the track is outside the valid geometry region.
0465  */
0466 CELER_FUNCTION bool OrangeTrackView::is_outside() const
0467 {
0468     // Zeroth volume in outermost universe is always the exterior by
0469     // construction in ORANGE
0470     auto lsa = this->make_lsa(LevelId{0});
0471     return lsa.vol() == orange_exterior_volume;
0472 }
0473 
0474 //---------------------------------------------------------------------------//
0475 /*!
0476  * Whether the track is exactly on a surface.
0477  */
0478 CELER_FORCEINLINE_FUNCTION bool OrangeTrackView::is_on_boundary() const
0479 {
0480     return static_cast<bool>(this->surface_level());
0481 }
0482 
0483 //---------------------------------------------------------------------------//
0484 /*!
0485  * Find the distance to the next geometric boundary.
0486  */
0487 CELER_FUNCTION Propagation OrangeTrackView::find_next_step()
0488 {
0489     if (CELER_UNLIKELY(this->boundary() == BoundaryResult::reentrant))
0490     {
0491         // On a boundary, headed back in: next step is zero
0492         return {0, true};
0493     }
0494 
0495     // Find intersection at the top level: always the first simple unit
0496     auto global_isect = [this] {
0497         SimpleUnitTracker t{params_, SimpleUnitId{0}};
0498         return t.intersect(this->make_local_state(LevelId{0}));
0499     }();
0500     // Find intersection for all lower levels
0501     return this->find_next_step_impl(global_isect);
0502 }
0503 
0504 //---------------------------------------------------------------------------//
0505 /*!
0506  * Find a nearby distance to the next geometric boundary up to a distance.
0507  *
0508  * This may reduce the number of surfaces needed to check, sort, or write to
0509  * temporary memory, thereby speeding up transport.
0510  */
0511 CELER_FUNCTION Propagation OrangeTrackView::find_next_step(real_type max_step)
0512 {
0513     CELER_EXPECT(max_step > 0);
0514 
0515     if (CELER_UNLIKELY(this->boundary() == BoundaryResult::reentrant))
0516     {
0517         // On a boundary, headed back in: next step is zero
0518         return {0, true};
0519     }
0520 
0521     // Find intersection at the top level: always the first simple unit
0522     auto global_isect = [this, &max_step] {
0523         SimpleUnitTracker t{params_, SimpleUnitId{0}};
0524         return t.intersect(this->make_local_state(LevelId{0}), max_step);
0525     }();
0526     // Find intersection for all lower levels
0527     auto result = this->find_next_step_impl(global_isect);
0528     CELER_ENSURE(result.distance <= max_step);
0529     return result;
0530 }
0531 
0532 //---------------------------------------------------------------------------//
0533 /*!
0534  * Move to the next straight-line boundary but do not change volume.
0535  */
0536 CELER_FUNCTION void OrangeTrackView::move_to_boundary()
0537 {
0538     CELER_EXPECT(this->boundary() != BoundaryResult::reentrant);
0539     CELER_EXPECT(this->has_next_step());
0540     CELER_EXPECT(this->has_next_surface());
0541 
0542     // Physically move next step
0543     real_type const dist = this->next_step();
0544     for (auto i : range(this->level() + 1))
0545     {
0546         auto lsa = this->make_lsa(LevelId{i});
0547         axpy(dist, lsa.dir(), &lsa.pos());
0548     }
0549 
0550     this->surface(this->next_surface_level(), this->next_surf());
0551     this->clear_next();
0552     CELER_ENSURE(this->is_on_boundary());
0553 }
0554 
0555 //---------------------------------------------------------------------------//
0556 /*!
0557  * Move within the current volume.
0558  *
0559  * The straight-line distance *must* be less than the distance to the
0560  * boundary.
0561  */
0562 CELER_FUNCTION void OrangeTrackView::move_internal(real_type dist)
0563 {
0564     CELER_EXPECT(this->has_next_step());
0565     CELER_EXPECT(dist > 0 && dist <= this->next_step());
0566     CELER_EXPECT(dist != this->next_step() || !this->has_next_surface());
0567 
0568     // Move and update the next step
0569     for (auto i : range(this->level() + 1))
0570     {
0571         auto lsa = this->make_lsa(LevelId{i});
0572         axpy(dist, lsa.dir(), &lsa.pos());
0573     }
0574     this->next_step(this->next_step() - dist);
0575     this->clear_surface();
0576 }
0577 
0578 //---------------------------------------------------------------------------//
0579 /*!
0580  * Move within the current volume to a nearby point.
0581  *
0582  * \todo Currently it's up to the caller to make sure that the position is
0583  * "nearby". We should actually test this with an "is inside" call.
0584  */
0585 CELER_FUNCTION void OrangeTrackView::move_internal(Real3 const& pos)
0586 {
0587     // Transform all nonlocal levels
0588     auto local_pos = pos;
0589     auto apply_transform = TransformVisitor{params_};
0590     auto translate_down
0591         = [&local_pos](auto&& t) { local_pos = t.transform_down(local_pos); };
0592     for (auto lev : range(LevelId{this->level()}))
0593     {
0594         auto lsa = this->make_lsa(lev);
0595         lsa.pos() = local_pos;
0596 
0597         // Apply "transform down" based on stored transform
0598         apply_transform(translate_down,
0599                         this->get_transform(this->get_daughter(lsa)));
0600     }
0601 
0602     // Save final level
0603     auto lsa = this->make_lsa();
0604     lsa.pos() = local_pos;
0605 
0606     // Clear surface state and next-step info
0607     this->clear_surface();
0608     this->clear_next();
0609 }
0610 
0611 //---------------------------------------------------------------------------//
0612 /*!
0613  * Cross from one side of the current surface to the other.
0614  *
0615  * The position *must* be on the boundary following a move-to-boundary. This
0616  * should only be called once per boundary crossing.
0617  */
0618 CELER_FUNCTION void OrangeTrackView::cross_boundary()
0619 {
0620     CELER_EXPECT(this->is_on_boundary());
0621     CELER_EXPECT(!this->has_next_step());
0622 
0623     if (CELER_UNLIKELY(this->boundary() == BoundaryResult::reentrant))
0624     {
0625         // Direction changed while on boundary leading to no change in
0626         // volume/surface. This is logically equivalent to a reflection.
0627         this->boundary(BoundaryResult::exiting);
0628         return;
0629     }
0630 
0631     // Cross surface by flipping the sense
0632     states_.sense[track_slot_] = flip_sense(this->sense());
0633     this->boundary(BoundaryResult::exiting);
0634 
0635     // Create local state from post-crossing level and updated sense
0636     LevelId level{this->surface_level()};
0637     UniverseId universe;
0638     LocalVolumeId volume;
0639     detail::LocalState local;
0640     {
0641         auto lsa = this->make_lsa(level);
0642         universe = lsa.universe();
0643         local.pos = lsa.pos();
0644         local.dir = lsa.dir();
0645         local.volume = lsa.vol();
0646         local.surface = {this->surf(), this->sense()};
0647         local.temp_sense = this->make_temp_sense();
0648     }
0649 
0650     TrackerVisitor visit_tracker{params_};
0651 
0652     // Update the post-crossing volume by crossing the boundary of the "surface
0653     // crossing" level universe
0654     volume = visit_tracker(
0655         [&local](auto&& t) { return t.cross_boundary(local).volume; },
0656         universe);
0657     if (CELER_UNLIKELY(!volume))
0658     {
0659         // Boundary crossing failure
0660 #if !CELER_DEVICE_COMPILE
0661         CELER_LOG_LOCAL(error)
0662             << "track failed to cross local surface "
0663             << this->surf().unchecked_get() << " in universe "
0664             << universe.unchecked_get() << " at local position "
0665             << repr(local.pos) << " along local direction " << repr(local.dir);
0666 #endif
0667         // Mark as failed and place in local "exterior" to end the search
0668         // but preserve the current level information
0669         failed_ = true;
0670         volume = orange_exterior_volume;
0671     }
0672     make_lsa(level).vol() = volume;
0673 
0674     // Clear local surface before diving into daughters
0675     // TODO: this is where we'd do inter-universe mapping
0676     local.volume = {};
0677     local.surface = {};
0678 
0679     // Starting with the current level (i.e., next_surface_level), iterate
0680     // down into the deepest level: *initializing* not *crossing*
0681     auto daughter_id = visit_tracker(
0682         [volume](auto&& t) { return t.daughter(volume); }, universe);
0683     while (daughter_id)
0684     {
0685         ++level;
0686         {
0687             // Update universe, local position/direction
0688             auto const& daughter = params_.daughters[daughter_id];
0689             TransformVisitor apply_transform{params_};
0690             auto transform_down_local = [&local](auto&& t) {
0691                 local.pos = t.transform_down(local.pos);
0692                 local.dir = t.rotate_down(local.dir);
0693             };
0694             apply_transform(transform_down_local, daughter.transform_id);
0695             universe = daughter.universe_id;
0696         }
0697 
0698         // Initialize in daughter and get IDs of volume and potential daughter
0699         volume = visit_tracker(
0700             [&local](auto&& t) { return t.initialize(local).volume; },
0701             universe);
0702 
0703         if (!volume)
0704         {
0705 #if !CELER_DEVICE_COMPILE
0706             auto msg = CELER_LOG_LOCAL(error);
0707             msg << "track failed to cross boundary: could not find associated "
0708                    "volume in universe "
0709                 << universe.unchecked_get() << " at local position "
0710                 << repr(local.pos);
0711 #endif
0712             // Mark as failed and place in local "exterior" to end the search
0713             // but preserve the current level information
0714             failed_ = true;
0715             volume = orange_exterior_volume;
0716         }
0717         daughter_id = visit_tracker(
0718             [volume](auto&& t) { return t.daughter(volume); }, universe);
0719 
0720         auto lsa = make_lsa(level);
0721         lsa.vol() = volume;
0722         lsa.pos() = local.pos;
0723         lsa.dir() = local.dir;
0724         lsa.universe() = universe;
0725     }
0726 
0727     // Save final level
0728     this->level(level);
0729 
0730     CELER_ENSURE(this->is_on_boundary());
0731 }
0732 
0733 //---------------------------------------------------------------------------//
0734 /*!
0735  * Change the track's direction.
0736  *
0737  * This happens after a scattering event or movement inside a magnetic field.
0738  * It resets the calculated distance-to-boundary. It is allowed to happen on
0739  * the boundary, but changing direction so that it goes from pointing outward
0740  * to inward (or vice versa) will mean that \c cross_boundary will be a
0741  * null-op.
0742  *
0743  * TODO: This needs to be updated to handle reflections through levels
0744  */
0745 CELER_FUNCTION void OrangeTrackView::set_dir(Real3 const& newdir)
0746 {
0747     CELER_EXPECT(is_soft_unit_vector(newdir));
0748 
0749     if (this->is_on_boundary())
0750     {
0751         // Changing direction on a boundary, which may result in not leaving
0752         // current volume upon the cross_surface call
0753         auto lsa = this->make_lsa(this->surface_level());
0754 
0755         TrackerVisitor visit_tracker{params_};
0756         auto normal = visit_tracker(
0757             [pos = lsa.pos(), local_surface = this->surf()](auto&& t) {
0758                 return t.normal(pos, local_surface);
0759             },
0760             lsa.universe());
0761 
0762         // Normal is in *local* coordinates but newdir is in *global*: rotate
0763         // up to check
0764         auto apply_transform = TransformVisitor{params_};
0765         auto rotate_up = [&normal](auto&& t) { normal = t.rotate_up(normal); };
0766         for (auto level : range<int>(this->level().unchecked_get()).step(-1))
0767         {
0768             apply_transform(rotate_up, this->get_transform(LevelId(level)));
0769         }
0770 
0771         // Evaluate whether the direction dotted with the surface normal
0772         // changes (i.e. heading from inside to outside or vice versa).
0773         if ((dot_product(normal, newdir) >= 0)
0774             != (dot_product(normal, this->dir()) >= 0))
0775         {
0776             // The boundary crossing direction has changed! Reverse our
0777             // plans to change the logical state and move to a new volume.
0778             this->boundary(flip_boundary(this->boundary()));
0779         }
0780     }
0781 
0782     // Complete direction setting by transforming direction all the way down
0783     Real3 localdir = newdir;
0784     auto apply_transform = TransformVisitor{params_};
0785     auto rotate_down
0786         = [&localdir](auto&& t) { localdir = t.rotate_down(localdir); };
0787     for (auto lev : range(this->level()))
0788     {
0789         auto lsa = this->make_lsa(lev);
0790         lsa.dir() = localdir;
0791         apply_transform(rotate_down,
0792                         this->get_transform(this->get_daughter(lsa)));
0793     }
0794     // Save final level
0795     this->make_lsa().dir() = localdir;
0796 
0797     this->clear_next();
0798 }
0799 
0800 //---------------------------------------------------------------------------//
0801 // STATE ACCESSORS
0802 //---------------------------------------------------------------------------//
0803 /*!
0804  * The current level.
0805  */
0806 CELER_FORCEINLINE_FUNCTION void OrangeTrackView::level(LevelId lev)
0807 {
0808     states_.level[track_slot_] = lev;
0809 }
0810 
0811 /*!
0812  * The boundary on the current surface level.
0813  */
0814 CELER_FORCEINLINE_FUNCTION void OrangeTrackView::boundary(BoundaryResult br)
0815 {
0816     states_.boundary[track_slot_] = br;
0817 }
0818 
0819 /*!
0820  * The next step distance.
0821  */
0822 CELER_FORCEINLINE_FUNCTION void OrangeTrackView::next_step(real_type dist)
0823 {
0824     states_.next_step[track_slot_] = dist;
0825 }
0826 
0827 /*!
0828  * The next surface to be encountered.
0829  */
0830 CELER_FORCEINLINE_FUNCTION void
0831 OrangeTrackView::next_surf(detail::OnLocalSurface const& s)
0832 {
0833     states_.next_surf[track_slot_] = s.id();
0834     states_.next_sense[track_slot_] = s.unchecked_sense();
0835 }
0836 
0837 /*!
0838  * The level of the next surface to be encounted.
0839  */
0840 CELER_FORCEINLINE_FUNCTION void OrangeTrackView::next_surface_level(LevelId lev)
0841 {
0842     states_.next_level[track_slot_] = lev;
0843 }
0844 
0845 //---------------------------------------------------------------------------//
0846 // CONST STATE ACCESSORS
0847 //---------------------------------------------------------------------------//
0848 /*!
0849  * The current level.
0850  */
0851 CELER_FORCEINLINE_FUNCTION LevelId const& OrangeTrackView::level() const
0852 {
0853     return states_.level[track_slot_];
0854 }
0855 
0856 /*!
0857  * The current surface level.
0858  */
0859 CELER_FORCEINLINE_FUNCTION LevelId const& OrangeTrackView::surface_level() const
0860 {
0861     return states_.surface_level[track_slot_];
0862 }
0863 
0864 /*!
0865  * The local surface on the current surface level.
0866  */
0867 CELER_FORCEINLINE_FUNCTION LocalSurfaceId const& OrangeTrackView::surf() const
0868 {
0869     return states_.surf[track_slot_];
0870 }
0871 
0872 /*!
0873  * The sense on the current surface level.
0874  */
0875 CELER_FORCEINLINE_FUNCTION Sense const& OrangeTrackView::sense() const
0876 {
0877     return states_.sense[track_slot_];
0878 }
0879 
0880 /*!
0881  * The boundary on the current surface level.
0882  */
0883 CELER_FORCEINLINE_FUNCTION BoundaryResult const&
0884 OrangeTrackView::boundary() const
0885 {
0886     return states_.boundary[track_slot_];
0887 }
0888 
0889 /*!
0890  * The next step distance.
0891  */
0892 CELER_FORCEINLINE_FUNCTION real_type const& OrangeTrackView::next_step() const
0893 {
0894     return states_.next_step[track_slot_];
0895 }
0896 
0897 /*!
0898  * The next surface to be encountered.
0899  */
0900 CELER_FORCEINLINE_FUNCTION detail::OnLocalSurface
0901 OrangeTrackView::next_surf() const
0902 {
0903     return {states_.next_surf[track_slot_], states_.next_sense[track_slot_]};
0904 }
0905 
0906 /*!
0907  * The level of the next surface to be encounted.
0908  */
0909 CELER_FORCEINLINE_FUNCTION LevelId const&
0910 OrangeTrackView::next_surface_level() const
0911 {
0912     return states_.next_level[track_slot_];
0913 }
0914 
0915 //---------------------------------------------------------------------------//
0916 // HELPER FUNCTIONS
0917 //---------------------------------------------------------------------------//
0918 /*!
0919  * Iterate over levels 1 to N to find the next step.
0920  *
0921  * Caller is responsible for finding the candidate next step on level 0, and
0922  * passing the resultant Intersection object as an argument.
0923  */
0924 CELER_FUNCTION Propagation
0925 OrangeTrackView::find_next_step_impl(detail::Intersection isect)
0926 {
0927     TrackerVisitor visit_tracker{params_};
0928 
0929     // The level with minimum distance to intersection
0930     LevelId min_level{0};
0931 
0932     // Find the nearest intersection from level 0 to current level
0933     // inclusive, prefering the shallowest level (i.e., lowest uid)
0934     for (auto levelid : range(LevelId{1}, this->level() + 1))
0935     {
0936         auto uid = this->make_lsa(levelid).universe();
0937         auto local_isect = visit_tracker(
0938             [local_state = this->make_local_state(levelid), &isect](auto&& t) {
0939                 return t.intersect(local_state, isect.distance);
0940             },
0941             uid);
0942 
0943         if (local_isect.distance < isect.distance)
0944         {
0945             isect = local_isect;
0946             min_level = levelid;
0947         }
0948     }
0949 
0950     this->next_step(isect.distance);
0951     this->next_surf(isect.surface);
0952     if (isect)
0953     {
0954         // Save level corresponding to the intersection
0955         this->next_surface_level(min_level);
0956     }
0957 
0958     Propagation result;
0959     result.distance = isect.distance;
0960     result.boundary = static_cast<bool>(isect);
0961     return result;
0962 }
0963 
0964 //---------------------------------------------------------------------------//
0965 /*!
0966  * Find the distance to the nearest boundary in any direction.
0967  *
0968  * The safety distance at a given point is the minimum safety distance over all
0969  * levels, since surface deduplication can potentionally elide bounding
0970  * surfaces at more deeply embedded levels.
0971  */
0972 CELER_FUNCTION real_type OrangeTrackView::find_safety()
0973 {
0974     CELER_EXPECT(!this->is_on_boundary());
0975 
0976     TrackerVisitor visit_tracker{params_};
0977 
0978     real_type min_safety_dist = numeric_limits<real_type>::infinity();
0979 
0980     for (auto lev : range(LevelId{this->level() + 1}))
0981     {
0982         auto lsa = this->make_lsa(lev);
0983         auto sd = visit_tracker(
0984             [&lsa](auto&& t) { return t.safety(lsa.pos(), lsa.vol()); },
0985             lsa.universe());
0986         min_safety_dist = celeritas::min(min_safety_dist, sd);
0987     }
0988     return min_safety_dist;
0989 }
0990 
0991 //---------------------------------------------------------------------------//
0992 /*!
0993  * Find the distance to the nearest nearby boundary.
0994  *
0995  * Since we currently support only "simple" safety distances, we can't
0996  * eliminate anything by checking only nearby surfaces.
0997  */
0998 CELER_FUNCTION real_type OrangeTrackView::find_safety(real_type)
0999 {
1000     return this->find_safety();
1001 }
1002 
1003 //---------------------------------------------------------------------------//
1004 /*!
1005  * Get a reference to the current volume, or to world volume if outside.
1006  */
1007 CELER_FUNCTION Span<Sense> OrangeTrackView::make_temp_sense() const
1008 {
1009     auto const max_faces = params_.scalars.max_faces;
1010     auto offset = track_slot_.get() * max_faces;
1011     return states_.temp_sense[AllItems<Sense, MemSpace::native>{}].subspan(
1012         offset, max_faces);
1013 }
1014 
1015 //---------------------------------------------------------------------------//
1016 /*!
1017  * Set up intersection scratch space.
1018  */
1019 CELER_FUNCTION detail::TempNextFace OrangeTrackView::make_temp_next() const
1020 {
1021     auto const max_isect = params_.scalars.max_intersections;
1022     auto offset = track_slot_.get() * max_isect;
1023 
1024     detail::TempNextFace result;
1025     result.face = states_.temp_face[AllItems<FaceId>{}].data() + offset;
1026     result.distance = states_.temp_distance[AllItems<real_type>{}].data()
1027                       + offset;
1028     result.isect = states_.temp_isect[AllItems<size_type>{}].data() + offset;
1029     result.size = max_isect;
1030 
1031     return result;
1032 }
1033 
1034 //---------------------------------------------------------------------------//
1035 /*!
1036  * Create a local state.
1037  */
1038 CELER_FUNCTION detail::LocalState
1039 OrangeTrackView::make_local_state(LevelId level) const
1040 {
1041     detail::LocalState local;
1042 
1043     auto lsa = this->make_lsa(level);
1044 
1045     local.pos = lsa.pos();
1046     local.dir = lsa.dir();
1047     local.volume = lsa.vol();
1048     if (level == this->surface_level())
1049     {
1050         local.surface = {this->surf(), this->sense()};
1051     }
1052     else
1053     {
1054         local.surface = {};
1055     }
1056     local.temp_sense = this->make_temp_sense();
1057     local.temp_next = this->make_temp_next();
1058     return local;
1059 }
1060 
1061 //---------------------------------------------------------------------------//
1062 /*!
1063  * Whether any next step has been calculated.
1064  */
1065 CELER_FORCEINLINE_FUNCTION bool OrangeTrackView::has_next_step() const
1066 {
1067     return this->next_step() != 0;
1068 }
1069 
1070 //---------------------------------------------------------------------------//
1071 /*!
1072  * Whether the next intersecting surface has been found.
1073  */
1074 CELER_FORCEINLINE_FUNCTION bool OrangeTrackView::has_next_surface() const
1075 {
1076     return static_cast<bool>(states_.next_surf[track_slot_]);
1077 }
1078 
1079 //---------------------------------------------------------------------------//
1080 /*!
1081  * Reset the next distance-to-boundary and surface.
1082  */
1083 CELER_FUNCTION void OrangeTrackView::clear_next()
1084 {
1085     this->next_step(0);
1086     states_.next_surf[track_slot_] = {};
1087     CELER_ENSURE(!this->has_next_step() && !this->has_next_surface());
1088 }
1089 
1090 //---------------------------------------------------------------------------//
1091 /*!
1092  * Assign the surface on the current level.
1093  */
1094 CELER_FUNCTION void
1095 OrangeTrackView::surface(LevelId level, detail::OnLocalSurface surf)
1096 {
1097     states_.surface_level[track_slot_] = level;
1098     states_.surf[track_slot_] = surf.id();
1099     states_.sense[track_slot_] = surf.unchecked_sense();
1100 }
1101 
1102 //---------------------------------------------------------------------------//
1103 /*!
1104  * Clear the surface on the current level.
1105  */
1106 CELER_FUNCTION void OrangeTrackView::clear_surface()
1107 {
1108     states_.surface_level[track_slot_] = {};
1109     CELER_ENSURE(!this->is_on_boundary());
1110 }
1111 
1112 //---------------------------------------------------------------------------//
1113 /*!
1114  * Make a LevelStateAccessor for the current thread and level.
1115  */
1116 CELER_FORCEINLINE_FUNCTION auto OrangeTrackView::make_lsa() const -> LSA
1117 {
1118     return this->make_lsa(this->level());
1119 }
1120 
1121 //---------------------------------------------------------------------------//
1122 /*!
1123  * Make a LevelStateAccessor for the current thread and a given level.
1124  */
1125 CELER_FORCEINLINE_FUNCTION auto
1126 OrangeTrackView::make_lsa(LevelId level) const -> LSA
1127 {
1128     return LSA(&states_, track_slot_, level);
1129 }
1130 
1131 //---------------------------------------------------------------------------//
1132 /*!
1133  * Get the daughter ID for the given volume in the given universe.
1134  *
1135  * \return DaughterId or {} if the current volume is a leaf.
1136  */
1137 CELER_FUNCTION DaughterId OrangeTrackView::get_daughter(LSA const& lsa)
1138 {
1139     TrackerVisitor visit_tracker{params_};
1140     return visit_tracker([&lsa](auto&& t) { return t.daughter(lsa.vol()); },
1141                          lsa.universe());
1142 }
1143 
1144 //---------------------------------------------------------------------------//
1145 /*!
1146  * Get the transform ID for the given daughter.
1147  */
1148 CELER_FUNCTION TransformId OrangeTrackView::get_transform(DaughterId daughter_id)
1149 {
1150     CELER_EXPECT(daughter_id);
1151     return params_.daughters[daughter_id].transform_id;
1152 }
1153 
1154 //---------------------------------------------------------------------------//
1155 /*!
1156  * Get the transform ID for the given daughter.
1157  */
1158 CELER_FUNCTION TransformId OrangeTrackView::get_transform(LevelId lev)
1159 {
1160     CELER_EXPECT(lev < this->level());
1161     LSA lsa(&states_, track_slot_, lev);
1162     return this->get_transform(this->get_daughter(lsa));
1163 }
1164 
1165 //---------------------------------------------------------------------------//
1166 }  // namespace celeritas