Back to home page

EIC code displayed by LXR

 
 

    


Warning, file /include/celeritas/phys/PhysicsData.hh was not indexed or was modified since last indexation (in which case cross-reference links may be missing, inaccurate or erroneous).

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/phys/PhysicsData.hh
0006 //---------------------------------------------------------------------------//
0007 #pragma once
0008 
0009 #include "corecel/cont/Array.hh"
0010 #include "corecel/data/Collection.hh"
0011 #include "corecel/data/CollectionBuilder.hh"
0012 #include "corecel/data/StackAllocatorData.hh"
0013 #include "celeritas/Quantities.hh"
0014 #include "celeritas/Types.hh"
0015 #include "celeritas/em/data/AtomicRelaxationData.hh"
0016 #include "celeritas/em/data/EPlusGGData.hh"
0017 #include "celeritas/em/data/LivermorePEData.hh"
0018 #include "celeritas/grid/XsGridData.hh"
0019 #include "celeritas/neutron/data/NeutronElasticData.hh"
0020 
0021 #include "Interaction.hh"
0022 #include "Secondary.hh"
0023 
0024 namespace celeritas
0025 {
0026 //---------------------------------------------------------------------------//
0027 // PARAMS
0028 //---------------------------------------------------------------------------//
0029 /*!
0030  * Set of value grids for all elements or materials.
0031  *
0032  * It is allowable for this to be "false" (i.e. no materials assigned)
0033  * indicating that the value table doesn't apply in the context -- for
0034  * example, an empty ValueTable macro_xs means that the process doesn't have a
0035  * discrete interaction.
0036  */
0037 template<class GridId>
0038 struct ValueTable
0039 {
0040     ItemRange<GridId> grids;  //!< Value grid by element or material index
0041 
0042     //! True if assigned
0043     explicit CELER_FUNCTION operator bool() const { return !grids.empty(); }
0044 };
0045 
0046 using UniformTable = ValueTable<UniformGridId>;
0047 
0048 //---------------------------------------------------------------------------//
0049 /*!
0050  * Set of cross section CDF tables for a model.
0051  *
0052  * Each material has a set of value grids for its constituent elements; these
0053  * are used to sample an element from a material when required by a discrete
0054  * interaction. A null \c ValueTableId means the material only has a single
0055  * element, so no cross sections need to be stored. An empty \c ModelCdfTable
0056  * means no element selection is required for the model.
0057  */
0058 struct ModelCdfTable
0059 {
0060     ItemRange<UniformTable> tables;  //!< Value table by material
0061 
0062     //! True if assigned
0063     explicit CELER_FUNCTION operator bool() const { return !tables.empty(); }
0064 };
0065 
0066 //---------------------------------------------------------------------------//
0067 /*!
0068  * Energy-dependent model IDs for a single process and particle type.
0069  *
0070  * For a given particle type, a single process should be divided into multiple
0071  * models as a function of energy. The \c ModelGroup represents this with an
0072  * energy grid, and each cell of the grid corresponding to a particular
0073  * \c ParticleModelId.
0074  */
0075 struct ModelGroup
0076 {
0077     using Energy = units::MevEnergy;
0078 
0079     ItemRange<real_type> energy;  //!< Energy grid bounds [MeV]
0080     ItemRange<ParticleModelId> model;  //!< Corresponding models
0081 
0082     //! True if assigned
0083     explicit CELER_FUNCTION operator bool() const
0084     {
0085         return (energy.size() >= 2) && (model.size() + 1 == energy.size());
0086     }
0087 };
0088 
0089 //---------------------------------------------------------------------------//
0090 /*!
0091  * Particle-process that uses MC integration to sample interaction length.
0092  *
0093  * This is needed for the integral approach for correctly sampling the discrete
0094  * interaction length after a particle loses energy along a step. An \c
0095  * IntegralXsProcess is stored for each particle-process. This will be "false"
0096  * (i.e. no energy_max assigned) if the particle associated with the process
0097  * does not have energy loss processes or if \c use_integral_xs is false.
0098  */
0099 struct IntegralXsProcess
0100 {
0101     ItemRange<real_type> energy_max_xs;  //!< Energy of the largest xs [mat]
0102 
0103     //! True if assigned
0104     explicit CELER_FUNCTION operator bool() const
0105     {
0106         return !energy_max_xs.empty();
0107     }
0108 };
0109 
0110 //---------------------------------------------------------------------------//
0111 /*!
0112  * Processes for a single particle type.
0113  *
0114  * Each index should be accessed with type ParticleProcessId. \c macro_xs
0115  * stores the cross section tables for each process, while \c energy_loss and
0116  * \c range are the process-integrated dE/dx and range for the particle.  \c
0117  * integral_xs will only be assigned if the integral approach is used and the
0118  * particle has continuous-discrete processes.
0119  *
0120  * \todo If it's possible for a particle to have multiple at-rest processes, \c
0121  * at_rest should be the process with the smallest lifetime. This is used in \c
0122  * select_discrete_interaction to choose the at-rest process with the smallest
0123  * time to interaction if the particle is stopped and has a process that
0124  * applies at rest.
0125  */
0126 struct ProcessGroup
0127 {
0128     ItemRange<ProcessId> processes;  //!< Processes that apply [ppid]
0129     ItemRange<ModelGroup> models;  //!< Model applicability [ppid]
0130     ItemRange<IntegralXsProcess> integral_xs;  //!< [ppid]
0131     ItemRange<ValueTable<XsGridId>> macro_xs;  //!< [ppid]
0132     UniformTable energy_loss;  //!< Process-integrated energy loss
0133     UniformTable range;  //!< Process-integrated range
0134     UniformTable inverse_range;  //!< Inverse process-integrated range
0135     ParticleProcessId at_rest;  //!< ID of the particle's at-rest process
0136 
0137     //! True if assigned and valid
0138     explicit CELER_FUNCTION operator bool() const
0139     {
0140         return !processes.empty() && models.size() == processes.size();
0141     }
0142 
0143     //! Number of processes that apply
0144     CELER_FUNCTION ParticleProcessId::size_type size() const
0145     {
0146         return processes.size();
0147     }
0148 };
0149 
0150 //---------------------------------------------------------------------------//
0151 /*!
0152  * IDs for models that do on-the-fly cross section calculation.
0153  */
0154 struct HardwiredIds
0155 {
0156     ProcessId annihilation;
0157     ModelId eplusgg;
0158 
0159     ProcessId photoelectric;
0160     ModelId livermore_pe;
0161 
0162     ProcessId neutron_elastic;
0163     ModelId chips;
0164 };
0165 
0166 //---------------------------------------------------------------------------//
0167 /*!
0168  * Model data for special hardwired cases (on-the-fly xs calculations).
0169  */
0170 template<Ownership W, MemSpace M>
0171 struct HardwiredModels
0172 {
0173     // Process and model IDs
0174     HardwiredIds ids;
0175 
0176     // Model data
0177     EPlusGGData eplusgg;
0178     LivermorePEData<W, M> livermore_pe;
0179     AtomicRelaxParamsData<W, M> relaxation;
0180     NeutronElasticData<W, M> chips;
0181 
0182     //! Assign from another set of data
0183     template<Ownership W2, MemSpace M2>
0184     HardwiredModels& operator=(HardwiredModels<W2, M2> const& other)
0185     {
0186         // Don't assign the references to model data
0187         ids = other.ids;
0188         eplusgg = other.eplusgg;
0189 
0190         return *this;
0191     }
0192 };
0193 
0194 //---------------------------------------------------------------------------//
0195 /*!
0196  * User-configurable particle-dependent physics constants.
0197  *
0198  * These scalar quantities can have different values for electrons/positrons
0199  * and muons/hadrons. They are described in \c PhysicsParams .
0200  */
0201 struct ParticleScalars
0202 {
0203     using Energy = units::MevEnergy;
0204 
0205     // Energy loss/range options
0206     real_type min_range{};  //!< rho [len]
0207     real_type max_step_over_range{};  //!< alpha [unitless]
0208     Energy lowest_energy{};  //!< Lowest kinetic energy
0209 
0210     // Multiple scattering options
0211     bool displaced{};  //!< Whether lateral displacement is enabled
0212     real_type range_factor{};
0213     MscStepLimitAlgorithm step_limit_algorithm{MscStepLimitAlgorithm::size_};
0214 
0215     //! True if assigned
0216     explicit CELER_FUNCTION operator bool() const
0217     {
0218         return min_range > 0 && max_step_over_range > 0
0219                && lowest_energy > zero_quantity() && range_factor > 0
0220                && range_factor < 1
0221                && step_limit_algorithm != MscStepLimitAlgorithm::size_;
0222     }
0223 };
0224 
0225 //---------------------------------------------------------------------------//
0226 /*!
0227  * Scalar (no template needed) quantities used by physics.
0228  *
0229  * The user-configurable constants and multiple scattering options are
0230  * described in \c PhysicsParams .
0231  *
0232  * The \c model_to_action value corresponds to the \c ActionId for the first \c
0233  * ModelId . Additionally it implies (by construction in physics_params) the
0234  * action IDs of several other physics actions.
0235  */
0236 struct PhysicsParamsScalars
0237 {
0238     using Energy = units::MevEnergy;
0239 
0240     //! Highest number of processes for any particle type
0241     ProcessId::size_type max_particle_processes{};
0242     //! Offset to create an ActionId from a ModelId
0243     ActionId::size_type model_to_action{};
0244     //! Number of physics models
0245     ModelId::size_type num_models{};
0246 
0247     // User-configurable constants
0248     real_type min_eprime_over_e{};  //!< xi [unitless]
0249     real_type linear_loss_limit{};  //!< For scaled range calculation
0250     real_type fixed_step_limiter{};  //!< Global charged step size limit [len]
0251 
0252     // User-configurable multiple scattering options
0253     real_type lambda_limit{};  //!< lambda limit
0254     real_type safety_factor{};  //!< safety factor
0255 
0256     // Particle-dependent user-configurable constants
0257     ParticleScalars light;
0258     ParticleScalars heavy;
0259 
0260     real_type secondary_stack_factor = 3;  //!< Secondary storage per state
0261                                            //!< size
0262     // When fixed step limiter is used, this is the corresponding action ID
0263     ActionId fixed_step_action{};
0264 
0265     //! True if assigned
0266     explicit CELER_FUNCTION operator bool() const
0267     {
0268         return max_particle_processes > 0 && model_to_action >= 4
0269                && num_models > 0 && min_eprime_over_e > 0
0270                && linear_loss_limit > 0 && secondary_stack_factor > 0
0271                && ((fixed_step_limiter > 0)
0272                    == static_cast<bool>(fixed_step_action))
0273                && lambda_limit > 0 && safety_factor >= 0.1 && light && heavy;
0274     }
0275 
0276     //! Stop early due to MSC limitation
0277     CELER_FORCEINLINE_FUNCTION ActionId msc_action() const
0278     {
0279         return ActionId{model_to_action - 4};
0280     }
0281 
0282     //! Stop early due to range limitation
0283     CELER_FORCEINLINE_FUNCTION ActionId range_action() const
0284     {
0285         return ActionId{model_to_action - 3};
0286     }
0287 
0288     //! Undergo a discrete interaction
0289     CELER_FORCEINLINE_FUNCTION ActionId discrete_action() const
0290     {
0291         return ActionId{model_to_action - 2};
0292     }
0293 
0294     //! Indicate a discrete interaction was rejected by the integral method
0295     CELER_FORCEINLINE_FUNCTION ActionId integral_rejection_action() const
0296     {
0297         return ActionId{model_to_action - 1};
0298     }
0299 
0300     //! Indicate an interaction failed to allocate memory
0301     CELER_FORCEINLINE_FUNCTION ActionId failure_action() const
0302     {
0303         return ActionId{model_to_action + num_models};
0304     }
0305 };
0306 
0307 //---------------------------------------------------------------------------//
0308 /*!
0309  * Persistent shared physics data.
0310  *
0311  * This includes macroscopic cross section tables ordered by
0312  * [particle][process][material][energy] and process-integrated energy loss and
0313  * range tables ordered by [particle][material][energy].
0314  */
0315 template<Ownership W, MemSpace M>
0316 struct PhysicsParamsData
0317 {
0318     //// TYPES ////
0319 
0320     template<class T>
0321     using Items = Collection<T, W, M>;
0322     template<class T>
0323     using ParticleItems = Collection<T, W, M, ParticleId>;
0324     template<class T>
0325     using ParticleModelItems = Collection<T, W, M, ParticleModelId>;
0326 
0327     //// DATA ////
0328 
0329     // Non-templated data
0330     PhysicsParamsScalars scalars;
0331 
0332     // Models that calculate cross sections on the fly
0333     HardwiredModels<Ownership::const_reference, M> hardwired;
0334 
0335     // Grid and table storage
0336     Items<XsGridId> xs_grid_ids;
0337     Items<XsGridRecord> xs_grids;
0338     Items<ValueTable<XsGridId>> xs_tables;
0339     Items<UniformGridId> uniform_grid_ids;
0340     Items<UniformGridRecord> uniform_grids;
0341     Items<UniformTable> uniform_tables;
0342     ParticleModelItems<ModelCdfTable> model_cdf;
0343 
0344     // Process and model storage
0345     Items<ModelGroup> model_groups;
0346     Items<IntegralXsProcess> integral_xs;
0347     ParticleItems<ProcessGroup> process_groups;
0348     ParticleModelItems<ModelId> model_ids;
0349     Items<ParticleModelId> pmodel_ids;
0350     Items<ProcessId> process_ids;
0351 
0352     // Backend storage
0353     Items<real_type> reals;
0354 
0355     //// METHODS ////
0356 
0357     //! True if assigned
0358     explicit CELER_FUNCTION operator bool() const
0359     {
0360         return !process_groups.empty() && !model_ids.empty() && scalars;
0361     }
0362 
0363     //! Assign from another set of data
0364     template<Ownership W2, MemSpace M2>
0365     PhysicsParamsData& operator=(PhysicsParamsData<W2, M2> const& other)
0366     {
0367         CELER_EXPECT(other);
0368 
0369         scalars = other.scalars;
0370 
0371         hardwired = other.hardwired;
0372 
0373         xs_grids = other.xs_grids;
0374         xs_grid_ids = other.xs_grid_ids;
0375         xs_tables = other.xs_tables;
0376         uniform_grids = other.uniform_grids;
0377         uniform_grid_ids = other.uniform_grid_ids;
0378         uniform_tables = other.uniform_tables;
0379         model_cdf = other.model_cdf;
0380 
0381         model_groups = other.model_groups;
0382         integral_xs = other.integral_xs;
0383         process_groups = other.process_groups;
0384         model_ids = other.model_ids;
0385         pmodel_ids = other.pmodel_ids;
0386         process_ids = other.process_ids;
0387 
0388         reals = other.reals;
0389 
0390         return *this;
0391     }
0392 };
0393 
0394 //---------------------------------------------------------------------------//
0395 // STATE
0396 //---------------------------------------------------------------------------//
0397 /*!
0398  * Physics state data for a single track.
0399  *
0400  * State that's persistent across steps:
0401  * - Remaining number of mean free paths to the next discrete interaction
0402  *
0403  * State that is reset at every step:
0404  * - Current macroscopic cross section
0405  * - Within-step energy deposition
0406  * - Within-step energy loss range
0407  * - Secondaries emitted from an interaction
0408  * - Discrete process element selection
0409  */
0410 struct PhysicsTrackState
0411 {
0412     real_type interaction_mfp;  //!< Remaining MFP to interaction
0413 
0414     // TEMPORARY STATE
0415     real_type macro_xs;  //!< Total cross section for discrete interactions
0416     real_type energy_deposition;  //!< Local energy deposition in a step [MeV]
0417     real_type dedx_range;  //!< Local energy loss range [len]
0418     MscRange msc_range;  //!< Range properties for multiple scattering
0419     Span<Secondary> secondaries;  //!< Emitted secondaries
0420     ElementComponentId element;  //!< Element sampled for interaction
0421 };
0422 
0423 //---------------------------------------------------------------------------//
0424 /*!
0425  * Initialize a physics track state.
0426  *
0427  * Currently no data is required at initialization -- it all must be evaluated
0428  * by the physics kernels itself.
0429  */
0430 struct PhysicsTrackInitializer
0431 {
0432 };
0433 
0434 //---------------------------------------------------------------------------//
0435 /*!
0436  * Dynamic physics (models, processes) state data.
0437  *
0438  * The "xs scratch space" is a 2D array of reals, indexed with
0439  * [track_id][el_component_id], where the fast-moving dimension has the
0440  * greatest number of element components of any material in the problem. This
0441  * can be used for the physics to calculate microscopic cross sections.
0442  */
0443 template<Ownership W, MemSpace M>
0444 struct PhysicsStateData
0445 {
0446     //// TYPES ////
0447 
0448     template<class T>
0449     using StateItems = celeritas::StateCollection<T, W, M>;
0450     template<class T>
0451     using Items = celeritas::Collection<T, W, M>;
0452 
0453     //// DATA ////
0454 
0455     StateItems<PhysicsTrackState> state;  //!< Track state [track]
0456     StateItems<MscStep> msc_step;  //!< Internal MSC data [track]
0457 
0458     Items<real_type> per_process_xs;  //!< XS [track][particle process]
0459 
0460     AtomicRelaxStateData<W, M> relaxation;  //!< Scratch data
0461     StackAllocatorData<Secondary, W, M> secondaries;  //!< Secondary stack
0462 
0463     //// METHODS ////
0464 
0465     //! True if assigned
0466     explicit CELER_FUNCTION operator bool() const
0467     {
0468         return !state.empty() && secondaries;
0469     }
0470 
0471     //! State size
0472     CELER_FUNCTION size_type size() const { return state.size(); }
0473 
0474     //! Assign from another set of states
0475     template<Ownership W2, MemSpace M2>
0476     PhysicsStateData& operator=(PhysicsStateData<W2, M2>& other)
0477     {
0478         CELER_EXPECT(other);
0479         state = other.state;
0480         msc_step = other.msc_step;
0481 
0482         per_process_xs = other.per_process_xs;
0483 
0484         relaxation = other.relaxation;
0485         secondaries = other.secondaries;
0486 
0487         return *this;
0488     }
0489 };
0490 
0491 //---------------------------------------------------------------------------//
0492 /*!
0493  * Resize the state in host code.
0494  */
0495 template<MemSpace M>
0496 inline void resize(PhysicsStateData<Ownership::value, M>* state,
0497                    HostCRef<PhysicsParamsData> const& params,
0498                    size_type size)
0499 {
0500     CELER_EXPECT(size > 0);
0501     CELER_EXPECT(params.scalars.max_particle_processes > 0);
0502     resize(&state->state, size);
0503     resize(&state->msc_step, size);
0504     resize(&state->per_process_xs,
0505            size * params.scalars.max_particle_processes);
0506     resize(&state->relaxation, params.hardwired.relaxation, size);
0507     resize(
0508         &state->secondaries,
0509         static_cast<size_type>(size * params.scalars.secondary_stack_factor));
0510 }
0511 
0512 //---------------------------------------------------------------------------//
0513 }  // namespace celeritas