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