celeritas/grid/UniformLogGridCalculator.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/UniformLogGridCalculator.hh
0006 //---------------------------------------------------------------------------//
0007 #pragma once
0008
0009 #include <cmath>
0010
0011 #include "corecel/grid/Interpolator.hh"
0012 #include "corecel/grid/SplineInterpolator.hh"
0013 #include "corecel/grid/UniformGrid.hh"
0014 #include "corecel/grid/UniformGridData.hh"
0015 #include "corecel/math/Quantity.hh"
0016 #include "celeritas/Quantities.hh"
0017
0018 namespace celeritas
0019 {
0020 //---------------------------------------------------------------------------//
0021 /*!
0022  * Find and interpolate values on a uniform log energy grid.
0023  *
0024  * Note that linear interpolation is applied with energy points, not log-energy
0025  * points.
0026  *
0027  * \code
0028     UniformLogGridCalculator calc(grid, params.reals);
0029     real_type y = calc(particle.energy());
0030    \endcode
0031  */
0032 class UniformLogGridCalculator
0033 {
0034   public:
0035     //!@{
0036     //! \name Type aliases
0037     using Energy = units::MevEnergy;
0038     using Values
0039         = Collection<real_type, Ownership::const_reference, MemSpace::native>;
0040     //!@}
0041
0042   public:
0043     // Construct from state-independent data
0044     inline CELER_FUNCTION
0045     UniformLogGridCalculator(UniformGridRecord const& grid,
0046                              Values const& reals);
0047
0048     // Find and interpolate from the energy
0049     inline CELER_FUNCTION real_type operator()(Energy energy) const;
0050
0051     // Get the tabulated value at the given index
0052     inline CELER_FUNCTION real_type operator[](size_type index) const;
0053
0054     // Get the minimum energy
0055     CELER_FUNCTION Energy energy_min() const
0056     {
0057         return Energy(std::exp(loge_grid_.front()));
0058     }
0059
0060     // Get the maximum energy
0061     CELER_FUNCTION Energy energy_max() const
0062     {
0063         return Energy(std::exp(loge_grid_.back()));
0064     }
0065
0066   private:
0067     UniformGridRecord const& data_;
0068     Values const& reals_;
0069     UniformGrid loge_grid_;
0070 };
0071
0072 //---------------------------------------------------------------------------//
0073 // INLINE DEFINITIONS
0074 //---------------------------------------------------------------------------//
0075 /*!
0076  * Construct from uniform grid data.
0077  */
0078 CELER_FUNCTION
0079 UniformLogGridCalculator::UniformLogGridCalculator(UniformGridRecord const& grid,
0080                                                    Values const& reals)
0081     : data_(grid), reals_(reals), loge_grid_(data_.grid)
0082 {
0083     CELER_EXPECT(data_);
0084 }
0085
0086 //---------------------------------------------------------------------------//
0087 /*!
0088  * Interpolate the value at the given energy.
0089  */
0090 CELER_FUNCTION real_type UniformLogGridCalculator::operator()(Energy energy) const
0091 {
0092     real_type const loge = std::log(value_as<Energy>(energy));
0093
0094     // Snap out-of-bounds values to closest grid points
0095     if (loge <= loge_grid_.front())
0096     {
0097         return (*this)[0];
0098     }
0099     if (loge >= loge_grid_.back())
0100     {
0101         return (*this)[loge_grid_.size() - 1];
0102     }
0103
0104     // Locate the energy bin
0105     size_type lower_idx = loge_grid_.find(loge);
0106     CELER_ASSERT(lower_idx + 1 < loge_grid_.size());
0107
0108     real_type lower_energy = std::exp(loge_grid_[lower_idx]);
0109     real_type upper_energy = std::exp(loge_grid_[lower_idx + 1]);
0110
0111     real_type result;
0112     if (data_.derivative.empty())
0113     {
0114         // Interpolate *linearly* on energy
0115         result = LinearInterpolator<real_type>(
0116             {lower_energy, (*this)[lower_idx]},
0117             {upper_energy, (*this)[lower_idx + 1]})(value_as<Energy>(energy));
0118     }
0119     else
0120     {
0121         // Use cubic spline interpolation
0122         real_type lower_deriv = reals_[data_.derivative[lower_idx]];
0123         real_type upper_deriv = reals_[data_.derivative[lower_idx + 1]];
0124
0125         result = SplineInterpolator<real_type>(
0126             {lower_energy, (*this)[lower_idx], lower_deriv},
0127             {upper_energy, (*this)[lower_idx + 1], upper_deriv})(
0128             value_as<Energy>(energy));
0129     }
0130     return result;
0131 }
0132
0133 //---------------------------------------------------------------------------//
0134 /*!
0135  * Get the tabulated value at the given index.
0136  */
0137 CELER_FUNCTION real_type UniformLogGridCalculator::operator[](size_type index) const
0138 {
0139     CELER_EXPECT(index < data_.value.size());
0140     return reals_[data_.value[index]];
0141 }
0142
0143 //---------------------------------------------------------------------------//
0144 }  // namespace celeritas