corecel/math/BisectionRootFinder.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/BisectionRootFinder.hh
0006 //---------------------------------------------------------------------------//
0007 #pragma once
0008
0009 #include <cmath>
0010 #include <type_traits>
0011
0012 #include "corecel/Config.hh"
0013
0014 #include "corecel/Assert.hh"
0015 #include "corecel/Macros.hh"
0016 #include "corecel/Types.hh"
0017
0018 #include "Algorithms.hh"
0019
0020 namespace celeritas
0021 {
0022 //---------------------------------------------------------------------------//
0023 /*!
0024  * Find a local root using the bisection algorithm.
0025  *
0026  * Using a \em left and \em right bound a Bisection approximates the \em
0027  * root as: \f[ root = 0.5 * (left + right) \f]
0028  *
0029  * Then value of \em func at the root is calculated compared to values of
0030  * \em func at the bounds. The \em root is then used update the bounds based on
0031  * the sign of \em func(root) and whether it matches the sign of \em func(left)
0032  * or \em func(right) . Performing this update of the bounds allows for the
0033  * iteration on root, using the convergence criteria based on \em func(root)
0034  * proximity to 0.
0035  */
0036 template<class F>
0037 class BisectionRootFinder
0038 {
0039   public:
0040     // Construct with function to solve and solution tolerance
0041     inline CELER_FUNCTION BisectionRootFinder(F&& func, real_type tol);
0042
0043     // Solve for a root between two points
0044     inline CELER_FUNCTION real_type operator()(real_type left, real_type right);
0045
0046   private:
0047     F func_;
0048     real_type tol_;
0049
0050     // Maximum amount of iterations
0051     static constexpr int max_iters_ = 50;
0052 };
0053
0054 //---------------------------------------------------------------------------//
0055 // TEMPLATE DEDUCTION
0056 //---------------------------------------------------------------------------//
0057
0058 template<class F, class... Args>
0059 CELER_CTAD_FUNCTION BisectionRootFinder(F&&, Args...)
0060     -> BisectionRootFinder<F>;
0061
0062 //---------------------------------------------------------------------------//
0063 // INLINE DEFINITIONS
0064 //---------------------------------------------------------------------------//
0065 /*!
0066  * Construct from function.
0067  */
0068 template<class F>
0069 CELER_FUNCTION BisectionRootFinder<F>::BisectionRootFinder(F&& func,
0070                                                            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     using namespace celeritas::literals;
0085
0086     // Initialize Iteration parameters
0087     real_type f_left = func_(left);
0088     real_type f_root = 1;
0089     real_type root = 0;
0090     int remaining_iters = max_iters_;
0091
0092     // Iterate on root
0093     do
0094     {
0095         // Estimate root and update value
0096         root = 0.5_r * (left + right);
0097         f_root = func_(root);
0098
0099         // Update the bound which produces the same sign as the root
0100         if (signum(f_left) == signum(f_root))
0101         {
0102             left = root;
0103             f_left = f_root;
0104         }
0105         else
0106         {
0107             right = root;
0108         }
0109     } while (std::fabs(f_root) > tol_ && --remaining_iters > 0);
0110
0111     CELER_ENSURE(remaining_iters > 0);
0112     return root;
0113 }
0114
0115 //---------------------------------------------------------------------------//
0116 }  // namespace celeritas