celeritas/grid/RangeCalculator.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/grid/RangeCalculator.hh
0006 //---------------------------------------------------------------------------//
0007 #pragma once
0008
0009 #include <cmath>
0010
0011 #include "corecel/data/Collection.hh"
0012 #include "corecel/grid/Interpolator.hh"
0013 #include "corecel/grid/SplineInterpolator.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  * \todo Construct with \c UniformGridRecord
0036  */
0037 class RangeCalculator
0038 {
0039   public:
0040     //!@{
0041     //! \name Type aliases
0042     using Energy = RealQuantity<XsGridRecord::EnergyUnits>;
0043     using Values
0044         = Collection<real_type, Ownership::const_reference, MemSpace::native>;
0045     //!@}
0046
0047   public:
0048     // Construct from state-independent data
0049     inline CELER_FUNCTION
0050     RangeCalculator(XsGridRecord const& grid, Values const& values);
0051
0052     // Find and interpolate from the energy
0053     inline CELER_FUNCTION real_type operator()(Energy energy) const;
0054
0055   private:
0056     UniformGridRecord const& data_;
0057     Values const& reals_;
0058
0059     CELER_FORCEINLINE_FUNCTION real_type get(size_type index) const;
0060 };
0061
0062 //---------------------------------------------------------------------------//
0063 // INLINE DEFINITIONS
0064 //---------------------------------------------------------------------------//
0065 /*!
0066  * Construct from cross section data.
0067  *
0068  * Range tables should be uniform in energy, without extra scaling.
0069  */
0070 CELER_FUNCTION
0071 RangeCalculator::RangeCalculator(XsGridRecord const& grid, Values const& values)
0072     : data_(grid.lower), reals_(values)
0073 {
0074     CELER_EXPECT(data_);
0075     CELER_EXPECT(!grid.upper);
0076 }
0077
0078 //---------------------------------------------------------------------------//
0079 /*!
0080  * Calculate the range.
0081  */
0082 CELER_FUNCTION real_type RangeCalculator::operator()(Energy energy) const
0083 {
0084     CELER_ASSERT(energy > zero_quantity());
0085     UniformGrid loge_grid(data_.grid);
0086     real_type const loge = std::log(energy.value());
0087
0088     if (loge <= loge_grid.front())
0089     {
0090         real_type result = this->get(0);
0091         // Scale by sqrt(E/Emin) = exp(.5 (log E - log Emin))
0092         result *= std::exp(real_type(.5) * (loge - loge_grid.front()));
0093         return result;
0094     }
0095     else if (loge >= loge_grid.back())
0096     {
0097         // Clip to highest range value
0098         return this->get(loge_grid.size() - 1);
0099     }
0100
0101     // Locate the energy bin
0102     auto idx = loge_grid.find(loge);
0103     CELER_ASSERT(idx + 1 < loge_grid.size());
0104
0105     real_type result;
0106     if (data_.derivative.empty())
0107     {
0108         // Interpolate *linearly* on energy
0109         result = LinearInterpolator<real_type>(
0110             {std::exp(loge_grid[idx]), this->get(idx)},
0111             {std::exp(loge_grid[idx + 1]), this->get(idx + 1)})(
0112             value_as<Energy>(energy));
0113     }
0114     else
0115     {
0116         // Use cubic spline interpolation
0117         real_type lower_deriv = reals_[data_.derivative[idx]];
0118         real_type upper_deriv = reals_[data_.derivative[idx + 1]];
0119
0120         result = SplineInterpolator<real_type>(
0121             {std::exp(loge_grid[idx]), this->get(idx), lower_deriv},
0122             {std::exp(loge_grid[idx + 1]), this->get(idx + 1), upper_deriv})(
0123             value_as<Energy>(energy));
0124     }
0125     return result;
0126 }
0127
0128 //---------------------------------------------------------------------------//
0129 /*!
0130  * Get the raw range data at a particular index.
0131  */
0132 CELER_FUNCTION real_type RangeCalculator::get(size_type index) const
0133 {
0134     CELER_EXPECT(index < reals_.size());
0135     return reals_[data_.value[index]];
0136 }
0137
0138 //---------------------------------------------------------------------------//
0139 }  // namespace celeritas