Merge pull request #4674 from OPM/feature-ensemble-rft

#4634 and #4635: Merge feature-ensemble-rft branch

ApplicationCode/ProjectDataModel/RimWellLogRftCurve.cpp

@@ -409,7 +409,7 @@ void RimWellLogRftCurve::onLoadDataAndUpdate(bool updateParentPlot)
 
         if (tvDepthVector.size() != measuredDepthVector.size())
         {
-            measuredDepthVector = interpolatedMeasuredDepthValuesFromObservedData(tvDepthVector);
+            measuredDepthVector = interpolatedMeasuredDepthValuesFromWellPathOrObservedData(tvDepthVector);
             if (measuredDepthVector.size() == tvDepthVector.size())
             {
                 m_derivingMDFromObservedData = true;
@@ -900,23 +900,34 @@ std::vector<double> RimWellLogRftCurve::measuredDepthValues()
 //--------------------------------------------------------------------------------------------------
 ///
 //--------------------------------------------------------------------------------------------------
-std::vector<double> RimWellLogRftCurve::interpolatedMeasuredDepthValuesFromObservedData(const std::vector<double>& tvDepthValues)
+std::vector<double> RimWellLogRftCurve::interpolatedMeasuredDepthValuesFromWellPathOrObservedData(const std::vector<double>& tvDepthValues)
 {
+    RimProject*  proj = RiaApplication::instance()->project();
+    RimWellPath* wellPath = proj->wellPathByName(m_wellName);
+
     std::vector<double> interpolatedMdValues;
-    if (m_observedFmuRftData)
+    if (wellPath)
+    {
+        const std::vector<double>& mdValuesOfWellPath  = wellPath->wellPathGeometry()->measureDepths();
+        std::vector<double>        tvdValuesOfWellPath = wellPath->wellPathGeometry()->trueVerticalDepths();
+
+        interpolatedMdValues = RigWellPathGeometryTools::interpolateMdFromTvd(mdValuesOfWellPath, tvdValuesOfWellPath, tvDepthValues);
+        CVF_ASSERT(interpolatedMdValues.size() == tvDepthValues.size());
+    }
+    else if (m_observedFmuRftData)
     {
         RifReaderRftInterface* reader = m_observedFmuRftData->rftReader();
         if (reader)
         {
-            std::vector<double>    tvdValuesOfObbservedData;
+            std::vector<double>    tvdValuesOfObservedData;
-            std::vector<double>    mdValuesOfObbservedData;
+            std::vector<double>    mdValuesOfObservedData;
 
             RifEclipseRftAddress tvdAddress(m_wellName(), m_timeStep, RifEclipseRftAddress::TVD);
             RifEclipseRftAddress mdAddress(m_wellName(), m_timeStep, RifEclipseRftAddress::MD);
 
-            reader->values(tvdAddress, &tvdValuesOfObbservedData);
+            reader->values(tvdAddress, &tvdValuesOfObservedData);
-            reader->values(mdAddress, &mdValuesOfObbservedData);
+            reader->values(mdAddress, &mdValuesOfObservedData);
-            interpolatedMdValues = RigWellPathGeometryTools::interpolateMdFromTvd(mdValuesOfObbservedData, tvdValuesOfObbservedData, tvDepthValues);
+            interpolatedMdValues = RigWellPathGeometryTools::interpolateMdFromTvd(mdValuesOfObservedData, tvdValuesOfObservedData, tvDepthValues);
             CVF_ASSERT(interpolatedMdValues.size() == tvDepthValues.size());
         }
     }

ApplicationCode/ProjectDataModel/RimWellLogRftCurve.h

@@ -96,7 +96,7 @@ private:
     std::vector<double>  xValues();
     std::vector<double>  tvDepthValues();
     std::vector<double>  measuredDepthValues();
-    std::vector<double>  interpolatedMeasuredDepthValuesFromObservedData(const std::vector<double>& tvDepthValues);
+    std::vector<double> interpolatedMeasuredDepthValuesFromWellPathOrObservedData(const std::vector<double>& tvDepthValues);
 
 private:
     caf::PdmPtrField<RimEclipseResultCase*>     m_eclipseResultCase;

ApplicationCode/ReservoirDataModel/RigWellPath.cpp

@@ -408,3 +408,16 @@ const std::vector<double>& RigWellPath::measureDepths() const
     return m_measuredDepths;
 }
 
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+std::vector<double> RigWellPath::trueVerticalDepths() const
+{
+    std::vector<double> tvds;
+    for (const cvf::Vec3d& point : m_wellPathPoints)
+    {
+        tvds.push_back(std::fabs(point.z()));
+    }
+    return tvds;
+}
+

ApplicationCode/ReservoirDataModel/RigWellPath.h

@@ -42,6 +42,7 @@ public:
 
     const std::vector<cvf::Vec3d>&  wellPathPoints() const;
     const std::vector<double>&      measureDepths() const;
+    std::vector<double>             trueVerticalDepths() const;
 
     RigWellPath();
     void                        setDatumElevation(double value);

ApplicationCode/ReservoirDataModel/RigWellPathGeometryTools.cpp

@@ -94,52 +94,89 @@ std::vector<double> RigWellPathGeometryTools::interpolateMdFromTvd(const std::ve
     std::vector<double> interpolatedMdValues;
     interpolatedMdValues.reserve(tvdValuesToInterpolateFrom.size());
 
-    std::vector<double>::const_iterator last_it = originalTvdValues.begin();
+    QwtSpline        spline              = createSpline(originalMdValues, originalTvdValues);
-    for (std::vector<double>::const_iterator it = tvdValuesToInterpolateFrom.begin(); it != tvdValuesToInterpolateFrom.end(); ++it)
+    std::vector<int> segmentStartIndices = findSplineSegmentsContainingRoots(spline, tvdValuesToInterpolateFrom);
+
+    for (size_t i = 0; i < segmentStartIndices.size(); ++i)
     {
-        double tvdValue = *it;
+        double currentTVDValue = tvdValuesToInterpolateFrom[i];
-        double sign = 0.0;
+        double startMD = spline.points().front().x();
-        if (it != tvdValuesToInterpolateFrom.begin())
+        double endMD = spline.points().back().y();
+        if (segmentStartIndices[i] != -1)
         {
-            sign = *it - *(it - 1);
+            int startIndex = segmentStartIndices[i];
+            int endIndex = startIndex + 1;
+
+            // Search interval for best MD value
+            startMD = spline.points()[startIndex].x();
+            endMD = spline.points().back().y();
+
+            if (endIndex < spline.points().size())
+            {
+                if (!interpolatedMdValues.empty())
+                {
+                    double mdDiff = 0.0;
+                    if (interpolatedMdValues.size() > 1)
+                    {
+                        mdDiff = interpolatedMdValues[i - 1] - interpolatedMdValues[i - 2];
+                    }
+                    startMD = std::max(startMD, interpolatedMdValues.back() + 0.1 * mdDiff);
+                }
+                endMD = spline.points()[endIndex].x();
+            }
+        }
+        double mdValue = solveForX(spline, startMD, endMD, currentTVDValue);
+        interpolatedMdValues.push_back(mdValue);
+    }
+    return interpolatedMdValues;
+}
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+std::vector<int> RigWellPathGeometryTools::findSplineSegmentsContainingRoots(const QwtSpline& spline, const std::vector<double>& tvdValuesToInterpolateFrom)
+{
+    std::vector<int> segmentStartIndices;
+    segmentStartIndices.reserve(tvdValuesToInterpolateFrom.size());
+
+    int lastSplineStartIndex = 0;
+    for (double tvdValue : tvdValuesToInterpolateFrom)
+    {
+        int currentSplineStartIndex = lastSplineStartIndex;
+       
+        bool foundMatch = false;
+        // Increment current_it until we find an interval containing our TVD
+        while (currentSplineStartIndex < spline.points().size() - 2)
+        {
+            double diffCurrent = spline.points()[currentSplineStartIndex].y() - tvdValue;
+            if (std::abs(diffCurrent) < 1.0e-8) // Current is matching the point
+            {
+                foundMatch = true;
+                break;
+            }
+
+            int nextStartIndex = currentSplineStartIndex + 1;
+
+            double diffNext = spline.points()[nextStartIndex].y() - tvdValue;
+            if (diffCurrent * diffNext < 0.0) // One is above, the other is below
+            {
+                foundMatch = true;
+                break;
+            }
+            currentSplineStartIndex = nextStartIndex;
+        }
+        if (foundMatch)
+        {
+            segmentStartIndices.push_back(currentSplineStartIndex);
+            lastSplineStartIndex = currentSplineStartIndex;
         }
         else
         {
-            sign = *(it + 1) - *it;
+            segmentStartIndices.push_back(-1);
         }
-        if (std::fabs(sign) < 1.0e-8)
-        {
-            continue;
     }
 
-        sign /= std::fabs(sign);
+    return segmentStartIndices;
-
-        auto current_it = last_it;
-        // Is incrementing current_it taking us closer to the TVD value we want?
-        while (current_it != originalTvdValues.end())
-        {
-            if (*current_it * sign >= tvdValue * sign)
-            {
-                break;
-            }
-
-            auto next_it = current_it + 1;
-            if (next_it != originalTvdValues.end())
-            {                
-                double originalDataSign = (*next_it - *current_it);
-                originalDataSign /= std::fabs(originalDataSign);
-                if (originalDataSign * sign < 0.0)
-                    break;
-            }
-            current_it = next_it;
-        }
-
-        int valueIndex = static_cast<int>(current_it - originalTvdValues.begin());
-        double mdValue = linearInterpolation(originalTvdValues, originalMdValues, valueIndex, tvdValue);
-        interpolatedMdValues.push_back(mdValue);
-        last_it = current_it;
-    }
-    return interpolatedMdValues;
 }
 
 std::vector<cvf::Vec3d> RigWellPathGeometryTools::interpolateUndefinedNormals(const cvf::Vec3d&              planeNormal,
@@ -224,20 +261,93 @@ cvf::Vec3d RigWellPathGeometryTools::estimateDominantDirectionInXYPlane(const st
     return directionSum.getNormalized();
 }
 
+//--------------------------------------------------------------------------------------------------
+/// Golden-section minimization: https://en.wikipedia.org/wiki/Golden-section_search
+//--------------------------------------------------------------------------------------------------
+double RigWellPathGeometryTools::solveForX(const QwtSpline& spline, double minX, double maxX, double y)
+{
+    const double phi = (1.0 + std::sqrt(5.0)) / 2.0;
+    const double tol = 1.0e-8;
+
+    double a = minX, b = maxX;
+    double c = b - (b - a) / phi;
+    double d = a + (b - a) / phi;
+
+    double fa = spline.value(a) - y;
+    double fb = spline.value(b) - y;
+    double fc = spline.value(c) - y;
+    double fd = spline.value(d) - y;
+
+    for (int n = 0; n < 100; ++n)
+    {
+        if (std::fabs(c - d) < tol)
+        {
+            break;
+        }
+        
+        if (std::fabs(fc) < std::fabs(fd))
+        {
+            b = d;
+            fb = fd;
+            d = c;
+            fd = fc;
+            c = b - (b - a) / phi;
+            fc = spline.value(c) - y;
+        }
+        else
+        {
+            a = c;
+            fa = fc;
+            c = d;
+            fc = fd;
+            d = a + (b - a) / phi;
+            fd = spline.value(d) - y;
+        }
+    }
+    return (a + b) / 2.0;
+}
+
 //--------------------------------------------------------------------------------------------------
 ///
 //--------------------------------------------------------------------------------------------------
-double RigWellPathGeometryTools::linearInterpolation(const std::vector<double>& xValues, const std::vector<double>& yValues, int valueIndex, double x)
+QwtSpline RigWellPathGeometryTools::createSpline(const std::vector<double>& originalMdValues, const std::vector<double>& originalTvdValues)
 {
-    int N = (int) xValues.size() - 1;
+    QPolygonF polygon;
-    int i = cvf::Math::clamp(valueIndex, 1, N);
+    for (size_t i = 0; i < originalMdValues.size(); ++i)
-        
+    {
-    double x1 = xValues[i - 1];
+        polygon << QPointF(originalMdValues[i], originalTvdValues[i]);
-    double x2 = xValues[i];
+    }
-    double y1 = yValues[i - 1];
+    QwtSplineCurveFitter curveFitter;
-    double y2 = yValues[i];
+    QPolygonF splinePoints = curveFitter.fitCurve(polygon);
 
+    // Extend spline from 0.0 to a large value for MD
+    // This is to force a specific and known extrapolation.
+    // Otherwise we get an undefined and unknown extrapolation.
+    {
+        double x1 = splinePoints[0].x();
+        double x2 = splinePoints[1].x();
+        double y1 = splinePoints[0].y();
+        double y2 = splinePoints[1].y();
         double M = (y2 - y1) / (x2 - x1);
 
-    return M * (x - x1) + y1;
+        QPointF startPoint(0.0f, M * (0.0f - x1) + y1);
+        splinePoints.push_front(startPoint);
+    }
+    {
+        int N = splinePoints.size() - 1;
+        double x1 = splinePoints[N - 1].x();
+        double x2 = splinePoints[N].x();
+        double y1 = splinePoints[N - 1].y();
+        double y2 = splinePoints[N].y();
+        double M = (y2 - y1) / (x2 - x1);
+        double endX = 2.0 * splinePoints[N].x();
+
+        QPointF endPoint(endX, M  * (endX - x1) + y1);
+        splinePoints.push_back(endPoint);
+    }
+
+    QwtSpline spline;
+    spline.setPoints(splinePoints);
+
+    return spline;
 }

ApplicationCode/ReservoirDataModel/RigWellPathGeometryTools.h

@@ -18,11 +18,14 @@
 
 #pragma once
 
-#include <vector>
-
 #include "cvfBase.h"
 #include "cvfVector3.h"
 
+#include <qwt_curve_fitter.h>
+#include <qwt_spline.h>
+
+#include <vector>
+
 class RigWellPath;
 
 //==================================================================================================
@@ -38,9 +41,10 @@ public:
     };
 
 public:
-    static std::vector<cvf::Vec3d> calculateLineSegmentNormals(const std::vector<cvf::Vec3d>& vertices,
+    static std::vector<cvf::Vec3d> calculateLineSegmentNormals(const std::vector<cvf::Vec3d>& vertices, double angle);
-                                                               double                         angle);
+    static std::vector<double>     interpolateMdFromTvd(const std::vector<double>& originalMdValues,
-    static std::vector<double> interpolateMdFromTvd(const std::vector<double>& originalMdValues, const std::vector<double>& originalTvdValues, const std::vector<double>& tvdValuesToInterpolateFrom);
+                                                        const std::vector<double>& originalTvdValues,
+                                                        const std::vector<double>& tvdValuesToInterpolateFrom);
 
 private:
     static std::vector<cvf::Vec3d> interpolateUndefinedNormals(const cvf::Vec3d&              planeNormal,
@@ -48,6 +52,9 @@ private:
                                                                const std::vector<cvf::Vec3d>& vertices);
     static cvf::Vec3d              estimateDominantDirectionInXYPlane(const std::vector<cvf::Vec3d>& vertices);
 
-    static double linearInterpolation(const std::vector<double>& xValues, const std::vector<double>& yValues, int valueIndex, double x);
+    static double solveForX(const QwtSpline& spline, double minX, double maxX, double y);
 
+    static QwtSpline createSpline(const std::vector<double>& originalMdValues, const std::vector<double>& originalTvdValues);
+    static std::vector<int> findSplineSegmentsContainingRoots(const QwtSpline&           spline,
+                                                              const std::vector<double>& tvdValuesToInterpolateFrom);
 };

ApplicationCode/UnitTests/RigWellPathGeometryTools-Test.cpp

@@ -20,9 +20,13 @@
 
 #include "RigWellPathGeometryTools.h"
 
+#include <algorithm>
 #include <vector>
 #include <complex>
 
+#define TOLERANCE 1.0e-7
+#define SOLVER_TOLERANCE 
+
 //--------------------------------------------------------------------------------------------------
 ///
 //--------------------------------------------------------------------------------------------------
@@ -36,7 +40,7 @@ TEST(RigWellPathGeometryTools, VerticalPath)
     EXPECT_EQ(fullTVDValues.size(), fullMDValues.size());
     for (size_t i = 0; i < fullTVDValues.size(); ++i)
     {
-        EXPECT_DOUBLE_EQ(fullTVDValues[i], fullMDValues[i]);
+        EXPECT_NEAR(fullTVDValues[i], fullMDValues[i], TOLERANCE);
     }
 }
 
@@ -45,12 +49,13 @@ TEST(RigWellPathGeometryTools, LinearPath)
     std::vector<double> mdValues = {100, 500, 1000};
     std::vector<double> tvdValues = {50, 250, 500};
     std::vector<double> fullTVDValues = {50, 100, 150, 200, 250, 300, 350, 400, 450, 500};
+
     std::vector<double> fullMDValues = RigWellPathGeometryTools::interpolateMdFromTvd(mdValues, tvdValues, fullTVDValues);
 
     EXPECT_EQ(fullTVDValues.size(), fullMDValues.size());
     for (size_t i = 0; i < fullTVDValues.size(); ++i)
     {
-        EXPECT_DOUBLE_EQ(2.0*fullTVDValues[i], fullMDValues[i]);
+        EXPECT_NEAR(2.0*fullTVDValues[i], fullMDValues[i], TOLERANCE);
     }
 }
 
@@ -73,11 +78,15 @@ TEST(RigWellPathGeometryTools, QuadraticPath)
     {
         fullTvdValues.push_back(quadraticFunction(md));
     }
+    
     std::vector<double> estimatedFullMDValues = RigWellPathGeometryTools::interpolateMdFromTvd(mdValues, tvdValues, fullTvdValues);
     EXPECT_EQ(estimatedFullMDValues.size(), fullMDValues.size());
     for (size_t i = 0; i < estimatedFullMDValues.size(); ++i)
     {
-        EXPECT_DOUBLE_EQ(fullMDValues[i], estimatedFullMDValues[i]);
+        if (std::find(mdValues.begin(), mdValues.end(), estimatedFullMDValues[i]) != mdValues.end())
+        {
+            EXPECT_NEAR(fullMDValues[i], estimatedFullMDValues[i], TOLERANCE);
+        }
     }
 }
 
@@ -87,6 +96,32 @@ double cubicFunction(double x)
 }
 
 TEST(RigWellPathGeometryTools, CubicPath)
+{
+    std::vector<double> mdValues = {100, 300, 700, 1000};
+    std::vector<double> tvdValues;
+    for (double md : mdValues)
+    {
+        tvdValues.push_back(cubicFunction(md));
+    }
+    std::vector<double> fullMDValues = {100, 200, 300, 400, 500, 600, 700, 800, 900, 1000};
+    std::vector<double> fullTvdValues;
+    for (double md : fullMDValues)
+    {
+        fullTvdValues.push_back(cubicFunction(md));
+    }
+
+    std::vector<double> estimatedFullMDValues = RigWellPathGeometryTools::interpolateMdFromTvd(mdValues, tvdValues, fullTvdValues);
+    EXPECT_EQ(estimatedFullMDValues.size(), fullMDValues.size());
+    for (size_t i = 0; i < estimatedFullMDValues.size(); ++i)
+    {
+        if (std::find(mdValues.begin(), mdValues.end(), estimatedFullMDValues[i]) != mdValues.end())
+        {
+            EXPECT_NEAR(fullMDValues[i], estimatedFullMDValues[i], TOLERANCE);
+        }
+    }
+}
+
+TEST(RigWellPathGeometryTools, CubicPathPoorSampling)
 {
     std::vector<double> mdValues = {100, 300, 600, 1000};
     std::vector<double> tvdValues;
@@ -100,10 +135,14 @@ TEST(RigWellPathGeometryTools, CubicPath)
     {
         fullTvdValues.push_back(cubicFunction(md));
     }
+
     std::vector<double> estimatedFullMDValues = RigWellPathGeometryTools::interpolateMdFromTvd(mdValues, tvdValues, fullTvdValues);
     EXPECT_EQ(estimatedFullMDValues.size(), fullMDValues.size());
     for (size_t i = 0; i < estimatedFullMDValues.size(); ++i)
     {
-        EXPECT_DOUBLE_EQ(fullMDValues[i], estimatedFullMDValues[i]);
+        if (std::find(mdValues.begin(), mdValues.end(), estimatedFullMDValues[i]) != mdValues.end())
+        {
+            EXPECT_NEAR(fullMDValues[i], estimatedFullMDValues[i], TOLERANCE);
+        }
     }
 }