celeritas/grid/XsCalculator.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/XsCalculator.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/math/Quantity.hh"
0015
0016 #include "SplineCalculator.hh"
0017 #include "UniformLogGridCalculator.hh"
0018 #include "XsGridData.hh"
0019
0020 namespace celeritas
0021 {
0022 //---------------------------------------------------------------------------//
0023 /*!
0024  * Find and interpolate scaled cross sections.
0025  *
0026  * This cross section calculator uses the same representation and interpolation
0027  * as Geant4's physics tables for EM physics:
0028  * - The energy grid is uniformly spaced in log(E),
0029  * - Values greater than or equal to an index i' are scaled by E and are stored
0030  *   on a separate energy grid also uniformly spaced in log(E) but not
0031  *   necessarily with the same spacing,
0032  * - Linear interpolation between energy points is used to calculate the final
0033  *   value, and
0034  * - If the energy is at or above the i' index, the final result is scaled by
0035  *   1/E.
0036  *
0037  * This scaling and interpolation exactly reproduces functions
0038  * \f$ f(E) \sim a E + b \f$ below the E' threshold and
0039  * \f$ f(E) \sim \frac{a'}{E} + b' \f$ above that threshold.
0040  *
0041  * Note that linear interpolation is applied with energy points, not log-energy
0042  * points.
0043  *
0044  * \code
0045     XsCalculator calc_xs(grid, params.reals);
0046     real_type xs = calc_xs(particle.energy());
0047    \endcode
0048  */
0049 class XsCalculator
0050 {
0051   public:
0052     //!@{
0053     //! \name Type aliases
0054     using Energy = RealQuantity<XsGridRecord::EnergyUnits>;
0055     using Values
0056         = Collection<real_type, Ownership::const_reference, MemSpace::native>;
0057     //!@}
0058
0059   public:
0060     // Construct from state-independent data
0061     inline CELER_FUNCTION
0062     XsCalculator(XsGridRecord const& grid, Values const& reals);
0063
0064     // Find and interpolate from the energy
0065     inline CELER_FUNCTION real_type operator()(Energy energy) const;
0066
0067     // Get the minimum energy
0068     inline CELER_FUNCTION Energy energy_min() const;
0069
0070     // Get the maximum energy
0071     inline CELER_FUNCTION Energy energy_max() const;
0072
0073   private:
0074     XsGridRecord const& data_;
0075     Values const& reals_;
0076
0077     inline CELER_FUNCTION bool use_scaled(Energy energy) const;
0078 };
0079
0080 //---------------------------------------------------------------------------//
0081 // INLINE DEFINITIONS
0082 //---------------------------------------------------------------------------//
0083 /*!
0084  * Construct from cross section data.
0085  */
0086 CELER_FUNCTION
0087 XsCalculator::XsCalculator(XsGridRecord const& grid, Values const& reals)
0088     : data_(grid), reals_(reals)
0089 {
0090     CELER_EXPECT(data_);
0091 }
0092
0093 //---------------------------------------------------------------------------//
0094 /*!
0095  * Calculate the cross section using linear or spline interpolation.
0096  */
0097 CELER_FUNCTION real_type XsCalculator::operator()(Energy energy) const
0098 {
0099     bool use_scaled = this->use_scaled(energy);
0100     auto const& grid = use_scaled ? data_.upper : data_.lower;
0101
0102     real_type result;
0103     if (grid.spline_order == 1)
0104     {
0105         // Linear or cubic spline interpolation
0106         result = UniformLogGridCalculator(grid, reals_)(energy);
0107     }
0108     else
0109     {
0110         // Spline interpolation without continuous derivatives
0111         result = SplineCalculator(grid, reals_)(energy);
0112     }
0113     if (use_scaled)
0114     {
0115         result /= value_as<Energy>(energy);
0116     }
0117     return result;
0118 }
0119
0120 //---------------------------------------------------------------------------//
0121 /*!
0122  * Get the minimum energy.
0123  */
0124 CELER_FUNCTION auto XsCalculator::energy_min() const -> Energy
0125 {
0126     return Energy(std::exp(data_.lower ? data_.lower.grid.front
0127                                        : data_.upper.grid.front));
0128 }
0129
0130 //---------------------------------------------------------------------------//
0131 /*!
0132  * Get the maximum energy.
0133  */
0134 CELER_FUNCTION auto XsCalculator::energy_max() const -> Energy
0135 {
0136     return Energy(
0137         std::exp(data_.upper ? data_.upper.grid.back : data_.lower.grid.back));
0138 }
0139
0140 //---------------------------------------------------------------------------//
0141 /*!
0142  * Whether to use the scaled cross section grid.
0143  */
0144 CELER_FUNCTION bool XsCalculator::use_scaled(Energy energy) const
0145 {
0146     return !data_.lower
0147            || (data_.upper
0148                && std::log(value_as<Energy>(energy)) >= data_.upper.grid.front);
0149 }
0150
0151 //---------------------------------------------------------------------------//
0152 }  // namespace celeritas