celeritas/random/ElementSelector.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/random/ElementSelector.hh
0006 //---------------------------------------------------------------------------//
0007 #pragma once
0008
0009 #include "corecel/Assert.hh"
0010 #include "corecel/Macros.hh"
0011 #include "corecel/Types.hh"
0012 #include "corecel/cont/Range.hh"
0013 #include "corecel/cont/Span.hh"
0014 #include "corecel/random/distribution/GenerateCanonical.hh"
0015 #include "celeritas/Types.hh"
0016 #include "celeritas/mat/MaterialView.hh"
0017
0018 namespace celeritas
0019 {
0020 //---------------------------------------------------------------------------//
0021 /*!
0022  * Make a weighted random selection of an element.
0023  *
0024  * The element chooser is for selecting an elemental component (atom) of a
0025  * material based on the microscopic cross section and the abundance fraction
0026  * of the element in the material.
0027  *
0028  * On construction, the element chooser uses the provided arguments to
0029  * precalculate all the microscopic cross sections in the given storage space.
0030  * The given function `calc_micro_xs` must accept a `ElementId` and return a
0031  * `real_type`, a non-negative microscopic cross section.
0032  *
0033  * The element chooser does \em not calculate macroscopic cross sections
0034  * because they're multiplied by fraction, not number density, and we only
0035  * care about the fractional abundances and cross section weighting.
0036  *
0037  * \code
0038     ElementSelector select_element(mat, calc_micro, storage);
0039     real_type total_macro_xs
0040         = select_element.material_micro_xs() * mat.number_density();
0041     ElementComponentId id = select_element(rng);
0042     real_type selected_micro_xs
0043         = select_element.elemental_micro_xs()[el.get()];
0044     ElementView el = mat.element_record(id);
0045     // use el.Z(), etc.
0046    \endcode
0047  *
0048  * Note that the units of the calculated microscopic cross section will be
0049  * identical to the units returned by `calc_micro_xs`. The macroscopic cross
0050  * section units (micro times \c mat.number_density() ) will be 1/len if and
0051  * only if calc_micro units are len^2.
0052  *
0053  * \todo Refactor to use Selector.
0054  */
0055 class ElementSelector
0056 {
0057   public:
0058     //!@{
0059     //! \name Type aliases
0060     using SpanReal = Span<real_type>;
0061     using SpanConstReal = LdgSpan<real_type const>;
0062     //!@}
0063
0064   public:
0065     // Construct with material, xs calculator, and storage.
0066     template<class MicroXsCalc>
0067     inline CELER_FUNCTION ElementSelector(MaterialView const& material,
0068                                           MicroXsCalc&& calc_micro_xs,
0069                                           SpanReal micro_xs_storage);
0070
0071     // Sample with the given RNG
0072     template<class Engine>
0073     inline CELER_FUNCTION ElementComponentId operator()(Engine& rng) const;
0074
0075     //! Weighted material microscopic cross section
0076     CELER_FUNCTION real_type material_micro_xs() const { return material_xs_; }
0077
0078     // Individual (unweighted) microscopic cross sections
0079     inline CELER_FUNCTION SpanConstReal elemental_micro_xs() const;
0080
0081   private:
0082     Span<MatElementComponent const> elements_;
0083     real_type material_xs_{0};
0084     real_type* elemental_xs_;
0085 };
0086
0087 //---------------------------------------------------------------------------//
0088 // INLINE DEFINITIONS
0089 //---------------------------------------------------------------------------//
0090 /*!
0091  * Construct with material, xs calculator, and storage.
0092  */
0093 template<class MicroXsCalc>
0094 CELER_FUNCTION ElementSelector::ElementSelector(MaterialView const& material,
0095                                                 MicroXsCalc&& calc_micro_xs,
0096                                                 SpanReal storage)
0097     : elements_(material.elements()), elemental_xs_(storage.data())
0098 {
0099     CELER_EXPECT(!elements_.empty());
0100     CELER_EXPECT(storage.size() >= material.num_elements());
0101     for (auto i : range<size_type>(elements_.size()))
0102     {
0103         real_type const micro_xs
0104             = (calc_micro_xs(elements_[i].element)).value();
0105         CELER_ASSERT(micro_xs >= 0);
0106         real_type const frac = elements_[i].fraction;
0107
0108         elemental_xs_[i] = micro_xs;
0109         material_xs_ += micro_xs * frac;
0110     }
0111     CELER_ENSURE(material_xs_ >= 0);
0112 }
0113
0114 //---------------------------------------------------------------------------//
0115 /*!
0116  * Sample the element with the given RNG.
0117  *
0118  * To reduce register usage, this function starts with the cumulative sums and
0119  * counts backward.
0120  */
0121 template<class Engine>
0122 CELER_FUNCTION ElementComponentId ElementSelector::operator()(Engine& rng) const
0123 {
0124     real_type accum_xs = -material_xs_ * generate_canonical(rng);
0125     size_type i = 0;
0126     size_type imax = elements_.size() - 1;
0127     for (; i != imax; ++i)
0128     {
0129         accum_xs += elements_[i].fraction * elemental_xs_[i];
0130         if (accum_xs > 0)
0131             break;
0132     }
0133     return ElementComponentId{i};
0134 }
0135
0136 //---------------------------------------------------------------------------//
0137 /*!
0138  * Get individual microscopic cross sections.
0139  */
0140 CELER_FUNCTION auto ElementSelector::elemental_micro_xs() const -> SpanConstReal
0141 {
0142     return {elemental_xs_, elements_.size()};
0143 }
0144
0145 //---------------------------------------------------------------------------//
0146 }  // namespace celeritas