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