celeritas/grid/InverseRangeCalculator.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/InverseRangeCalculator.hh
0006 //---------------------------------------------------------------------------//
0007 #pragma once
0008
0009 #include <cmath>
0010
0011 #include "corecel/Assert.hh"
0012 #include "corecel/data/Collection.hh"
0013 #include "corecel/grid/Interpolator.hh"
0014 #include "corecel/grid/NonuniformGrid.hh"
0015 #include "corecel/grid/SplineInterpolator.hh"
0016 #include "corecel/grid/UniformGrid.hh"
0017 #include "corecel/math/Algorithms.hh"
0018 #include "corecel/math/Quantity.hh"
0019
0020 #include "XsGridData.hh"
0021
0022 namespace celeritas
0023 {
0024 //---------------------------------------------------------------------------//
0025 /*!
0026  * Calculate the energy that would limit a particle to a particular range.
0027  *
0028  * This should provide the inverse of the result of \c RangeCalculator. The
0029  * given \c range is not allowed to be greater than the maximum range in the
0030  * physics data.
0031  *
0032  * The range must be monotonically increasing in energy, since it's defined as
0033  * the integral of the inverse of the stopping power (which is always
0034  * positive). For ranges shorter than the minimum energy in the table, the
0035  * resulting energy is scaled:
0036  * \f[
0037     E = E_\mathrm{min} \left( \frac{r}{r_\mathrm{min}} \right)^2
0038  * \f]
0039  * This scaling is the inverse of the off-the-end energy scaling in the
0040  * RangeCalculator.
0041  *
0042  * \todo Construct with \c UniformGridRecord
0043  */
0044 class InverseRangeCalculator
0045 {
0046   public:
0047     //!@{
0048     //! \name Type aliases
0049     using Energy = RealQuantity<XsGridRecord::EnergyUnits>;
0050     using Values
0051         = Collection<real_type, Ownership::const_reference, MemSpace::native>;
0052     //!@}
0053
0054   public:
0055     // Construct from state-independent data
0056     inline CELER_FUNCTION
0057     InverseRangeCalculator(XsGridRecord const& grid, Values const& values);
0058
0059     // Find and interpolate from the energy
0060     inline CELER_FUNCTION Energy operator()(real_type range) const;
0061
0062   private:
0063     UniformGrid log_energy_;
0064     NonuniformGrid<real_type> range_;
0065     Span<real_type const> deriv_;
0066 };
0067
0068 //---------------------------------------------------------------------------//
0069 // INLINE DEFINITIONS
0070 //---------------------------------------------------------------------------//
0071 /*!
0072  * Construct from range data.
0073  *
0074  * The range is expected to be monotonically increaing with energy.
0075  * Lower-energy particles have shorter ranges.
0076  */
0077 CELER_FUNCTION
0078 InverseRangeCalculator::InverseRangeCalculator(XsGridRecord const& grid,
0079                                                Values const& values)
0080     : log_energy_(grid.lower.grid)
0081     , range_(grid.lower.value, values)
0082     , deriv_(values[grid.lower.derivative])
0083 {
0084     CELER_EXPECT(range_.size() == log_energy_.size());
0085     CELER_EXPECT(!grid.upper);
0086 }
0087
0088 //---------------------------------------------------------------------------//
0089 /*!
0090  * Calculate the energy of a particle that has the given range.
0091  */
0092 CELER_FUNCTION auto
0093 InverseRangeCalculator::operator()(real_type range) const -> Energy
0094 {
0095     CELER_EXPECT(range >= 0 && range <= range_.back());
0096
0097     if (range < range_.front())
0098     {
0099         // Very short range:  this corresponds to "energy < emin" for range
0100         // calculation: range = r[0] * sqrt(E / E[0])
0101         return Energy{std::exp(log_energy_.front())
0102                       * ipow<2>(range / range_.front())};
0103     }
0104     // Range should *never* exceed the longest range (highest energy) since
0105     // that should have limited the step
0106     if (CELER_UNLIKELY(range >= range_.back()))
0107     {
0108         CELER_ASSERT(range == range_.back());
0109         return Energy{std::exp(log_energy_.back())};
0110     }
0111
0112     // Search for lower bin index
0113     auto idx = range_.find(range);
0114     CELER_ASSERT(idx + 1 < log_energy_.size());
0115
0116     real_type result;
0117     if (deriv_.empty())
0118     {
0119         // Interpolate: 'x' = range, y = log energy
0120         result = LinearInterpolator<real_type>(
0121             {range_[idx], std::exp(log_energy_[idx])},
0122             {range_[idx + 1], std::exp(log_energy_[idx + 1])})(range);
0123     }
0124     else
0125     {
0126         // Use cubic spline interpolation
0127         result = SplineInterpolator<real_type>(
0128             {range_[idx], std::exp(log_energy_[idx]), deriv_[idx]},
0129             {range_[idx + 1], std::exp(log_energy_[idx + 1]), deriv_[idx + 1]})(
0130             range);
0131     }
0132     return Energy{result};
0133 }
0134
0135 //---------------------------------------------------------------------------//
0136 }  // namespace celeritas