spline: some fixes to make Spline::intersect() work correctly

most of these issues are fallout from the conversion of the spline
class from a moments (second derivative) based approach to hermite
base functions:

- the second and third derivatives where incorrect, and the third
  derivative was not available in the public interface.
- The unit test did not check the derivatives
- The coefficients for the monomial basis were off by the factors
  stemming from the derivatives
- The intersectIntervall_() method used std::max() instead of
  std::min() at one place and still added the base offset for the x
  values as indicated by Stoer

opm/core/utility/Spline.hpp

@@ -858,6 +858,27 @@ public:
         return evalDerivative2_(x, segmentIdx_(x));
     }
 
+    /*!
+     * \brief Evaluate the spline's third derivative at a given position.
+     *
+     * \param x The value on the abscissa where the spline's
+     *          derivative ought to be evaluated
+     *
+     * \param extrapolate If this parameter is set to true, the spline
+     *                    will be extended beyond its range by
+     *                    straight lines, if false calling extrapolate
+     *                    for \f$ x \not [x_{min}, x_{max}]\f$ will
+     *                    cause a failed assertation.
+     */
+    Scalar evalThirdDerivative(Scalar x, bool extrapolate=false) const
+    {
+        assert(extrapolate || applies(x));
+        if (extrapolate)
+            return 0.0;
+
+        return evalDerivative3_(x, segmentIdx_(x));
+    }
+
     /*!
      * \brief Find the intersections of the spline with a cubic
      *        polynomial in the whole intervall, throws
@@ -1493,7 +1514,7 @@ protected:
         Scalar alpha = 1 / delta;
 
         return
-            alpha
+            alpha*alpha
             *(h00_prime2_(t) * y_(i)
               + h10_prime2_(t) * slope_(i)*delta
               + h01_prime2_(t) * y_(i + 1)
@@ -1505,15 +1526,16 @@ protected:
     Scalar evalDerivative3_(Scalar x, int i) const
     {
         // See http://en.wikipedia.org/wiki/Cubic_Hermite_spline
-        Scalar t = (x - x_(i))/h_(i + 1);
-        Scalar alpha = 1 / h_(i + 1);
+        Scalar delta = h_(i + 1);
+        Scalar t = (x - x_(i))/delta;
+        Scalar alpha = 1 / delta;
 
         return
-            alpha
+            alpha*alpha*alpha
             *(h00_prime3_(t)*y_(i)
-              + h10_prime3_(t)*slope_(i)
+              + h10_prime3_(t)*slope_(i)*delta
               + h01_prime3_(t)*y_(i + 1)
-              + h11_prime3_(t)*slope_(i + 1));
+              + h11_prime3_(t)*slope_(i + 1)*delta);
     }
 
     // hermite basis functions
@@ -1578,8 +1600,8 @@ protected:
     // -1: spline is monotonously decreasing in the specified interval
     int monotonic_(int i, Scalar x0, Scalar x1) const
     {
-        Scalar a = 3*a_(i);
-        Scalar b = 2*b_(i);
+        Scalar a = a_(i);
+        Scalar b = b_(i);
         Scalar c = c_(i);
 
         if (std::abs(a) < 1e-20 && std::abs(b) < 1e-20 && std::abs(c) < 1e-20)
@@ -1631,12 +1653,11 @@ protected:
                                            c_(segIdx) - c,
                                            d_(segIdx) - d);
         x0 = std::max(x_(segIdx), x0);
-        x1 = std::max(x_(segIdx+1), x1);
+        x1 = std::min(x_(segIdx+1), x1);
 
         // filter the intersections outside of the specified intervall
         int k = 0;
         for (int j = 0; j < n; ++j) {
-            sol[j] += x_(segIdx); // add the offset of the intervall. For details see Stoer
             if (x0 <= sol[j] && sol[j] <= x1) {
                 sol[k] = sol[j];
                 ++k;
@@ -1691,19 +1712,20 @@ protected:
     Scalar slope_(int i) const
     { return slopeVec_[i]; }
 
-    // returns the coefficient in front of the x^0 term. In Stoer this
+    // returns the coefficient in front of the x^3 term. In Stoer this
     // is delta.
     Scalar a_(int i) const
-    { return evalDerivative3_(/*x=*/0, i); }
+    { return evalDerivative3_(/*x=*/0, i)/6.0; }
 
     // returns the coefficient in front of the x^2 term In Stoer this
     // is gamma.
     Scalar b_(int i) const
-    { return evalDerivative2_(/*x=*/0, i); }
+    { return evalDerivative2_(/*x=*/0, i)/2.0; }
 
     // returns the coefficient in front of the x^1 term. In Stoer this
     // is beta.
-    Scalar c_(int i) const { return evalDerivative_(/*x=*/0, i); }
+    Scalar c_(int i) const
+    { return evalDerivative_(/*x=*/0, i); }
 
     // returns the coefficient in front of the x^0 term. In Stoer this
     // is alpha.

tests/test_spline.cpp

@@ -50,6 +50,7 @@ void testCommon(const Spline &sp,
                 const double *y)
 {
     static double eps = 1e-10;
+    static double epsFD = 1e-7;
 
     int n = sp.numSamples();
     for (int i = 0; i < n; ++i) {
@@ -73,7 +74,49 @@ void testCommon(const Spline &sp,
 
         if (std::abs(d1 - d0) > 1000*eps || std::abs(d2 - d0) > 1000*eps)
             OPM_THROW(std::runtime_error,
-                       "Spline seems to exhibit a discontinuous derivative at sampling point " << i << "!");
+                      "Spline seems to exhibit a discontinuous derivative at sampling point " << i << "!");
+    }
+
+    // make sure the derivatives are consistent with the curve
+    int np = 3*n;
+    for (int i = 0; i < np; ++i) {
+        double x = sp.xMin() + (sp.xMax() - sp.xMin())*i/np;
+
+        // first derivative
+        double y1 = sp.eval(x+epsFD);
+        double y0 = sp.eval(x);
+
+        double mFD = (y1 - y0)/epsFD;
+        double m = sp.evalDerivative(x);
+
+        if (std::abs( mFD - m ) > 1000*epsFD)
+            OPM_THROW(std::runtime_error,
+                      "Derivative of spline seems to be inconsistent with cuve"
+                      " (" << mFD << " - " << m << " = " << mFD - m << ")!");
+
+        // second derivative
+        y1 = sp.evalDerivative(x+epsFD);
+        y0 = sp.evalDerivative(x);
+
+        mFD = (y1 - y0)/epsFD;
+        m = sp.evalSecondDerivative(x);
+
+        if (std::abs( mFD - m ) > 1000*epsFD)
+            OPM_THROW(std::runtime_error,
+                      "Second derivative of spline seems to be inconsistent with cuve"
+                      " (" << mFD << " - " << m << " = " << mFD - m << ")!");
+
+        // Third derivative
+        y1 = sp.evalSecondDerivative(x+epsFD);
+        y0 = sp.evalSecondDerivative(x);
+
+        mFD = (y1 - y0)/epsFD;
+        m = sp.evalThirdDerivative(x);
+
+        if (std::abs( mFD - m ) > 1000*epsFD)
+            OPM_THROW(std::runtime_error,
+                      "Third derivative of spline seems to be inconsistent with cuve"
+                      " (" << mFD << " - " << m << " = " << mFD - m << ")!");
     }
 }