celeritas/phys/InteractionApplier.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/phys/InteractionApplier.hh
0006 //---------------------------------------------------------------------------//
0007 #pragma once
0008
0009 #include "corecel/Macros.hh"
0010 #include "corecel/cont/Span.hh"
0011 #include "corecel/sys/KernelTraits.hh"
0012 #include "celeritas/Quantities.hh"
0013 #include "celeritas/Types.hh"
0014 #include "celeritas/geo/GeoFwd.hh"
0015 #include "celeritas/global/CoreTrackView.hh"
0016 #include "celeritas/track/SimTrackView.hh"
0017
0018 #include "CutoffView.hh"
0019 #include "Interaction.hh"
0020 #include "ParticleTrackView.hh"
0021 #include "ParticleView.hh"
0022 #include "PhysicsData.hh"
0023 #include "PhysicsStepView.hh"
0024 #include "PhysicsTrackView.hh"
0025 #include "Secondary.hh"
0026
0027 namespace celeritas
0028 {
0029 //---------------------------------------------------------------------------//
0030 /*!
0031  * Wrap an Interaction executor to apply it to a track.
0032  *
0033  * The function F must take a \c CoreTrackView and return a \c Interaction
0034  */
0035 template<class F>
0036 struct InteractionApplierBaseImpl
0037 {
0038     //// DATA ////
0039
0040     F sample_interaction;
0041
0042     //// METHODS ////
0043
0044     CELER_FUNCTION void operator()(celeritas::CoreTrackView const&);
0045 };
0046
0047 //---------------------------------------------------------------------------//
0048 /*!
0049  *
0050  * This class is partially specialized with a second template argument to
0051  * extract any launch bounds from the functor class. TODO: we could probably
0052  * inherit from a helper class to pull in those constants (if available).
0053  */
0054 template<class F, typename = void>
0055 struct InteractionApplier : public InteractionApplierBaseImpl<F>
0056 {
0057     CELER_FUNCTION InteractionApplier(F&& f)
0058         : InteractionApplierBaseImpl<F>{celeritas::forward<F>(f)}
0059     {
0060     }
0061 };
0062
0063 template<class F>
0064 struct InteractionApplier<F, std::enable_if_t<kernel_max_blocks_min_warps<F>>>
0065     : public InteractionApplierBaseImpl<F>
0066 {
0067     static constexpr int max_block_size = F::max_block_size;
0068     static constexpr int min_warps_per_eu = F::min_warps_per_eu;
0069
0070     CELER_FUNCTION InteractionApplier(F&& f)
0071         : InteractionApplierBaseImpl<F>{celeritas::forward<F>(f)}
0072     {
0073     }
0074 };
0075
0076 template<class F>
0077 struct InteractionApplier<F, std::enable_if_t<kernel_max_blocks<F>>>
0078     : public InteractionApplierBaseImpl<F>
0079 {
0080     static constexpr int max_block_size = F::max_block_size;
0081
0082     CELER_FUNCTION InteractionApplier(F&& f)
0083         : InteractionApplierBaseImpl<F>{celeritas::forward<F>(f)}
0084     {
0085     }
0086 };
0087
0088 //---------------------------------------------------------------------------//
0089 // DEDUCTION GUIDES
0090 //---------------------------------------------------------------------------//
0091 template<class F>
0092 CELER_FUNCTION InteractionApplier(F&&) -> InteractionApplier<F>;
0093
0094 //---------------------------------------------------------------------------//
0095 // INLINE DEFINITIONS
0096 //---------------------------------------------------------------------------//
0097 /*!
0098  * Sample an interaction and apply to the track view.
0099  *
0100  * The given track *must* be an active track with the correct step limit action
0101  * ID.
0102  */
0103 template<class F>
0104 CELER_FUNCTION void
0105 InteractionApplierBaseImpl<F>::operator()(celeritas::CoreTrackView const& track)
0106 {
0107     Interaction result = this->sample_interaction(track);
0108
0109     auto sim = track.sim();
0110     if (CELER_UNLIKELY(result.action == Interaction::Action::failed))
0111     {
0112         auto phys = track.physics();
0113         // Particle already moved to the collision site, but an out-of-memory
0114         // (allocation failure) occurred. Someday we can add error handling,
0115         // but for now use the "failure" action in the physics and set the step
0116         // limit to zero since it needs to interact again at this location.
0117         sim.step_limit({0, phys.scalars().failure_action()});
0118         return;
0119     }
0120     else if (!result.changed())
0121     {
0122         return;
0123     }
0124
0125     // Scattered or absorbed
0126     {
0127         // Update post-step energy
0128         auto particle = track.particle();
0129         particle.energy(result.energy);
0130     }
0131
0132     if (result.action != Interaction::Action::absorbed)
0133     {
0134         // Update direction
0135         auto geo = track.geometry();
0136         geo.set_dir(result.direction);
0137     }
0138     else
0139     {
0140         // Mark particle as dead
0141         sim.status(TrackStatus::killed);
0142     }
0143
0144     real_type deposition = result.energy_deposition.value();
0145     auto cutoff = track.cutoff();
0146     if (cutoff.apply_post_interaction())
0147     {
0148         // Kill secondaries with energies below the production cut
0149         for (auto& secondary : result.secondaries)
0150         {
0151             if (cutoff.apply(secondary))
0152             {
0153                 // Secondary is an electron, positron or gamma with energy
0154                 // below the production cut -- deposit the energy locally
0155                 // and clear the secondary
0156                 deposition += secondary.energy.value();
0157                 auto sec_par = track.particle_record(secondary.particle_id);
0158                 if (sec_par.is_antiparticle())
0159                 {
0160                     // Conservation of energy for positrons
0161                     deposition += 2 * sec_par.mass().value();
0162                 }
0163                 secondary = {};
0164             }
0165         }
0166     }
0167     auto phys = track.physics_step();
0168     phys.deposit_energy(units::MevEnergy{deposition});
0169     phys.secondaries(result.secondaries);
0170 }
0171
0172 //---------------------------------------------------------------------------//
0173 }  // namespace celeritas