corecel/math/Integrator.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 corecel/math/Integrator.hh
0006 //---------------------------------------------------------------------------//
0007 #pragma once
0008
0009 #include <cmath>
0010
0011 #include "corecel/Assert.hh"
0012 #include "corecel/Types.hh"
0013 #include "corecel/cont/Range.hh"
0014
0015 #include "Algorithms.hh"
0016
0017 namespace celeritas
0018 {
0019 //---------------------------------------------------------------------------//
0020 //! Solver options
0021 struct IntegratorOptions
0022 {
0023     using DepthInt = short unsigned int;
0024
0025     real_type epsilon{1e-3};  //!< Convergence criterion
0026     DepthInt max_depth{20};  //!< Maximum number of outer iterations
0027
0028     //! Whether the options are valid
0029     explicit operator bool() const
0030     {
0031         return epsilon > 0 && max_depth > 0
0032                && static_cast<std::size_t>(max_depth) < 8 * sizeof(size_type);
0033     }
0034 };
0035
0036 //---------------------------------------------------------------------------//
0037 /*!
0038  * Perform numerical integration of a generic 1-D function.
0039  *
0040  * Currently this is a simple nonrecursive Newton-Coates integrator extracted
0041  * from NuclearZoneBuilder. It should be improved for robustness, accuracy, and
0042  * efficiency, probably by using a Gauss-Legendre quadrature.
0043  *
0044  * This class is to be used \em only during setup.
0045  */
0046 template<class F>
0047 class Integrator
0048 {
0049   public:
0050     //!@{
0051     //! \name Type aliases
0052     using argument_type = real_type;
0053     using result_type = real_type;
0054     using Options = IntegratorOptions;
0055     //!@}
0056
0057   public:
0058     // Construct with the function and options
0059     inline Integrator(F&& func, Options opts);
0060
0061     //! Construct with the function and default options
0062     explicit Integrator(F&& func) : Integrator{std::forward<F>(func), {}} {}
0063
0064     // Calculate the integral over the given integral
0065     inline result_type operator()(argument_type lo, argument_type hi);
0066
0067   private:
0068     F eval_;
0069     Options options_;
0070
0071     //! Whether convergence is achieved (simple soft equality)
0072     bool converged(real_type prev, real_type cur) const
0073     {
0074         CELER_EXPECT(!std::isinf(cur) && !std::isnan(cur));
0075         return std::fabs(cur - prev) <= options_.epsilon * std::fabs(cur);
0076     }
0077 };
0078
0079 //---------------------------------------------------------------------------//
0080 // TEMPLATE DEDUCTION
0081 //---------------------------------------------------------------------------//
0082
0083 template<class F, class... Args>
0084 Integrator(F&&, Args...) -> Integrator<F>;
0085
0086 //---------------------------------------------------------------------------//
0087 // INLINE DEFINITIONS
0088 //---------------------------------------------------------------------------//
0089 /*!
0090  * Construct with the function and options.
0091  */
0092 template<class F>
0093 Integrator<F>::Integrator(F&& func, Options opts)
0094     : eval_{std::forward<F>(func)}, options_{opts}
0095 {
0096     CELER_EXPECT(options_);
0097 }
0098
0099 //---------------------------------------------------------------------------//
0100 /*!
0101  * Calculate the integral over the given dx.
0102  */
0103 template<class F>
0104 auto Integrator<F>::operator()(argument_type lo, argument_type hi) -> result_type
0105 {
0106     CELER_EXPECT(lo < hi);
0107     constexpr real_type half{0.5};
0108
0109     size_type num_points = 1;
0110     real_type dx = (hi - lo);
0111
0112     // Initial estimate is the endpoints of the function
0113     real_type result = half * dx * (eval_(lo) + eval_(hi));
0114     real_type prev{};
0115
0116     auto remaining_trials = options_.max_depth;
0117     do
0118     {
0119         // Accumulate the sum of midpoints along the grid spacing
0120         real_type accum = 0;
0121         for (auto i : range(num_points))
0122         {
0123             real_type x = std::fma((half + static_cast<real_type>(i)), dx, lo);
0124             accum += eval_(x);
0125         }
0126
0127         // Average previous and new integrations, i.e. combining all the
0128         // existing and current grid points
0129         prev = result;
0130         result = half * (prev + accum * dx);
0131
0132         // Increase number of intervals (and decrease dx) for next iteration
0133         num_points *= 2;
0134         dx *= half;
0135     } while (!this->converged(prev, result) && --remaining_trials > 0);
0136
0137     return result;
0138 }
0139
0140 //---------------------------------------------------------------------------//
0141 }  // namespace celeritas