corecel/math/Integrator.hh

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