orange/univ/RectArrayTracker.hh

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