Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://mooseframework.inl.gov
       3             : //*
       4             : //* All rights reserved, see COPYRIGHT for full restrictions
       5             : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
       6             : //*
       7             : //* Licensed under LGPL 2.1, please see LICENSE for details
       8             : //* https://www.gnu.org/licenses/lgpl-2.1.html
       9             : 
      10             : #include "BrentsMethod.h"
      11             : #include "MooseError.h"
      12             : #include "Conversion.h"
      13             : 
      14             : namespace BrentsMethod
      15             : {
      16             : void
      17     8440156 : bracket(std::function<Real(Real)> const & f, Real & x1, Real & x2)
      18             : {
      19             :   Real f1, f2;
      20             :   // Factor to scale the interval by
      21             :   Real factor = 1.6;
      22             :   // Maximum number of times to attempt bracketing the interval
      23             :   unsigned int n = 50;
      24             :   // Small positive value to keep guess above
      25             :   Real eps = 1.0e-10;
      26             : 
      27             :   // If the initial guesses are identical
      28     8440156 :   if (x1 == x2)
      29           1 :     throw MooseException("Bad initial range (0) used in BrentsMethod::bracket");
      30             : 
      31     8440155 :   f1 = f(x1);
      32     8440155 :   f2 = f(x2);
      33             : 
      34     8440155 :   if (f1 * f2 > 0.0)
      35             :   {
      36             :     unsigned int iter = 0;
      37     8108275 :     std::stringstream debug_ss;
      38    16216601 :     while (f1 * f2 > 0.0)
      39             :     {
      40             : #ifdef DEBUG
      41             :       debug_ss << "  iteration " << iter << ": (x1,x2) = (" << x1 << "," << x2 << "), (f1,f2) = ("
      42             :                << f1 << "," << f2 << ")\n";
      43             : #endif
      44     8108327 :       if (std::abs(f1) < std::abs(f2))
      45             :       {
      46     8107825 :         x1 += factor * (x1 - x2);
      47     8107825 :         x1 = (x1 < eps ? eps : x1);
      48     8107825 :         f1 = f(x1);
      49             :       }
      50             :       else
      51             :       {
      52         502 :         x2 += factor * (x2 - x1);
      53         502 :         x2 = (x2 < eps ? eps : x2);
      54         502 :         f2 = f(x2);
      55             :       }
      56             :       /// Increment counter
      57     8108327 :       iter++;
      58     8108327 :       if (iter >= n)
      59           1 :         throw MooseException("No bracketing interval found by BrentsMethod::bracket after " +
      60           3 :                              Moose::stringify(n) + " iterations.\n" + debug_ss.str());
      61             :     }
      62     8108275 :   }
      63     8440154 : }
      64             : 
      65             : Real
      66     8440155 : root(std::function<Real(Real)> const & f, Real x1, Real x2, Real tol)
      67             : {
      68             :   Real a = x1, b = x2, c = x2, d = 0.0, e = 0.0, min1, min2;
      69     8440155 :   Real fa = f(a);
      70     8440155 :   Real fb = f(b);
      71             :   Real fc, p, q, r, s, tol1, xm = 0;
      72             :   unsigned int iter_max = 100;
      73             :   Real eps = 1.0e-12;
      74             : 
      75     8440155 :   if (fa * fb > 0.0)
      76           1 :     throw MooseException("Root must be bracketed in BrentsMethod::root");
      77             : 
      78             :   fc = fb;
      79     8440154 :   std::stringstream debug_ss;
      80    58900125 :   for (unsigned int i = 1; i <= iter_max; ++i)
      81             :   {
      82             : #ifdef DEBUG
      83             :     debug_ss << "  iteration " << i << ": dx = " << xm << ", x = " << b << ", f(x) = " << fb
      84             :              << "\n";
      85             : #endif
      86    58900125 :     if (fb * fc > 0.0)
      87             :     {
      88             :       // Rename a,b and c and adjust bounding interval d
      89             :       c = a;
      90             :       fc = fa;
      91    20052736 :       d = b - a;
      92             :       e = d;
      93             :     }
      94    58900125 :     if (std::abs(fc) < std::abs(fb))
      95             :     {
      96             :       a = b;
      97             :       b = c;
      98             :       c = a;
      99             :       fa = fb;
     100             :       fb = fc;
     101             :       fc = fa;
     102             :     }
     103             :     // Convergence check tolerance
     104    58900125 :     tol1 = 2.0 * eps * std::abs(b) + 0.5 * tol;
     105    58900125 :     xm = 0.5 * (c - b);
     106             : 
     107    58900125 :     if (std::abs(xm) <= tol1 || fb == 0.0)
     108             :       return b;
     109             : 
     110    50459971 :     if (std::abs(e) >= tol1 && std::abs(fa) > std::abs(fb))
     111             :     {
     112             :       // Attempt inverse quadratic interpolation
     113    50459904 :       s = fb / fa;
     114    50459904 :       if (a == c)
     115             :       {
     116    12922717 :         p = 2.0 * xm * s;
     117    12922717 :         q = 1.0 - s;
     118             :       }
     119             :       else
     120             :       {
     121    37537187 :         q = fa / fc;
     122    37537187 :         r = fb / fc;
     123    37537187 :         p = s * (2.0 * xm * q * (q - r) - (b - a) * (r - 1.0));
     124    37537187 :         q = (q - 1.0) * (r - 1.0) * (s - 1.0);
     125             :       }
     126             :       // Check whether in bounds
     127    50459904 :       if (p > 0.0)
     128    38814730 :         q = -q;
     129             :       p = std::abs(p);
     130    50459904 :       min1 = 3.0 * xm * q - std::abs(tol1 * q);
     131    50459904 :       min2 = std::abs(e * q);
     132             : 
     133   100628954 :       if (2.0 * p < (min1 < min2 ? min1 : min2))
     134             :       {
     135             :         // Accept interpolation
     136             :         e = d;
     137    50446478 :         d = p / q;
     138             :       }
     139             :       else
     140             :       {
     141             :         // Interpolation failed, use bisection
     142             :         d = xm;
     143             :         e = d;
     144             :       }
     145             :     }
     146             :     else
     147             :     {
     148             :       // Bounds decreasing too slowly, use bisection
     149             :       d = xm;
     150             :       e = d;
     151             :     }
     152             :     // Move last best guess to a
     153             :     a = b;
     154             :     fa = fb;
     155             :     // Evaluate new trial root
     156    50459971 :     if (std::abs(d) > tol1)
     157    43468724 :       b += d;
     158             :     else
     159             :     {
     160     6991247 :       Real sgn = (xm >= 0.0 ? std::abs(tol1) : -std::abs(tol1));
     161     6991247 :       b += sgn;
     162             :     }
     163             : 
     164    50459971 :     fb = f(b);
     165             :   }
     166             : 
     167           0 :   throw MooseException("Maximum number of iterations exceeded in BrentsMethod::root.\n" +
     168           0 :                        debug_ss.str());
     169             :   return 0.0; // Should never get here
     170     8440154 : }
     171             : } // namespace BrentsMethod