celeritas/grid/NonuniformGridCalculator.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/NonuniformGridCalculator.hh
0006 //---------------------------------------------------------------------------//
0007 #pragma once
0008
0009 #include <cmath>
0010
0011 #include "corecel/Macros.hh"
0012 #include "corecel/Types.hh"
0013 #include "corecel/data/Collection.hh"
0014 #include "corecel/grid/Interpolator.hh"
0015 #include "corecel/grid/NonuniformGrid.hh"
0016 #include "corecel/grid/NonuniformGridData.hh"
0017 #include "corecel/grid/SplineInterpolator.hh"
0018
0019 namespace celeritas
0020 {
0021 //---------------------------------------------------------------------------//
0022 /*!
0023  * Find and interpolate real numbers on a nonuniform grid.
0024  *
0025  * The end points of the grid are extrapolated outward as constant values.
0026  *
0027  * \todo Template on value type and/or units?
0028  */
0029 class NonuniformGridCalculator
0030 {
0031   public:
0032     //@{
0033     //! Type aliases
0034     using Values
0035         = Collection<real_type, Ownership::const_reference, MemSpace::native>;
0036     using Grid = NonuniformGrid<real_type>;
0037     //@}
0038
0039   public:
0040     // Construct by *inverting* a monotonicially increasing nonuniform grid
0041     static inline CELER_FUNCTION NonuniformGridCalculator
0042     from_inverse(NonuniformGridRecord const& grid, Values const& reals);
0043
0044     // Construct from grid data and backend storage
0045     inline CELER_FUNCTION
0046     NonuniformGridCalculator(NonuniformGridRecord const& grid,
0047                              Values const& reals);
0048
0049     // Find and interpolate the y value from the given x value
0050     inline CELER_FUNCTION real_type operator()(real_type x) const;
0051
0052     // Get the tabulated y value at a particular index.
0053     inline CELER_FUNCTION real_type operator[](size_type index) const;
0054
0055     // Get the tabulated x grid
0056     inline CELER_FUNCTION Grid const& grid() const;
0057
0058     // Make a calculator with x and y flipped
0059     inline CELER_FUNCTION NonuniformGridCalculator make_inverse() const;
0060
0061   private:
0062     //// TYPES ////
0063
0064     using RealIds = ItemRange<real_type>;
0065
0066     //// DATA ////
0067
0068     Values const& reals_;
0069     Grid x_grid_;
0070     RealIds y_offset_;
0071     RealIds deriv_offset_;
0072
0073     //// HELPER FUNCTIONS ////
0074
0075     // Private constructor implementation
0076     inline CELER_FUNCTION NonuniformGridCalculator(Values const& reals,
0077                                                    RealIds x_grid,
0078                                                    RealIds y_grid,
0079                                                    RealIds deriv);
0080 };
0081
0082 //---------------------------------------------------------------------------//
0083 // INLINE DEFINITIONS
0084 //---------------------------------------------------------------------------//
0085 /*!
0086  * Construct by \em inverting a monotonicially increasing nonuniform grid.
0087  */
0088 CELER_FUNCTION NonuniformGridCalculator NonuniformGridCalculator::from_inverse(
0089     NonuniformGridRecord const& grid, Values const& reals)
0090 {
0091     CELER_EXPECT(grid.derivative.empty());
0092     return NonuniformGridCalculator{reals, grid.value, grid.grid, {}};
0093 }
0094
0095 //---------------------------------------------------------------------------//
0096 /*!
0097  * Construct from grid data and backend storage.
0098  */
0099 CELER_FUNCTION
0100 NonuniformGridCalculator::NonuniformGridCalculator(
0101     NonuniformGridRecord const& grid, Values const& reals)
0102     : NonuniformGridCalculator{reals, grid.grid, grid.value, grid.derivative}
0103 {
0104     CELER_EXPECT(grid);
0105 }
0106
0107 //---------------------------------------------------------------------------//
0108 /*!
0109  * Calculate the y value at the given x value.
0110  */
0111 CELER_FUNCTION real_type NonuniformGridCalculator::operator()(real_type x) const
0112 {
0113     // Snap out-of-bounds values to closest grid points
0114     if (x <= x_grid_.front())
0115     {
0116         return (*this)[0];
0117     }
0118     if (x >= x_grid_.back())
0119     {
0120         return (*this)[x_grid_.size() - 1];
0121     }
0122
0123     // Locate the x bin
0124     size_type lower_idx = x_grid_.find(x);
0125     CELER_ASSERT(lower_idx + 1 < x_grid_.size());
0126
0127     real_type result;
0128     if (deriv_offset_.empty())
0129     {
0130         // Interpolate *linearly* on x using the bin data.
0131         result = LinearInterpolator<real_type>(
0132             {x_grid_[lower_idx], (*this)[lower_idx]},
0133             {x_grid_[lower_idx + 1], (*this)[lower_idx + 1]})(x);
0134     }
0135     else
0136     {
0137         // Use cubic spline interpolation
0138         real_type lower_deriv = reals_[deriv_offset_[lower_idx]];
0139         real_type upper_deriv = reals_[deriv_offset_[lower_idx + 1]];
0140
0141         result = SplineInterpolator<real_type>(
0142             {x_grid_[lower_idx], (*this)[lower_idx], lower_deriv},
0143             {x_grid_[lower_idx + 1], (*this)[lower_idx + 1], upper_deriv})(x);
0144     }
0145     return result;
0146 }
0147
0148 //---------------------------------------------------------------------------//
0149 /*!
0150  * Get the tabulated y value at a particular index.
0151  */
0152 CELER_FUNCTION real_type NonuniformGridCalculator::operator[](size_type index) const
0153 {
0154     CELER_EXPECT(index < y_offset_.size());
0155     return reals_[y_offset_[index]];
0156 }
0157
0158 //---------------------------------------------------------------------------//
0159 /*!
0160  * Get the tabulated x values.
0161  */
0162 CELER_FORCEINLINE_FUNCTION NonuniformGrid<real_type> const&
0163 NonuniformGridCalculator::grid() const
0164 {
0165     return x_grid_;
0166 }
0167
0168 //---------------------------------------------------------------------------//
0169 /*!
0170  * Make a calculator with x and y flipped.
0171  *
0172  * \pre The y values must be monotonic increasing.
0173  *
0174  * \note Spline interpolation will not be used on an inverse grid.
0175  */
0176 CELER_FUNCTION NonuniformGridCalculator
0177 NonuniformGridCalculator::make_inverse() const
0178 {
0179     CELER_EXPECT(deriv_offset_.empty());
0180     return NonuniformGridCalculator{reals_, y_offset_, x_grid_.offset(), {}};
0181 }
0182
0183 //---------------------------------------------------------------------------//
0184 // PRIVATE HELPERS
0185 //---------------------------------------------------------------------------//
0186 /*!
0187  * Construct from grid data and backend storage.
0188  */
0189 CELER_FUNCTION
0190 NonuniformGridCalculator::NonuniformGridCalculator(Values const& reals,
0191                                                    RealIds x_grid,
0192                                                    RealIds y_grid,
0193                                                    RealIds deriv)
0194     : reals_{reals}
0195     , x_grid_{x_grid, reals_}
0196     , y_offset_{y_grid}
0197     , deriv_offset_(deriv)
0198 {
0199     CELER_EXPECT(!x_grid.empty() && x_grid.size() == y_grid.size());
0200     CELER_EXPECT(*x_grid.end() <= reals_.size()
0201                  && *y_grid.end() <= reals_.size());
0202     CELER_EXPECT(deriv.empty() || deriv.size() == x_grid.size());
0203 }
0204
0205 //---------------------------------------------------------------------------//
0206 }  // namespace celeritas