orange/univ/RectArrayTracker.hh

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/univ/RectArrayTracker.hh
0007 //---------------------------------------------------------------------------//
0008 #pragma once
0009
0010 #include "corecel/Assert.hh"
0011 #include "corecel/data/HyperslabIndexer.hh"
0012 #include "corecel/grid/NonuniformGrid.hh"
0013 #include "corecel/math/Algorithms.hh"
0014 #include "orange/OrangeData.hh"
0015
0016 #include "detail/RaggedRightIndexer.hh"
0017 #include "detail/Types.hh"
0018 #include "detail/Utils.hh"
0019
0020 namespace celeritas
0021 {
0022 //---------------------------------------------------------------------------//
0023 /*!
0024  * Track a particle within an axes-aligned rectilinear grid.
0025  */
0026 class RectArrayTracker
0027 {
0028   public:
0029     //!@{
0030     //! \name Type aliases
0031     using ParamsRef = NativeCRef<OrangeParamsData>;
0032     using Initialization = detail::Initialization;
0033     using Intersection = detail::Intersection;
0034     using LocalState = detail::LocalState;
0035     using Grid = NonuniformGrid<real_type>;
0036     using VolumeIndexer = HyperslabIndexer<3>;
0037     using VolumeInverseIndexer = HyperslabInverseIndexer<3>;
0038     using SurfaceIndexer = detail::RaggedRightIndexer<3>;
0039     using SurfaceInverseIndexer = detail::RaggedRightInverseIndexer<3>;
0040     using Coords = Array<size_type, 3>;
0041     //!@}
0042
0043   public:
0044     // Construct with parameters (unit definitions and this one's ID)
0045     inline CELER_FUNCTION
0046     RectArrayTracker(ParamsRef const& params, RectArrayId rid);
0047
0048     //// ACCESSORS ////
0049
0050     //! Number of local volumes
0051     CELER_FUNCTION LocalVolumeId::size_type num_volumes() const
0052     {
0053         return record_.daughters.size();
0054     }
0055
0056     //! Number of local surfaces
0057     CELER_FUNCTION LocalSurfaceId::size_type num_surfaces() const
0058     {
0059         size_type num_surfs = 0;
0060         for (auto ax : range(Axis::size_))
0061         {
0062             num_surfs += record_.dims[to_int(ax)] + 1;
0063         }
0064         return num_surfs;
0065     }
0066
0067     // DaughterId of universe embedded in a given volume
0068     inline CELER_FUNCTION DaughterId daughter(LocalVolumeId vol) const;
0069
0070     ////// OPERATIONS ////
0071
0072     // Find the local volume from a position
0073     inline CELER_FUNCTION Initialization
0074     initialize(LocalState const& state) const;
0075
0076     // Calculate distance-to-intercept for the next surface
0077     inline CELER_FUNCTION Intersection intersect(LocalState const& state) const;
0078
0079     // Calculate distance-to-intercept for the next surface, with max distance
0080     inline CELER_FUNCTION Intersection intersect(LocalState const& state,
0081                                                  real_type max_dist) const;
0082
0083     // Find the local volume given a post-crossing state
0084     inline CELER_FUNCTION Initialization
0085     cross_boundary(LocalState const& state) const;
0086
0087     // Calculate closest distance to a surface in any direction
0088     inline CELER_FUNCTION real_type safety(Real3 const& pos,
0089                                            LocalVolumeId vol) const;
0090
0091     // Calculate the local surface normal
0092     inline CELER_FUNCTION Real3 normal(Real3 const& pos,
0093                                        LocalSurfaceId surf) const;
0094
0095   private:
0096     //// DATA ////
0097     ParamsRef const& params_;
0098     RectArrayRecord const& record_;
0099
0100     //// METHODS ////
0101
0102     // Calculate distance-to-intercept for the next surface.
0103     template<class F>
0104     inline CELER_FUNCTION Intersection intersect_impl(LocalState const&,
0105                                                       F) const;
0106
0107     // Find the index of axis (x/y/z) we are about to cross
0108     inline CELER_FUNCTION size_type find_surface_axis_idx(LocalSurfaceId s) const;
0109
0110     // Create Grid object for a given axis
0111     inline CELER_FUNCTION Grid make_grid(Axis ax) const;
0112 };
0113
0114 //---------------------------------------------------------------------------//
0115 // INLINE DEFINITIONS
0116 //---------------------------------------------------------------------------//
0117 /*!
0118  * Construct with reference to persistent parameter data.
0119  */
0120 CELER_FUNCTION
0121 RectArrayTracker::RectArrayTracker(ParamsRef const& params, RectArrayId rid)
0122     : params_(params), record_(params.rect_arrays[rid])
0123 {
0124     CELER_EXPECT(params_);
0125 }
0126
0127 //---------------------------------------------------------------------------//
0128 /*!
0129  * Find the local volume from a position.
0130  *
0131  * To avoid edge cases and inconsistent logical/physical states, it is
0132  * prohibited to initialize from an arbitrary point directly onto a surface.
0133  */
0134 CELER_FUNCTION auto
0135 RectArrayTracker::initialize(LocalState const& state) const -> Initialization
0136 {
0137     CELER_EXPECT(params_);
0138     CELER_EXPECT(!state.surface && !state.volume);
0139
0140     Coords coords;
0141
0142     for (auto ax : range(Axis::size_))
0143     {
0144         auto const& pos = state.pos[to_int(ax)];
0145         auto grid = this->make_grid(ax);
0146
0147         if (pos < grid.front() || pos > grid.back())
0148         {
0149             // Outside the rect array
0150             return {};
0151         }
0152         else
0153         {
0154             size_type index = grid.find(pos);
0155             if (grid[index] == pos)
0156             {
0157                 // Initialization exactly on an edge is prohibited
0158                 return {};
0159             }
0160             coords[to_int(ax)] = index;
0161         }
0162     }
0163
0164     VolumeIndexer to_index(record_.dims);
0165     return {LocalVolumeId{to_index(coords)}, {}};
0166 }
0167
0168 //---------------------------------------------------------------------------//
0169 /*!
0170  * Find the local volume given a post-crossing state.
0171  */
0172 CELER_FUNCTION auto
0173 RectArrayTracker::cross_boundary(LocalState const& state) const -> Initialization
0174 {
0175     CELER_EXPECT(state.surface && state.volume);
0176
0177     // Find the coords of the current volume
0178     VolumeIndexer to_index(record_.dims);
0179     VolumeInverseIndexer to_coords(record_.dims);
0180     auto coords = to_coords(state.volume.unchecked_get());
0181     auto ax_idx = this->find_surface_axis_idx(state.surface.id());
0182
0183     // Due to surface deduplication, crossing out of rect arrays is not
0184     // possible
0185     CELER_ASSERT(
0186         !(coords[ax_idx] == 0 && state.surface.sense() == Sense::inside)
0187         && !(coords[ax_idx] == record_.dims[ax_idx] - 1
0188              && state.surface.sense() == Sense::outside));
0189
0190     // Value for incrementing the axial coordinate upon crossing
0191     // NOTE: that surface sense is currently the POST crossing value
0192     int inc = (state.surface.sense() == Sense::outside) ? 1 : -1;
0193
0194     coords[ax_idx] += inc;
0195     return {LocalVolumeId(to_index(coords)), state.surface};
0196 }
0197
0198 //---------------------------------------------------------------------------//
0199 /*!
0200  * Calculate distance-to-intercept for the next surface.
0201  */
0202 CELER_FUNCTION auto
0203 RectArrayTracker::intersect(LocalState const& state) const -> Intersection
0204 {
0205     Intersection result = this->intersect_impl(state, detail::IsFinite{});
0206     return result;
0207 }
0208
0209 //---------------------------------------------------------------------------//
0210 /*!
0211  * Calculate distance-to-intercept for the next surface, with max distance.
0212  */
0213 CELER_FUNCTION auto
0214 RectArrayTracker::intersect(LocalState const& state,
0215                             real_type max_dist) const -> Intersection
0216 {
0217     CELER_EXPECT(max_dist > 0);
0218     Intersection result
0219         = this->intersect_impl(state, detail::IsNotFurtherThan{max_dist});
0220     if (!result)
0221     {
0222         result.distance = max_dist;
0223     }
0224     return result;
0225 }
0226
0227 //---------------------------------------------------------------------------//
0228 /*!
0229  * Calculate nearest distance to a surface in any direction.
0230  *
0231  * On an axis-aligned rectlinear grid the minimum distance to any surface is
0232  * always always occurs along a line parallel to an axis.
0233  */
0234 CELER_FUNCTION real_type RectArrayTracker::safety(Real3 const& pos,
0235                                                   LocalVolumeId volid) const
0236 {
0237     CELER_EXPECT(volid && volid.get() < this->num_volumes());
0238
0239     VolumeInverseIndexer to_coords(record_.dims);
0240     auto coords = to_coords(volid.unchecked_get());
0241
0242     real_type min_dist = numeric_limits<real_type>::infinity();
0243
0244     for (auto ax : range(Axis::size_))
0245     {
0246         auto grid = this->make_grid(ax);
0247         for (auto i : range(2))
0248         {
0249             auto target_coord = coords[to_int(ax)] + i;
0250             real_type target = grid[target_coord];
0251             min_dist = min(min_dist, std::fabs(pos[to_int(ax)] - target));
0252         }
0253     }
0254
0255     CELER_ENSURE(min_dist >= 0
0256                  && min_dist < numeric_limits<real_type>::infinity());
0257     return min_dist;
0258 }
0259
0260 //---------------------------------------------------------------------------//
0261 /*!
0262  * Calculate the local surface normal.
0263  */
0264 CELER_FUNCTION auto
0265 RectArrayTracker::normal(Real3 const&, LocalSurfaceId surf) const -> Real3
0266 {
0267     CELER_EXPECT(surf && surf.get() < this->num_surfaces());
0268     size_type ax = this->find_surface_axis_idx(surf);
0269
0270     Real3 normal{0, 0, 0};
0271     normal[ax] = 1;
0272
0273     return normal;
0274 }
0275
0276 //---------------------------------------------------------------------------//
0277 /*!
0278  * DaughterId of universe embedded in a given volume.
0279  */
0280 CELER_FORCEINLINE_FUNCTION DaughterId
0281 RectArrayTracker::daughter(LocalVolumeId vol) const
0282 {
0283     CELER_EXPECT(vol && vol.get() < this->num_volumes());
0284     return record_.daughters[vol];
0285 }
0286
0287 //---------------------------------------------------------------------------//
0288 // PRIVATE INLINE DEFINITIONS
0289 //---------------------------------------------------------------------------//
0290 /*!
0291  * Calculate distance-to-intercept for the next surface.
0292  */
0293 template<class F>
0294 CELER_FUNCTION auto
0295 RectArrayTracker::intersect_impl(LocalState const& state,
0296                                  F is_valid) const -> Intersection
0297 {
0298     CELER_EXPECT(state.volume && !state.temp_sense.empty());
0299
0300     auto coords
0301         = VolumeInverseIndexer{record_.dims}(state.volume.unchecked_get());
0302
0303     Intersection result;
0304     Sense sense;
0305     SurfaceIndexer to_index(record_.surface_indexer_data);
0306
0307     for (auto ax : range(Axis::size_))
0308     {
0309         auto dir = state.dir[to_int(ax)];
0310
0311         // Ignore any stationary axis
0312         if (dir == 0)
0313         {
0314             continue;
0315         }
0316
0317         auto target_coord = coords[to_int(ax)] + static_cast<int>(dir > 0);
0318
0319         auto target_value = this->make_grid(ax)[target_coord];
0320
0321         real_type dist
0322             = (target_value - static_cast<real_type>(state.pos[to_int(ax)]))
0323               / state.dir[to_int(ax)];
0324
0325         if (dist > 0 && is_valid(dist) && dist < result.distance)
0326         {
0327             result.distance = dist;
0328
0329             sense = dir > 0 ? Sense::inside : Sense::outside;
0330             auto local_surface = LocalSurfaceId(
0331                 to_index({static_cast<size_type>(to_int(ax)), target_coord}));
0332             result.surface = detail::OnLocalSurface(local_surface, sense);
0333         }
0334     }
0335
0336     return result;
0337 }
0338
0339 //---------------------------------------------------------------------------//
0340 /*!
0341  * Find the index of axis (x/y/z) we are about to cross:
0342  */
0343 CELER_FUNCTION size_type
0344 RectArrayTracker::find_surface_axis_idx(LocalSurfaceId s) const
0345 {
0346     SurfaceInverseIndexer to_axis(record_.surface_indexer_data);
0347     return to_axis(s.unchecked_get())[0];
0348 }
0349
0350 //---------------------------------------------------------------------------//
0351 /*!
0352  * Create Grid object for a given axis.
0353  */
0354 CELER_FUNCTION RectArrayTracker::Grid RectArrayTracker::make_grid(Axis ax) const
0355 {
0356     return Grid(record_.grid[to_int(ax)], params_.reals);
0357 }
0358
0359 //---------------------------------------------------------------------------//
0360 }  // namespace celeritas