track/detail/ProcessSecondariesExecutor.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/track/detail/ProcessSecondariesExecutor.hh
0006 //---------------------------------------------------------------------------//
0007 #pragma once
0008
0009 #include "corecel/Assert.hh"
0010 #include "corecel/Macros.hh"
0011 #include "corecel/cont/Span.hh"
0012 #include "corecel/sys/ThreadId.hh"
0013 #include "celeritas/Quantities.hh"
0014 #include "celeritas/Types.hh"
0015 #include "celeritas/global/CoreTrackData.hh"
0016 #include "celeritas/global/CoreTrackView.hh"
0017 #include "celeritas/phys/ParticleData.hh"
0018 #include "celeritas/phys/Secondary.hh"
0019
0020 #include "../CoreStateCounters.hh"
0021
0022 namespace celeritas
0023 {
0024 namespace detail
0025 {
0026 //---------------------------------------------------------------------------//
0027 /*!
0028  * Create track initializers from secondaries.
0029  */
0030 struct ProcessSecondariesExecutor
0031 {
0032     //// TYPES ////
0033
0034     using ParamsPtr = CRefPtr<CoreParamsData, MemSpace::native>;
0035     using StatePtr = RefPtr<CoreStateData, MemSpace::native>;
0036
0037     //// DATA ////
0038
0039     ParamsPtr params;
0040     StatePtr state;
0041     CoreStateCounters counters;
0042
0043     //// FUNCTIONS ////
0044
0045     // Determine which tracks are alive and count secondaries
0046     inline CELER_FUNCTION void operator()(TrackSlotId tid) const;
0047
0048     CELER_FORCEINLINE_FUNCTION void operator()(ThreadId tid) const
0049     {
0050         // The grid size should be equal to the state size and no thread/slot
0051         // remapping should be performed
0052         return (*this)(TrackSlotId{tid.unchecked_get()});
0053     }
0054 };
0055
0056 //---------------------------------------------------------------------------//
0057 /*!
0058  * Create track initializers from secondaries.
0059  *
0060  * This kernel is executed with a grid size equal to the number of track
0061  * slots, so ThreadId should be equal to TrackSlotId. No remapping should be
0062  * done.
0063  */
0064 CELER_FUNCTION void
0065 ProcessSecondariesExecutor::operator()(TrackSlotId tid) const
0066 {
0067     CELER_EXPECT(tid < state->size());
0068
0069     CoreTrackView track(*params, *state, tid);
0070
0071     auto sim = track.sim();
0072     if (sim.status() == TrackStatus::inactive)
0073     {
0074         // Do not create secondaries from stale data on inactive tracks
0075         return;
0076     }
0077
0078     // Offset in the vector of track initializers
0079     auto& data = state->init;
0080     CELER_ASSERT(data.secondary_counts[tid] <= counters.num_secondaries);
0081     size_type offset = counters.num_secondaries - data.secondary_counts[tid];
0082
0083     // Save the parent ID since it will be overwritten if a secondary is
0084     // initialized in this slot
0085     TrackId const track_id{sim.track_id()};
0086
0087     for (auto const& secondary : track.physics_step().secondaries())
0088     {
0089         if (secondary)
0090         {
0091             CELER_ASSERT(secondary.energy > zero_quantity()
0092                          && is_soft_unit_vector(secondary.direction));
0093
0094             // Particles should not be making secondaries while crossing a
0095             // surface
0096             auto geo = track.geometry();
0097             CELER_ASSERT(!geo.is_on_boundary());
0098
0099             // Create a track initializer from the secondary
0100             TrackInitializer ti;
0101             ti.sim.track_id = make_track_id(params->init, data, sim.event_id());
0102             ti.sim.primary_id = sim.primary_id();
0103             ti.sim.parent_id = track_id;
0104             ti.sim.event_id = sim.event_id();
0105             ti.sim.time = sim.time();
0106             ti.sim.weight = sim.weight();
0107             ti.geo.pos = geo.pos();
0108             ti.geo.dir = secondary.direction;
0109             ti.particle.particle_id = secondary.particle_id;
0110             ti.particle.energy = secondary.energy;
0111             CELER_ASSERT(ti);
0112
0113             if (sim.track_id() == track_id && sim.status() != TrackStatus::alive
0114                 && params->init.track_order != TrackOrder::init_charge)
0115             {
0116                 /*!
0117                  * Skip in-place initialization when tracks are partitioned by
0118                  * charge to reduce the amount of mixing
0119                  *
0120                  * \todo Consider allowing this if the parent's charge is the
0121                  * same as the secondary's
0122                  */
0123
0124                 // The parent was killed, so initialize the first secondary in
0125                 // the parent's track slot. Keep the parent's geometry state
0126                 // but get the direction from the secondary.
0127                 ti.geo.parent = tid;
0128                 track = ti;
0129             }
0130             else
0131             {
0132                 CELER_ASSERT(offset > 0 && offset <= counters.num_initializers);
0133
0134                 if (offset <= min(counters.num_secondaries,
0135                                   counters.num_vacancies)
0136                     && (params->init.track_order != TrackOrder::init_charge
0137                         || sim.status() == TrackStatus::alive))
0138                 {
0139                     // Store the thread ID of the secondary's parent if the
0140                     // secondary will be initialized in the next step. If the
0141                     // initializers are partitioned by charge, the in-place
0142                     // initialization of the secondary is skipped, so another
0143                     // track might overwrite this state during initialization
0144                     // unless the track is alive.
0145                     ti.geo.parent = tid;
0146                 }
0147
0148                 // Store the track initializer
0149                 data.initializers[ItemId<TrackInitializer>{
0150                     counters.num_initializers - offset}]
0151                     = ti;
0152
0153                 --offset;
0154             }
0155         }
0156     }
0157
0158     if (sim.track_id() == track_id && sim.status() == TrackStatus::killed)
0159     {
0160         // Track is no longer used as part of transport
0161         sim.status(TrackStatus::inactive);
0162     }
0163     CELER_ENSURE(sim.status() != TrackStatus::killed);
0164 }
0165
0166 //---------------------------------------------------------------------------//
0167 }  // namespace detail
0168 }  // namespace celeritas