celeritas/grid/RangeCalculator.hh

0001 //----------------------------------*-C++-*----------------------------------//
0002 // Copyright 2021-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/grid/RangeCalculator.hh
0007 //---------------------------------------------------------------------------//
0008 #pragma once
0009
0010 #include <cmath>
0011
0012 #include "corecel/data/Collection.hh"
0013 #include "corecel/grid/Interpolator.hh"
0014 #include "corecel/grid/UniformGrid.hh"
0015 #include "corecel/math/Quantity.hh"
0016
0017 #include "XsGridData.hh"
0018
0019 namespace celeritas
0020 {
0021 //---------------------------------------------------------------------------//
0022 /*!
0023  * Find and interpolate range on a uniform log grid.
0024  *
0025  * \code
0026     RangeCalculator calc_range(xs_grid, xs_params.reals);
0027     real_type range = calc_range(particle);
0028    \endcode
0029  *
0030  * Below the minimum tabulated energy, the range is scaled:
0031  * \f[
0032     r = r_\mathrm{min} \sqrt{\frac{E}{E_\mathrm{min}}}
0033  * \f]
0034  */
0035 class RangeCalculator
0036 {
0037   public:
0038     //!@{
0039     //! \name Type aliases
0040     using Energy = Quantity<XsGridData::EnergyUnits>;
0041     using Values
0042         = Collection<real_type, Ownership::const_reference, MemSpace::native>;
0043     //!@}
0044
0045   public:
0046     // Construct from state-independent data
0047     inline CELER_FUNCTION
0048     RangeCalculator(XsGridData const& grid, Values const& values);
0049
0050     // Find and interpolate from the energy
0051     inline CELER_FUNCTION real_type operator()(Energy energy) const;
0052
0053   private:
0054     XsGridData const& data_;
0055     Values const& reals_;
0056
0057     CELER_FORCEINLINE_FUNCTION real_type get(size_type index) const;
0058 };
0059
0060 //---------------------------------------------------------------------------//
0061 // INLINE DEFINITIONS
0062 //---------------------------------------------------------------------------//
0063 /*!
0064  * Construct from cross section data.
0065  *
0066  * Range tables should be uniform in energy, without extra scaling.
0067  */
0068 CELER_FUNCTION
0069 RangeCalculator::RangeCalculator(XsGridData const& grid, Values const& values)
0070     : data_(grid), reals_(values)
0071 {
0072     CELER_EXPECT(data_);
0073     CELER_EXPECT(data_.prime_index == XsGridData::no_scaling());
0074 }
0075
0076 //---------------------------------------------------------------------------//
0077 /*!
0078  * Calculate the range.
0079  */
0080 CELER_FUNCTION real_type RangeCalculator::operator()(Energy energy) const
0081 {
0082     CELER_ASSERT(energy > zero_quantity());
0083     UniformGrid loge_grid(data_.log_energy);
0084     real_type const loge = std::log(energy.value());
0085
0086     if (loge <= loge_grid.front())
0087     {
0088         real_type result = this->get(0);
0089         // Scale by sqrt(E/Emin) = exp(.5 (log E - log Emin))
0090         result *= std::exp(real_type(.5) * (loge - loge_grid.front()));
0091         return result;
0092     }
0093     else if (loge >= loge_grid.back())
0094     {
0095         // Clip to highest range value
0096         return this->get(loge_grid.size() - 1);
0097     }
0098
0099     // Locate the energy bin
0100     auto idx = loge_grid.find(loge);
0101     CELER_ASSERT(idx + 1 < loge_grid.size());
0102
0103     // Interpolate *linearly* on energy
0104     LinearInterpolator<real_type> interpolate_xs(
0105         {std::exp(loge_grid[idx]), this->get(idx)},
0106         {std::exp(loge_grid[idx + 1]), this->get(idx + 1)});
0107     return interpolate_xs(energy.value());
0108 }
0109
0110 //---------------------------------------------------------------------------//
0111 /*!
0112  * Get the raw range data at a particular index.
0113  */
0114 CELER_FUNCTION real_type RangeCalculator::get(size_type index) const
0115 {
0116     CELER_EXPECT(index < reals_.size());
0117     return reals_[data_.value[index]];
0118 }
0119
0120 //---------------------------------------------------------------------------//
0121 }  // namespace celeritas