From 288679ff3a735ca16a1676ea7eb3c317dcf9621e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 24 May 2012 10:02:14 +0200 Subject: [PATCH] Policy-based design enables custom error handling in scalar regula falsi. --- .../reorder/TransportModelTwophase.cpp | 4 +- opm/core/utility/RootFinders.hpp | 84 +++++++++++++------ 2 files changed, 59 insertions(+), 29 deletions(-) diff --git a/opm/core/transport/reorder/TransportModelTwophase.cpp b/opm/core/transport/reorder/TransportModelTwophase.cpp index 1c7dbb53..c920deb5 100644 --- a/opm/core/transport/reorder/TransportModelTwophase.cpp +++ b/opm/core/transport/reorder/TransportModelTwophase.cpp @@ -168,7 +168,7 @@ namespace Opm // } int iters_used; // saturation_[cell] = modifiedRegulaFalsi(res, smin_[2*cell], smax_[2*cell], maxit_, tol_, iters_used); - saturation_[cell] = modifiedRegulaFalsi(res, saturation_[cell], 0.0, 1.0, maxit_, tol_, iters_used); + saturation_[cell] = RegulaFalsi<>::solve(res, saturation_[cell], 0.0, 1.0, maxit_, tol_, iters_used); fractionalflow_[cell] = fracFlow(saturation_[cell], cell); } @@ -541,7 +541,7 @@ namespace Opm GravityResidual res(*this, cells, pos, gravflux); if (std::fabs(res(saturation_[cell])) > tol_) { int iters_used; - saturation_[cell] = modifiedRegulaFalsi(res, smin_[2*cell], smax_[2*cell], maxit_, tol_, iters_used); + saturation_[cell] = RegulaFalsi<>::solve(res, smin_[2*cell], smax_[2*cell], maxit_, tol_, iters_used); } saturation_[cell] = std::min(std::max(saturation_[cell], smin_[2*cell]), smax_[2*cell]); mobility(saturation_[cell], cell, &mob_[2*cell]); diff --git a/opm/core/utility/RootFinders.hpp b/opm/core/utility/RootFinders.hpp index cbd0f294..6c2c3a14 100644 --- a/opm/core/utility/RootFinders.hpp +++ b/opm/core/utility/RootFinders.hpp @@ -46,14 +46,33 @@ namespace Opm { - inline 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); - } + struct ThrowOnError + { + static double handleBracketingFailure(const double x0, const double x1, const double f0, const double f1) + { + THROW("Error in parameters, zero not bracketed: [a, b] = [" + << x0 << ", " << x1 << "] f(a) = " << f0 << " f(b) = " << f1); + return -1e100; // Never reached. + } + static double handleTooManyIterations(const double x0, const double x1, const int maxiter) + { + THROW("Maximum number of iterations exceeded: " << maxiter << "\n" + << "Current interval is [" << std::min(x0, x1) << ", " + << std::max(x0, x1) << "]"); + } + }; + + + struct ContinueOnError + { + }; + + + + template + class RegulaFalsi + { + public: /// Implements a modified regula falsi method as described in @@ -61,13 +80,13 @@ namespace Opm /// solution of nonlinear equations" /// by J. A. Ford. /// Current variant is the 'Pegasus' method. - template - inline double modifiedRegulaFalsi(const Functor& f, - const double a, - const double b, - const int max_iter, - const double tolerance, - int& iterations_used) + template + 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::epsilon(); @@ -85,8 +104,7 @@ namespace Opm return x1; } if (f0*f1 > 0.0) { - THROW("Error in parameters, zero not bracketed: [a, b] = [" - << a << ", " << b << "] fa = " << f0 << " fb = " << f1); + return ErrorPolicy::handleBracketingFailure(a, b, f0, f1); } iterations_used = 0; // In every iteraton, x1 is the last point computed, @@ -97,9 +115,7 @@ namespace Opm // cout << "xnew = " << xnew << " fnew = " << fnew << endl; ++iterations_used; if (iterations_used > max_iter) { - THROW("Maximum number of iterations exceeded.\n" - << "Current interval is [" << min(x0, x1) << ", " - << max(x0, x1) << "]"); + return ErrorPolicy::handleTooManyIterations(x0, x1, max_iter); } if (fabs(fnew) < epsF) { return xnew; @@ -144,13 +160,13 @@ namespace Opm /// Current variant is the 'Pegasus' method. /// This version takes an extra parameter for the initial guess. template - inline double modifiedRegulaFalsi(const Functor& f, - const double initial_guess, - const double a, - const double b, - const int max_iter, - const double tolerance, - int& iterations_used) + 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::epsilon(); @@ -237,6 +253,20 @@ namespace Opm } + private: + 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); + } + + + }; + + template