corecel/math/BisectionRootFinder.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/BisectionRootFinder.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 Bisection iterations given a root function \em func and
0025  * tolerance \em tol .
0026  *
0027  * Using a \em left and \em right bound a Bisection approximates the \em
0028  * root as: \f[ root = 0.5 * (left + right) \f]
0029  *
0030  * Then value of \em func at the root is calculated compared to values of
0031  * \em func at the bounds. The \em root is then used update the bounds based on
0032  * the sign of \em func(root) and whether it matches the sign of \em func(left)
0033  * or \em func(right) . Performing this update of the bounds allows for the
0034  * iteration on root, using the convergence criteria based on \em func(root)
0035  * proximity to 0.
0036  */
0037 template<class F>
0038 class BisectionRootFinder
0039 {
0040   public:
0041     // Contruct with function to solve and solution tolerance
0042     inline CELER_FUNCTION BisectionRootFinder(F&& func, real_type tol);
0043
0044     // Solve for a root between two points
0045     inline CELER_FUNCTION real_type operator()(real_type left, real_type right);
0046
0047   private:
0048     F func_;
0049     real_type tol_;
0050
0051     // Maximum amount of iterations
0052     static constexpr inline int max_iters_ = 50;
0053 };
0054
0055 //---------------------------------------------------------------------------//
0056 // TEMPLATE DEDUCTION
0057 //---------------------------------------------------------------------------//
0058
0059 template<class F, class... Args>
0060 CELER_FUNCTION BisectionRootFinder(F&&, Args...) -> BisectionRootFinder<F>;
0061
0062 //---------------------------------------------------------------------------//
0063 // INLINE DEFINITIONS
0064 //---------------------------------------------------------------------------//
0065 /*!
0066  * Construct from function.
0067  */
0068 template<class F>
0069 CELER_FUNCTION
0070 BisectionRootFinder<F>::BisectionRootFinder(F&& func, real_type tol)
0071     : func_{celeritas::forward<F>(func)}, tol_{tol}
0072 {
0073     CELER_EXPECT(tol_ > 0);
0074 }
0075
0076 //---------------------------------------------------------------------------//
0077 /*!
0078  * Solve for a root between the two points.
0079  */
0080 template<class F>
0081 CELER_FUNCTION real_type BisectionRootFinder<F>::operator()(real_type left,
0082                                                             real_type right)
0083 {
0084     // Initialize Iteration parameters
0085     real_type f_left = func_(left);
0086     real_type f_root = 1;
0087     real_type root = 0;
0088     int remaining_iters = max_iters_;
0089
0090     // Iterate on root
0091     do
0092     {
0093         // Estimate root and update value
0094         root = real_type(0.5) * (left + right);
0095         f_root = func_(root);
0096
0097         // Update the bound which produces the same sign as the root
0098         if (signum(f_left) == signum(f_root))
0099         {
0100             left = root;
0101             f_left = f_root;
0102         }
0103         else
0104         {
0105             right = root;
0106         }
0107     } while (std::fabs(f_root) > tol_ && --remaining_iters > 0);
0108
0109     CELER_ENSURE(remaining_iters > 0);
0110     return root;
0111 }
0112
0113 //---------------------------------------------------------------------------//
0114 }  // namespace celeritas