corecel/math/IllinoisRootFinder.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/IllinoisRootFinder.hh
0007 //---------------------------------------------------------------------------//
0008 #pragma once
0009
0010 #include <cmath>
0011 #include <type_traits>
0012
0013 #include "corecel/Config.hh"
0014
0015 #include "corecel/Assert.hh"
0016 #include "corecel/Macros.hh"
0017 #include "corecel/Types.hh"
0018 #include "corecel/math/Algorithms.hh"
0019
0020 namespace celeritas
0021 {
0022 //---------------------------------------------------------------------------//
0023 /*!
0024  * Perform Regula Falsi (see RegulaFalsi for more details) iterations given a
0025  * root function \em func and tolerance \em tol using the Illinois method.
0026  *
0027  * Illinois method modifies the standard approach by comparing the sign of
0028  * \em func(root) approximation in the current iteration with the previous
0029  * approximation. If both iterations are on the same side then the \em func at
0030  * the bound on the other side is halved.
0031  */
0032 template<class F>
0033 class IllinoisRootFinder
0034 {
0035   public:
0036     // Contruct with function to solve and solution tolerance
0037     inline CELER_FUNCTION IllinoisRootFinder(F&& func, real_type tol);
0038
0039     // Solve for a root between two points
0040     inline CELER_FUNCTION real_type operator()(real_type left, real_type right);
0041
0042   private:
0043     F func_;
0044     real_type tol_;
0045
0046     // Maximum amount of iterations
0047     static constexpr inline int max_iters_ = 50;
0048 };
0049
0050 //---------------------------------------------------------------------------//
0051 // TEMPLATE DEDUCTION
0052 //---------------------------------------------------------------------------//
0053
0054 template<class F, class... Args>
0055 CELER_FUNCTION IllinoisRootFinder(F&&, Args...) -> IllinoisRootFinder<F>;
0056
0057 //---------------------------------------------------------------------------//
0058 // INLINE DEFINITIONS
0059 //---------------------------------------------------------------------------//
0060 /*!
0061  * Construct from function.
0062  */
0063 template<class F>
0064 CELER_FUNCTION
0065 IllinoisRootFinder<F>::IllinoisRootFinder(F&& func, real_type tol)
0066     : func_{celeritas::forward<F>(func)}, tol_{tol}
0067 {
0068     CELER_EXPECT(tol_ > 0);
0069 }
0070
0071 //---------------------------------------------------------------------------//
0072 /*!
0073  * Solve for a root between the two points.
0074  */
0075 template<class F>
0076 CELER_FUNCTION real_type IllinoisRootFinder<F>::operator()(real_type left,
0077                                                            real_type right)
0078 {
0079     //! Enum defining side of aproximated root to true root
0080     enum class Side
0081     {
0082         left = -1,
0083         init = 0,
0084         right = 1
0085     };
0086
0087     // Initialize Iteration parameters
0088     real_type f_left = func_(left);
0089     real_type f_right = func_(right);
0090     real_type f_root = 1;
0091     real_type root = 0;
0092     Side side = Side::init;
0093     int remaining_iters = max_iters_;
0094
0095     // Iterate on root
0096     do
0097     {
0098         // Estimate root and update value
0099         root = (left * f_right - right * f_left) / (f_right - f_left);
0100         f_root = func_(root);
0101
0102         // Update the bound which produces the same sign as the root
0103         if (signum(f_left) == signum(f_root))
0104         {
0105             left = root;
0106             f_left = f_root;
0107             if (side == Side::left)
0108             {
0109                 f_right *= real_type(0.5);
0110             }
0111             side = Side::left;
0112         }
0113         else
0114         {
0115             right = root;
0116             f_right = f_root;
0117             if (side == Side::right)
0118             {
0119                 f_left *= real_type(0.5);
0120             }
0121             side = Side::right;
0122         }
0123     } while (std::fabs(f_root) > tol_ && --remaining_iters > 0);
0124
0125     CELER_ENSURE(remaining_iters > 0);
0126     return root;
0127 }
0128
0129 //---------------------------------------------------------------------------//
0130 }  // namespace celeritas