track/detail/ProcessSecondariesExecutor.hh

0001 //----------------------------------*-C++-*----------------------------------//
0002 // Copyright 2020-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/track/detail/ProcessSecondariesExecutor.hh
0007 //---------------------------------------------------------------------------//
0008 #pragma once
0009
0010 #include "corecel/Assert.hh"
0011 #include "corecel/Macros.hh"
0012 #include "corecel/cont/Span.hh"
0013 #include "corecel/math/Atomics.hh"
0014 #include "corecel/sys/ThreadId.hh"
0015 #include "celeritas/Quantities.hh"
0016 #include "celeritas/Types.hh"
0017 #include "celeritas/geo/GeoTrackView.hh"
0018 #include "celeritas/global/CoreTrackData.hh"
0019 #include "celeritas/phys/ParticleData.hh"
0020 #include "celeritas/phys/ParticleTrackView.hh"
0021 #include "celeritas/phys/PhysicsStepView.hh"
0022 #include "celeritas/phys/PhysicsTrackView.hh"
0023 #include "celeritas/phys/Secondary.hh"
0024
0025 #include "../CoreStateCounters.hh"
0026 #include "../SimTrackView.hh"
0027
0028 namespace celeritas
0029 {
0030 namespace detail
0031 {
0032 //---------------------------------------------------------------------------//
0033 /*!
0034  * Create track initializers from secondaries.
0035  */
0036 struct ProcessSecondariesExecutor
0037 {
0038     //// TYPES ////
0039
0040     using ParamsPtr = CRefPtr<CoreParamsData, MemSpace::native>;
0041     using StatePtr = RefPtr<CoreStateData, MemSpace::native>;
0042
0043     //// DATA ////
0044
0045     ParamsPtr params;
0046     StatePtr state;
0047     CoreStateCounters counters;
0048
0049     //// FUNCTIONS ////
0050
0051     // Determine which tracks are alive and count secondaries
0052     inline CELER_FUNCTION void operator()(TrackSlotId tid) const;
0053
0054     CELER_FORCEINLINE_FUNCTION void operator()(ThreadId tid) const
0055     {
0056         // The grid size should be equal to the state size and no thread/slot
0057         // remapping should be performed
0058         return (*this)(TrackSlotId{tid.unchecked_get()});
0059     }
0060 };
0061
0062 //---------------------------------------------------------------------------//
0063 /*!
0064  * Create track initializers from secondaries.
0065  *
0066  * This kernel is executed with a grid size equal to the number of track
0067  * slots, so ThreadId should be equal to TrackSlotId. No remapping should be
0068  * done.
0069  */
0070 CELER_FUNCTION void
0071 ProcessSecondariesExecutor::operator()(TrackSlotId tid) const
0072 {
0073     CELER_EXPECT(tid < state->size());
0074
0075     SimTrackView sim(params->sim, state->sim, tid);
0076     if (sim.status() == TrackStatus::inactive)
0077     {
0078         // Do not create secondaries from stale data on inactive tracks
0079         return;
0080     }
0081
0082     // Offset in the vector of track initializers
0083     auto const& data = state->init;
0084     CELER_ASSERT(data.secondary_counts[tid] <= counters.num_secondaries);
0085     size_type offset = counters.num_secondaries - data.secondary_counts[tid];
0086
0087     // A new track was initialized from a secondary in the parent's track slot
0088     bool initialized = false;
0089
0090     // Save the parent ID since it will be overwritten if a secondary is
0091     // initialized in this slot
0092     TrackId const parent_id{sim.track_id()};
0093
0094     PhysicsStepView const phys_step(params->physics, state->physics, tid);
0095     for (auto const& secondary : phys_step.secondaries())
0096     {
0097         if (secondary)
0098         {
0099             CELER_ASSERT(secondary.energy > zero_quantity()
0100                          && is_soft_unit_vector(secondary.direction));
0101
0102             // Particles should not be making secondaries while crossing a
0103             // surface
0104             GeoTrackView geo(params->geometry, state->geometry, tid);
0105             CELER_ASSERT(!geo.is_on_boundary());
0106
0107             /*!
0108              * Increment the total number of tracks created for this event and
0109              * calculate the track ID of the secondary.
0110              *
0111              * \todo This is nondeterministic; we need to calculate the track
0112              * ID in a reproducible way.
0113              */
0114             CELER_ASSERT(sim.event_id() < data.track_counters.size());
0115             TrackId::size_type track_id = atomic_add(
0116                 &data.track_counters[sim.event_id()], size_type{1});
0117
0118             // Create a track initializer from the secondary
0119             TrackInitializer ti;
0120             ti.sim.track_id = TrackId{track_id};
0121             ti.sim.parent_id = parent_id;
0122             ti.sim.event_id = sim.event_id();
0123             ti.sim.time = sim.time();
0124             ti.geo.pos = geo.pos();
0125             ti.geo.dir = secondary.direction;
0126             ti.particle.particle_id = secondary.particle_id;
0127             ti.particle.energy = secondary.energy;
0128             CELER_ASSERT(ti);
0129
0130             if (!initialized && sim.status() != TrackStatus::alive
0131                 && params->init.track_order != TrackOrder::init_charge)
0132             {
0133                 /*!
0134                  * Skip in-place initialization when tracks are partitioned by
0135                  * charge to reduce the amount of mixing
0136                  *
0137                  * \todo Consider allowing this if the parent's charge is the
0138                  * same as the secondary's
0139                  */
0140
0141                 ParticleTrackView particle(
0142                     params->particles, state->particles, tid);
0143                 PhysicsTrackView phys(
0144                     params->physics, state->physics, {}, {}, tid);
0145
0146                 // The parent was killed, so initialize the first secondary in
0147                 // the parent's track slot. Keep the parent's geometry state
0148                 // but get the direction from the secondary. Reset the physics
0149                 // state so the multiple scattering range properties are
0150                 // cleared. The material state will be the same as the
0151                 // parent's.
0152                 sim = ti.sim;
0153                 geo = GeoTrackView::DetailedInitializer{geo, ti.geo.dir};
0154                 particle = ti.particle;
0155                 phys = {};
0156                 initialized = true;
0157
0158                 /*!
0159                  * \todo make it easier to determine what states need to be
0160                  * reset: the physics MFP, for example, is OK to preserve
0161                  */
0162             }
0163             else
0164             {
0165                 // Store the track initializer
0166                 CELER_ASSERT(offset > 0 && offset <= counters.num_initializers);
0167                 data.initializers[ItemId<TrackInitializer>{
0168                     counters.num_initializers - offset}]
0169                     = ti;
0170
0171                 if (offset <= data.parents.size()
0172                     && (params->init.track_order != TrackOrder::init_charge
0173                         || sim.status() == TrackStatus::alive))
0174                 {
0175                     // Store the thread ID of the secondary's parent if the
0176                     // secondary could be initialized in the next step. If the
0177                     // tracks are partitioned by charge we skip in-place
0178                     // initialization of the secondary, so the parent track
0179                     // must still be alive to ensure the state isn't
0180                     // overwritten
0181                     data.parents[TrackSlotId(data.parents.size() - offset)]
0182                         = tid;
0183                 }
0184                 --offset;
0185             }
0186         }
0187     }
0188
0189     if (!initialized && sim.status() == TrackStatus::killed)
0190     {
0191         // Track is no longer used as part of transport
0192         sim.status(TrackStatus::inactive);
0193     }
0194     CELER_ENSURE(sim.status() != TrackStatus::killed);
0195 }
0196
0197 //---------------------------------------------------------------------------//
0198 }  // namespace detail
0199 }  // namespace celeritas