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 RectArrayTracker::initialize(LocalState const& state) const
0134     -> 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
0173     -> 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 RectArrayTracker::intersect(LocalState const& state) const
0203     -> 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, real_type max_dist) const
0215     -> 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 vol_id) const
0236 {
0237     CELER_EXPECT(vol_id && vol_id.get() < this->num_volumes());
0238
0239     VolumeInverseIndexer to_coords(record_.dims);
0240     auto coords = to_coords(vol_id.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, F is_valid) const
0296     -> 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