alongstep/detail/PropagationApplier.hh

0001 //----------------------------------*-C++-*----------------------------------//
0002 // Copyright 2023-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 celeritas/global/alongstep/detail/PropagationApplier.hh
0007 //---------------------------------------------------------------------------//
0008 #pragma once
0009
0010 #include <type_traits>
0011
0012 #include "corecel/math/Algorithms.hh"
0013 #include "corecel/sys/KernelTraits.hh"
0014 #include "celeritas/global/CoreTrackView.hh"
0015
0016 #if !CELER_DEVICE_COMPILE
0017 #    include "corecel/io/Logger.hh"
0018 #    include "corecel/io/Repr.hh"
0019 #endif
0020
0021 namespace celeritas
0022 {
0023 namespace detail
0024 {
0025 //---------------------------------------------------------------------------//
0026 /*!
0027  * Apply propagation over the step (implementation).
0028  */
0029 template<class MP>
0030 struct PropagationApplierBaseImpl
0031 {
0032     inline CELER_FUNCTION void operator()(CoreTrackView& track);
0033
0034     MP make_propagator;
0035 };
0036
0037 //---------------------------------------------------------------------------//
0038 /*!
0039  * Apply propagation over the step.
0040  *
0041  * \tparam MP Propagator factory
0042  *
0043  * MP should be a function-like object:
0044  * \code Propagator(*)(CoreTrackView const&) \endcode
0045  *
0046  * This class is partially specialized with a second template argument to
0047  * extract any launch bounds from the MP class. TODO: we could probably inherit
0048  * from a helper class to pull in those constants (if available).
0049  */
0050 template<class MP, typename = void>
0051 struct PropagationApplier : public PropagationApplierBaseImpl<MP>
0052 {
0053     CELER_FUNCTION PropagationApplier(MP&& mp)
0054         : PropagationApplierBaseImpl<MP>{celeritas::forward<MP>(mp)}
0055     {
0056     }
0057 };
0058
0059 template<class MP>
0060 struct PropagationApplier<MP, std::enable_if_t<kernel_max_blocks_min_warps<MP>>>
0061     : public PropagationApplierBaseImpl<MP>
0062 {
0063     static constexpr int max_block_size = MP::max_block_size;
0064     static constexpr int min_warps_per_eu = MP::min_warps_per_eu;
0065
0066     CELER_FUNCTION PropagationApplier(MP&& mp)
0067         : PropagationApplierBaseImpl<MP>{celeritas::forward<MP>(mp)}
0068     {
0069     }
0070 };
0071
0072 template<class MP>
0073 struct PropagationApplier<MP, std::enable_if_t<kernel_max_blocks<MP>>>
0074     : public PropagationApplierBaseImpl<MP>
0075 {
0076     static constexpr int max_block_size = MP::max_block_size;
0077
0078     CELER_FUNCTION PropagationApplier(MP&& mp)
0079         : PropagationApplierBaseImpl<MP>{celeritas::forward<MP>(mp)}
0080     {
0081     }
0082 };
0083
0084 //---------------------------------------------------------------------------//
0085 // DEDUCTION GUIDES
0086 //---------------------------------------------------------------------------//
0087 template<class MP>
0088 CELER_FUNCTION PropagationApplier(MP&&) -> PropagationApplier<MP>;
0089
0090 //---------------------------------------------------------------------------//
0091 // INLINE DEFINITIONS
0092 //---------------------------------------------------------------------------//
0093 template<class MP>
0094 CELER_FUNCTION void
0095 PropagationApplierBaseImpl<MP>::operator()(CoreTrackView& track)
0096 {
0097     auto sim = track.make_sim_view();
0098     if (sim.step_length() == 0)
0099     {
0100         // Track is stopped: no movement or energy loss will happen
0101         // (could be a stopped positron waiting for annihilation, or a
0102         // particle waiting to decay?)
0103         CELER_ASSERT(track.make_particle_view().is_stopped());
0104         CELER_ASSERT(sim.post_step_action()
0105                      == track.make_physics_view().scalars().discrete_action());
0106         CELER_ASSERT(track.make_physics_view().has_at_rest());
0107         return;
0108     }
0109
0110     bool tracks_can_loop;
0111     Propagation p;
0112     {
0113 #if CELERITAS_DEBUG
0114         Real3 const orig_pos = track.make_geo_view().pos();
0115 #endif
0116         auto propagate = make_propagator(track);
0117         p = propagate(sim.step_length());
0118         tracks_can_loop = propagate.tracks_can_loop();
0119         CELER_ASSERT(p.distance > 0);
0120 #if CELERITAS_DEBUG
0121         if (CELER_UNLIKELY(track.make_geo_view().pos() == orig_pos))
0122         {
0123             // This unusual case happens when the step length is less than
0124             // machine epsilon compared to the actual position. This case seems
0125             // to occur in vecgeom when "stuck" on a boundary, and when in a
0126             // field taking a small step when the track's position has a large
0127             // magnitude.
0128 #    if !CELER_DEVICE_COMPILE
0129             CELER_LOG_LOCAL(error)
0130                 << "Propagation of step length " << repr(sim.step_length())
0131                 << " due to post-step action "
0132                 << sim.post_step_action().unchecked_get()
0133                 << " leading to distance " << repr(p.distance)
0134                 << (p.boundary  ? " (boundary hit)"
0135                     : p.looping ? " (**LOOPING**)"
0136                                 : "")
0137                 << " failed to change position";
0138 #    endif
0139             track.apply_errored();
0140             return;
0141         }
0142 #endif
0143     }
0144
0145     if (tracks_can_loop)
0146     {
0147         sim.update_looping(p.looping);
0148     }
0149     if (tracks_can_loop && p.looping)
0150     {
0151         // The track is looping, i.e. progressing little over many
0152         // integration steps in the field propagator (likely a low energy
0153         // particle in a low density material/strong magnetic field).
0154         sim.step_length(p.distance);
0155
0156         // Kill the track if it's stable and below the threshold energy or
0157         // above the threshold number of steps allowed while looping.
0158         sim.post_step_action([&track, &sim] {
0159             auto particle = track.make_particle_view();
0160             if (particle.is_stable()
0161                 && sim.is_looping(particle.particle_id(), particle.energy()))
0162             {
0163 #if !CELER_DEVICE_COMPILE
0164                 CELER_LOG_LOCAL(debug)
0165                     << "Track (pid=" << particle.particle_id().get()
0166                     << ", E=" << particle.energy().value() << ' '
0167                     << ParticleTrackView::Energy::unit_type::label()
0168                     << ") is looping after " << sim.num_looping_steps()
0169                     << " steps";
0170 #endif
0171                 return track.tracking_cut_action();
0172             }
0173             return track.propagation_limit_action();
0174         }());
0175     }
0176     else if (p.boundary)
0177     {
0178         // Stopped at a geometry boundary: this is the new step action.
0179         CELER_ASSERT(p.distance <= sim.step_length());
0180         sim.step_length(p.distance);
0181         sim.post_step_action(track.boundary_action());
0182     }
0183     else if (p.distance < sim.step_length())
0184     {
0185         // Some tracks may get stuck on a boundary and fail to move at
0186         // all in the field propagator, and will get bumped a small
0187         // distance. This primarily occurs with reentrant tracks on a
0188         // boundary with VecGeom.
0189         sim.step_length(p.distance);
0190         sim.post_step_action(track.propagation_limit_action());
0191     }
0192 }
0193
0194 //---------------------------------------------------------------------------//
0195 }  // namespace detail
0196 }  // namespace celeritas