celeritas/grid/InverseRangeCalculator.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/InverseRangeCalculator.hh
0007 //---------------------------------------------------------------------------//
0008 #pragma once
0009
0010 #include <cmath>
0011
0012 #include "corecel/Assert.hh"
0013 #include "corecel/data/Collection.hh"
0014 #include "corecel/grid/Interpolator.hh"
0015 #include "corecel/grid/NonuniformGrid.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 class InverseRangeCalculator
0043 {
0044   public:
0045     //!@{
0046     //! \name Type aliases
0047     using Energy = Quantity<XsGridData::EnergyUnits>;
0048     using Values
0049         = Collection<real_type, Ownership::const_reference, MemSpace::native>;
0050     //!@}
0051
0052   public:
0053     // Construct from state-independent data
0054     inline CELER_FUNCTION
0055     InverseRangeCalculator(XsGridData const& grid, Values const& values);
0056
0057     // Find and interpolate from the energy
0058     inline CELER_FUNCTION Energy operator()(real_type range) const;
0059
0060   private:
0061     UniformGrid log_energy_;
0062     NonuniformGrid<real_type> range_;
0063 };
0064
0065 //---------------------------------------------------------------------------//
0066 // INLINE DEFINITIONS
0067 //---------------------------------------------------------------------------//
0068 /*!
0069  * Construct from range data.
0070  *
0071  * The range is expected to be monotonically increaing with energy.
0072  * Lower-energy particles have shorter ranges.
0073  */
0074 CELER_FUNCTION
0075 InverseRangeCalculator::InverseRangeCalculator(XsGridData const& grid,
0076                                                Values const& values)
0077     : log_energy_(grid.log_energy), range_(grid.value, values)
0078 {
0079     CELER_EXPECT(range_.size() == log_energy_.size());
0080 }
0081
0082 //---------------------------------------------------------------------------//
0083 /*!
0084  * Calculate the energy of a particle that has the given range.
0085  */
0086 CELER_FUNCTION auto
0087 InverseRangeCalculator::operator()(real_type range) const -> Energy
0088 {
0089     CELER_EXPECT(range >= 0 && range <= range_.back());
0090
0091     if (range < range_.front())
0092     {
0093         // Very short range:  this corresponds to "energy < emin" for range
0094         // calculation: range = r[0] * sqrt(E / E[0])
0095         return Energy{std::exp(log_energy_.front())
0096                       * ipow<2>(range / range_.front())};
0097     }
0098     // Range should *never* exceed the longest range (highest energy) since
0099     // that should have limited the step
0100     if (CELER_UNLIKELY(range >= range_.back()))
0101     {
0102         CELER_ASSERT(range == range_.back());
0103         return Energy{std::exp(log_energy_.back())};
0104     }
0105
0106     // Search for lower bin index
0107     auto idx = range_.find(range);
0108     CELER_ASSERT(idx + 1 < log_energy_.size());
0109
0110     // Interpolate: 'x' = range, y = log energy
0111     LinearInterpolator<real_type> interpolate_log_energy(
0112         {range_[idx], std::exp(log_energy_[idx])},
0113         {range_[idx + 1], std::exp(log_energy_[idx + 1])});
0114     auto loge = interpolate_log_energy(range);
0115     return Energy{loge};
0116 }
0117
0118 //---------------------------------------------------------------------------//
0119 }  // namespace celeritas