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 #include "celeritas/track/Utils.hh"
0020
0021 #include "../CoreStateCounters.hh"
0022
0023 namespace celeritas
0024 {
0025 namespace detail
0026 {
0027 //---------------------------------------------------------------------------//
0028 /*!
0029  * Create track initializers from secondaries.
0030  */
0031 struct ProcessSecondariesExecutor
0032 {
0033     //// TYPES ////
0034
0035     using ParamsPtr = CRefPtr<CoreParamsData, MemSpace::native>;
0036     using StatePtr = RefPtr<CoreStateData, MemSpace::native>;
0037
0038     //// DATA ////
0039
0040     ParamsPtr params;
0041     StatePtr state;
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 ProcessSecondariesExecutor::operator()(
0065     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     auto const& counters = track.counters();
0081     CELER_ASSERT(data.secondary_counts[tid] <= counters.num_secondaries);
0082     size_type offset = counters.num_secondaries - data.secondary_counts[tid];
0083
0084     // Save the parent ID since it will be overwritten if a secondary is
0085     // initialized in this slot
0086     TrackId const track_id{sim.track_id()};
0087
0088     for (auto const& secondary : track.physics_step().secondaries())
0089     {
0090         if (secondary)
0091         {
0092             CELER_ASSERT(secondary.energy > zero_quantity()
0093                          && is_soft_unit_vector(secondary.direction));
0094
0095             // Particles should not be making secondaries while crossing a
0096             // surface
0097             auto geo = track.geometry();
0098             CELER_ASSERT(!geo.is_on_boundary());
0099
0100             // Create a track initializer from the secondary
0101             TrackInitializer ti;
0102             ti.sim.track_id = make_track_id(params->init, data, sim.event_id());
0103             ti.sim.primary_id = sim.primary_id();
0104             ti.sim.parent_id = track_id;
0105             ti.sim.event_id = sim.event_id();
0106             ti.sim.time = sim.time();
0107             ti.sim.weight = sim.weight();
0108             ti.geo.pos = geo.pos();
0109             ti.geo.dir = secondary.direction;
0110             ti.particle.particle_id = secondary.particle_id;
0111             ti.particle.energy = secondary.energy;
0112             if constexpr (CELERITAS_RESEED == CELERITAS_RESEED_TRACK)
0113             {
0114                 // Branch the parent's RNG regardless of whether the
0115                 // secondary reuses the parent's slot, so that results are
0116                 // reproducible independent of the number of track slots
0117                 ti.rng = track.rng().branch();
0118             }
0119             CELER_ASSERT(ti);
0120
0121             if (sim.track_id() == track_id && sim.status() != TrackStatus::alive
0122                 && params->init.track_order != TrackOrder::init_charge)
0123             {
0124                 /*!
0125                  * Skip in-place initialization when tracks are partitioned by
0126                  * charge to reduce the amount of mixing
0127                  *
0128                  * \todo Consider allowing this if the parent's charge is the
0129                  * same as the secondary's
0130                  */
0131
0132                 // The parent was killed, so initialize the first secondary in
0133                 // the parent's track slot. Keep the parent's geometry state
0134                 // but get the direction from the secondary.
0135                 ti.geo.parent = tid;
0136                 track = ti;
0137             }
0138             else
0139             {
0140                 CELER_ASSERT(offset > 0 && offset <= counters.num_initializers);
0141
0142                 if (offset <= min(counters.num_secondaries,
0143                                   counters.num_vacancies)
0144                     && (params->init.track_order != TrackOrder::init_charge
0145                         || sim.status() == TrackStatus::alive))
0146                 {
0147                     // Store the thread ID of the secondary's parent if the
0148                     // secondary will be initialized in the next step. If the
0149                     // initializers are partitioned by charge, the in-place
0150                     // initialization of the secondary is skipped, so another
0151                     // track might overwrite this state during initialization
0152                     // unless the track is alive.
0153                     ti.geo.parent = tid;
0154                 }
0155
0156                 // Store the track initializer
0157                 data.initializers[ItemId<TrackInitializer>{
0158                     counters.num_initializers - offset}] = ti;
0159
0160                 --offset;
0161             }
0162         }
0163     }
0164
0165     if (sim.track_id() == track_id && sim.status() == TrackStatus::killed)
0166     {
0167         // Track is no longer used as part of transport
0168         sim.status(TrackStatus::inactive);
0169     }
0170     CELER_ENSURE(sim.status() != TrackStatus::killed);
0171 }
0172
0173 //---------------------------------------------------------------------------//
0174 }  // namespace detail
0175 }  // namespace celeritas