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