OilfieldToolsNet/opm-core
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Whitespace cleanup.

Browse Source
This commit is contained in:
Atgeirr Flø Rasmussen 2012-05-25 10:12:54 +02:00
parent 4714c30cb0
commit b8b66d6f35
1 changed files with 125 additions and 125 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

250
opm/core/utility/RootFinders.hpp
View File

@ -103,110 +103,110 @@ namespace Opm
public:
/// Implements a modified regula falsi method as described in
/// "Improved algorithms of Illinois-type for the numerical
/// solution of nonlinear equations"
/// by J. A. Ford.
/// Current variant is the 'Pegasus' method.
/// Implements a modified regula falsi method as described in
/// "Improved algorithms of Illinois-type for the numerical
/// solution of nonlinear equations"
/// by J. A. Ford.
/// Current variant is the 'Pegasus' method.
template <class Functor>
inline static double solve(const Functor& f,
inline static double solve(const Functor& f,
const double a,
const double b,
const int max_iter,
const double tolerance,
int& iterations_used)
{
using namespace std;
const double macheps = numeric_limits<double>::epsilon();
const double eps = tolerance + macheps*max(max(fabs(a), fabs(b)), 1.0);
{
using namespace std;
const double macheps = numeric_limits<double>::epsilon();
const double eps = tolerance + macheps*max(max(fabs(a), fabs(b)), 1.0);
double x0 = a;
double x1 = b;
double f0 = f(x0);
const double epsF = tolerance + macheps*max(fabs(f0), 1.0);
if (fabs(f0) < epsF) {
return x0;
}
double f1 = f(x1);
if (fabs(f1) < epsF) {
return x1;
}
if (f0*f1 > 0.0) {
double x0 = a;
double x1 = b;
double f0 = f(x0);
const double epsF = tolerance + macheps*max(fabs(f0), 1.0);
if (fabs(f0) < epsF) {
return x0;
}
double f1 = f(x1);
if (fabs(f1) < epsF) {
return x1;
}
if (f0*f1 > 0.0) {
return ErrorPolicy::handleBracketingFailure(a, b, f0, f1);
}
iterations_used = 0;
// In every iteraton, x1 is the last point computed,
// and x0 is the last point computed that makes it a bracket.
while (fabs(x1 - x0) >= 1e-9*eps) {
double xnew = regulaFalsiStep(x0, x1, f0, f1);
double fnew = f(xnew);
}
iterations_used = 0;
// In every iteraton, x1 is the last point computed,
// and x0 is the last point computed that makes it a bracket.
while (fabs(x1 - x0) >= 1e-9*eps) {
double xnew = regulaFalsiStep(x0, x1, f0, f1);
double fnew = f(xnew);
// cout << "xnew = " << xnew << " fnew = " << fnew << endl;
++iterations_used;
if (iterations_used > max_iter) {
++iterations_used;
if (iterations_used > max_iter) {
return ErrorPolicy::handleTooManyIterations(x0, x1, max_iter);
}
if (fabs(fnew) < epsF) {
return xnew;
}
// Now we must check which point we must replace.
if ((fnew > 0.0) == (f0 > 0.0)) {
// We must replace x0.
x0 = x1;
f0 = f1;
} else {
// We must replace x1, this is the case where
// the modification to regula falsi kicks in,
// by modifying f0.
// 1. The classic Illinois method
// const double gamma = 0.5;
// @afr: The next two methods do not work??!!?
// 2. The method called 'Method 3' in the paper.
// const double phi0 = f1/f0;
// const double phi1 = fnew/f1;
// const double gamma = 1.0 - phi1/(1.0 - phi0);
// 3. The method called 'Method 4' in the paper.
// const double phi0 = f1/f0;
// const double phi1 = fnew/f1;
// const double gamma = 1.0 - phi0 - phi1;
// cout << "phi0 = " << phi0 <<" phi1 = " << phi1 <<
// " gamma = " << gamma << endl;
// 4. The 'Pegasus' method
const double gamma = f1/(f1 + fnew);
f0 *= gamma;
}
x1 = xnew;
f1 = fnew;
}
return 0.5*(x0 + x1);
}
}
if (fabs(fnew) < epsF) {
return xnew;
}
// Now we must check which point we must replace.
if ((fnew > 0.0) == (f0 > 0.0)) {
// We must replace x0.
x0 = x1;
f0 = f1;
} else {
// We must replace x1, this is the case where
// the modification to regula falsi kicks in,
// by modifying f0.
// 1. The classic Illinois method
// const double gamma = 0.5;
// @afr: The next two methods do not work??!!?
// 2. The method called 'Method 3' in the paper.
// const double phi0 = f1/f0;
// const double phi1 = fnew/f1;
// const double gamma = 1.0 - phi1/(1.0 - phi0);
// 3. The method called 'Method 4' in the paper.
// const double phi0 = f1/f0;
// const double phi1 = fnew/f1;
// const double gamma = 1.0 - phi0 - phi1;
// cout << "phi0 = " << phi0 <<" phi1 = " << phi1 <<
// " gamma = " << gamma << endl;
// 4. The 'Pegasus' method
const double gamma = f1/(f1 + fnew);
f0 *= gamma;
}
x1 = xnew;
f1 = fnew;
}
return 0.5*(x0 + x1);
}
/// Implements a modified regula falsi method as described in
/// "Improved algorithms of Illinois-type for the numerical
/// solution of nonlinear equations"
/// by J. A. Ford.
/// Current variant is the 'Pegasus' method.
/// Implements a modified regula falsi method as described in
/// "Improved algorithms of Illinois-type for the numerical
/// solution of nonlinear equations"
/// by J. A. Ford.
/// Current variant is the 'Pegasus' method.
/// This version takes an extra parameter for the initial guess.
template <class Functor>
inline static double solve(const Functor& f,
template <class Functor>
inline static double solve(const Functor& f,
const double initial_guess,
const double a,
const double b,
const int max_iter,
const double tolerance,
int& iterations_used)
{
using namespace std;
const double macheps = numeric_limits<double>::epsilon();
const double eps = tolerance + macheps*max(max(fabs(a), fabs(b)), 1.0);
{
using namespace std;
const double macheps = numeric_limits<double>::epsilon();
const double eps = tolerance + macheps*max(max(fabs(a), fabs(b)), 1.0);
double f_initial = f(initial_guess);
const double epsF = tolerance + macheps*max(fabs(f_initial), 1.0);
const double epsF = tolerance + macheps*max(fabs(f_initial), 1.0);
if (fabs(f_initial) < epsF) {
return initial_guess;
}
double x0 = a;
double x1 = b;
double x0 = a;
double x1 = b;
double f0 = f_initial;
double f1 = f_initial;
if (x0 != initial_guess) {
@ -214,7 +214,7 @@ namespace Opm
if (fabs(f0) < epsF) {
return x0;
}
}
}
if (x1 != initial_guess) {
f1 = f(x1);
if (fabs(f1) < epsF) {
@ -228,65 +228,65 @@ namespace Opm
x0 = initial_guess;
f0 = f_initial;
}
if (f0*f1 > 0.0) {
if (f0*f1 > 0.0) {
return ErrorPolicy::handleBracketingFailure(a, b, f0, f1);
}
iterations_used = 0;
// In every iteraton, x1 is the last point computed,
// and x0 is the last point computed that makes it a bracket.
while (fabs(x1 - x0) >= 1e-9*eps) {
double xnew = regulaFalsiStep(x0, x1, f0, f1);
double fnew = f(xnew);
}
iterations_used = 0;
// In every iteraton, x1 is the last point computed,
// and x0 is the last point computed that makes it a bracket.
while (fabs(x1 - x0) >= 1e-9*eps) {
double xnew = regulaFalsiStep(x0, x1, f0, f1);
double fnew = f(xnew);
// cout << "xnew = " << xnew << " fnew = " << fnew << endl;
++iterations_used;
if (iterations_used > max_iter) {
++iterations_used;
if (iterations_used > max_iter) {
return ErrorPolicy::handleTooManyIterations(x0, x1, max_iter);
}
if (fabs(fnew) < epsF) {
return xnew;
}
// Now we must check which point we must replace.
if ((fnew > 0.0) == (f0 > 0.0)) {
// We must replace x0.
x0 = x1;
f0 = f1;
} else {
// We must replace x1, this is the case where
// the modification to regula falsi kicks in,
// by modifying f0.
// 1. The classic Illinois method
// const double gamma = 0.5;
// @afr: The next two methods do not work??!!?
// 2. The method called 'Method 3' in the paper.
}
if (fabs(fnew) < epsF) {
return xnew;
}
// Now we must check which point we must replace.
if ((fnew > 0.0) == (f0 > 0.0)) {
// We must replace x0.
x0 = x1;
f0 = f1;
} else {
// We must replace x1, this is the case where
// the modification to regula falsi kicks in,
// by modifying f0.
// 1. The classic Illinois method
// const double gamma = 0.5;
// @afr: The next two methods do not work??!!?
// 2. The method called 'Method 3' in the paper.
// const double phi0 = f1/f0;
// const double phi1 = fnew/f1;
// const double gamma = 1.0 - phi1/(1.0 - phi0);
// 3. The method called 'Method 4' in the paper.
// 3. The method called 'Method 4' in the paper.
// const double phi0 = f1/f0;
// const double phi1 = fnew/f1;
// const double gamma = 1.0 - phi0 - phi1;
// cout << "phi0 = " << phi0 <<" phi1 = " << phi1 <<
// " gamma = " << gamma << endl;
// 4. The 'Pegasus' method
const double gamma = f1/(f1 + fnew);
f0 *= gamma;
}
x1 = xnew;
f1 = fnew;
}
return 0.5*(x0 + x1);
}
// 4. The 'Pegasus' method
const double gamma = f1/(f1 + fnew);
f0 *= gamma;
}
x1 = xnew;
f1 = fnew;
}
return 0.5*(x0 + x1);
}
private:
inline static double regulaFalsiStep(const double a,
inline static double regulaFalsiStep(const double a,
const double b,
const double fa,
const double fb)
{
ASSERT(fa*fb < 0.0);
return (b*fa - a*fb)/(fa - fb);
}
{
ASSERT(fa*fb < 0.0);
return (b*fa - a*fb)/(fa - fb);
}
};
@ -294,8 +294,8 @@ namespace Opm
template <class Functor>
inline void bracketZero(const Functor& f,
template <class Functor>
inline void bracketZero(const Functor& f,
const double x0,
const double dx,
double& a,