Policy-based design enables custom error handling in scalar regula falsi.

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]);

opm/core/utility/RootFinders.hpp

@@ -46,14 +46,33 @@
 namespace Opm
 {
 
-	inline double regulaFalsiStep(const double a,
+    struct ThrowOnError
-                                      const double b,
+    {
-                                      const double fa,
+        static double handleBracketingFailure(const double x0, const double x1, const double f0, const double f1)
-                                      const double fb)
+        {
-	{
+            THROW("Error in parameters, zero not bracketed: [a, b] = ["
-	    ASSERT(fa*fb < 0.0);
+                  << x0 << ", " << x1 << "]    f(a) = " << f0 << "   f(b) = " << f1);
-	    return (b*fa - a*fb)/(fa - fb);
+            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 ErrorPolicy = ThrowOnError>
+    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 <class Functor>
+        template <class Functor>
-	inline double modifiedRegulaFalsi(const Functor& f,
+	inline static double solve(const Functor& f,
-                                          const double a,
+                            const double a,
-                                          const double b,
+                            const double b,
-                                          const int max_iter,
+                            const int max_iter,
-                                          const double tolerance,
+                            const double tolerance,
-                                          int& iterations_used)
+                            int& iterations_used)
 	{
 	    using namespace std;
 	    const double macheps = numeric_limits<double>::epsilon();
@@ -85,8 +104,7 @@ namespace Opm
 		return x1;
 	    }
 	    if (f0*f1 > 0.0) {
-		THROW("Error in parameters, zero not bracketed: [a, b] = ["
+                return ErrorPolicy::handleBracketingFailure(a, b, f0, f1);
-		      << a << ", " << b << "]    fa = " << f0 << "   fb = " << 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"
+                    return ErrorPolicy::handleTooManyIterations(x0, x1, max_iter);
-			  << "Current interval is [" << min(x0, x1) << ", "
-			  << max(x0, x1) << "]");
 		}
 		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 <class Functor>
-	inline double modifiedRegulaFalsi(const Functor& f,
+	inline static double solve(const Functor& f,
-                                          const double initial_guess,
+                            const double initial_guess,
-                                          const double a,
+                            const double a,
-                                          const double b,
+                            const double b,
-                                          const int max_iter,
+                            const int max_iter,
-                                          const double tolerance,
+                            const double tolerance,
-                                          int& iterations_used)
+                            int& iterations_used)
 	{
 	    using namespace std;
 	    const double macheps = numeric_limits<double>::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 <class Functor>