Interpolating MDvalues

ApplicationCode/ProjectDataModel/Flow/RimWellRftPlot.cpp

@@ -450,6 +450,7 @@ void RimWellRftPlot::updateCurvesInPlot(const std::set<RiaRftPltCurveDefinition>
             auto rftCase = curveDefToAdd.address().summaryCase();
             curve->setSummaryCase(rftCase);
             curve->setEnsemble(curveDefToAdd.address().ensemble());
+            curve->setObservedFmuRftData(this->findObservedFmuData(m_wellPathNameOrSimWellName, curveDefToAdd.timeStep()));
             RifEclipseRftAddress address(m_wellPathNameOrSimWellName, curveDefToAdd.timeStep(), RifEclipseRftAddress::PRESSURE);
             curve->setRftAddress(address);
             curve->setZOrder(1);
@@ -472,6 +473,7 @@ void RimWellRftPlot::updateCurvesInPlot(const std::set<RiaRftPltCurveDefinition>
                     plotTrack->addCurve(curve);
                     curve->setEnsemble(ensemble);
                     curve->setRftAddress(rftAddress);
+                    curve->setObservedFmuRftData(this->findObservedFmuData(m_wellPathNameOrSimWellName, curveDefToAdd.timeStep()));
                     curve->setZOrder(RiuQwtPlotCurve::Z_ENSEMBLE_STAT_CURVE);
                     applyCurveAppearance(curve);
                     auto                        symbol   = statisticsCurveSymbolFromAddress(rftAddress);
@@ -1124,3 +1126,19 @@ void RimWellRftPlot::defineCurveColorsAndSymbols(const std::set<RiaRftPltCurveDe
         }
     }
 }
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+RimObservedFmuRftData* RimWellRftPlot::findObservedFmuData(const QString& wellPathName, const QDateTime& timeStep) const
+{
+    auto allObservedDataForWell = RimWellPlotTools::observedFmuRftDataForWell(wellPathName);
+    for (auto observedData : allObservedDataForWell)
+    {
+        if (observedData->rftReader()->availableTimeSteps(wellPathName).count(timeStep))
+        {
+            return observedData;
+        }
+    }
+    return nullptr;
+}

ApplicationCode/ProjectDataModel/Flow/RimWellRftPlot.h

@@ -116,6 +116,8 @@ private:
     void syncCurvesFromUiSelection();
     void assignWellPathToExtractionCurves();
 
+    RimObservedFmuRftData* findObservedFmuData(const QString& wellPathName, const QDateTime& timeStep) const;
+
     std::set<RiaRftPltCurveDefinition> selectedCurveDefs() const;
     std::set<RiaRftPltCurveDefinition> curveDefsFromCurves() const;
 
@@ -161,5 +163,4 @@ private:
 	std::map<RifDataSourceForRftPlt, cvf::Color3f>     m_dataSourceColors;
     std::map<QDateTime, RiuQwtSymbol::PointSymbolEnum> m_timeStepSymbols;
     bool m_isOnLoad;
-
 };

ApplicationCode/ProjectDataModel/RimObservedFmuRftData.cpp

@@ -72,12 +72,12 @@ RifReaderRftInterface* RimObservedFmuRftData::rftReader()
 //--------------------------------------------------------------------------------------------------
 ///
 //--------------------------------------------------------------------------------------------------
-bool RimObservedFmuRftData::hasWell(const QString& simWellName) const
+bool RimObservedFmuRftData::hasWell(const QString& wellPathName) const
 {
     std::vector<QString> allWells = wells();
 	for (const QString& well : allWells)
 	{
-		if (well == simWellName)
+		if (well == wellPathName)
 		{
             return true;
 		}

ApplicationCode/ProjectDataModel/RimObservedFmuRftData.h

@@ -36,9 +36,11 @@ public:
 	void createRftReaderInterface();
     RifReaderRftInterface* rftReader();
 
-	bool hasWell(const QString& simWellName) const;
+	bool hasWell(const QString& wellPathName) const;
 	std::vector<QString> wells() const;
 
+    std::vector<std::pair<double, double>> derivedWellPathTvdMd(const QString& wellPathName) const;
+
 private:
     cvf::ref<RifReaderFmuRft> m_fmuRftReader;
 

ApplicationCode/ProjectDataModel/RimWellLogRftCurve.cpp

@@ -31,6 +31,7 @@
 #include "RigMainGrid.h"
 #include "RigWellLogCurveData.h"
 #include "RigWellPath.h"
+#include "RigWellPathGeometryTools.h"
 #include "RigWellPathIntersectionTools.h"
 
 #include "RimEclipseResultCase.h"
@@ -119,6 +120,7 @@ RimWellLogRftCurve::RimWellLogRftCurve()
     CAF_PDM_InitFieldNoDefault(&m_wellLogChannelName, "WellLogChannelName", "Well Property", "", "", "");
 
     m_isUsingPseudoLength = false;
+    m_derivingMDFromObservedData = false;
 }
 
 //--------------------------------------------------------------------------------------------------
@@ -405,6 +407,15 @@ void RimWellLogRftCurve::onLoadDataAndUpdate(bool updateParentPlot)
 			CVF_ASSERT(false && "Need to have either an eclipse result case, a summary case or an ensemble");
 		}
 
+        if (tvDepthVector.size() != measuredDepthVector.size())
+        {
+            measuredDepthVector = interpolatedMeasuredDepthValuesFromObservedData(tvDepthVector);
+            if (measuredDepthVector.size() == tvDepthVector.size())
+            {
+                m_derivingMDFromObservedData = true;
+            }
+        }
+
         if (tvDepthVector.size() != measuredDepthVector.size())
         {
             m_isUsingPseudoLength = true;
@@ -432,9 +443,10 @@ void RimWellLogRftCurve::onLoadDataAndUpdate(bool updateParentPlot)
                                        m_curveData->trueDepthPlotValues(displayUnit).data(),
                                        static_cast<int>(m_curveData->xPlotValues().size()));
             m_isUsingPseudoLength = false;
+            m_derivingMDFromObservedData = false;
         }
 
-        if (m_isUsingPseudoLength)
+        if (m_isUsingPseudoLength || m_derivingMDFromObservedData)
         {
             RimWellLogTrack* wellLogTrack;
             firstAncestorOrThisOfType(wellLogTrack);
@@ -443,6 +455,9 @@ void RimWellLogRftCurve::onLoadDataAndUpdate(bool updateParentPlot)
             RiuWellLogTrack* viewer = wellLogTrack->viewer();
             if (viewer)
             {
+                if (m_derivingMDFromObservedData)
+                    viewer->setDepthTitle("OBS/" + wellLogPlot->depthPlotTitle());
+                else
                     viewer->setDepthTitle("PL/" + wellLogPlot->depthPlotTitle());
             }
         }
@@ -835,7 +850,7 @@ std::vector<double> RimWellLogRftCurve::tvDepthValues()
 //--------------------------------------------------------------------------------------------------
 std::vector<double> RimWellLogRftCurve::measuredDepthValues()
 {
-	if (m_observedFmuRftData)
+	if (m_observedFmuRftData && !m_ensemble && !m_summaryCase)
 	{
         RifReaderRftInterface* reader = rftReader();
         std::vector<double>    values;
@@ -881,3 +896,29 @@ std::vector<double> RimWellLogRftCurve::measuredDepthValues()
 
     return measuredDepthForCellsWhichHasRftData;
 }
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+std::vector<double> RimWellLogRftCurve::interpolatedMeasuredDepthValuesFromObservedData(const std::vector<double>& tvDepthValues)
+{
+    std::vector<double> interpolatedMdValues;
+    if (m_observedFmuRftData)
+    {
+        RifReaderRftInterface* reader = m_observedFmuRftData->rftReader();
+        if (reader)
+        {
+            std::vector<double>    tvdValuesOfObbservedData;
+            std::vector<double>    mdValuesOfObbservedData;
+
+            RifEclipseRftAddress tvdAddress(m_wellName(), m_timeStep, RifEclipseRftAddress::TVD);
+            RifEclipseRftAddress mdAddress(m_wellName(), m_timeStep, RifEclipseRftAddress::MD);
+
+            reader->values(tvdAddress, &tvdValuesOfObbservedData);
+            reader->values(mdAddress, &mdValuesOfObbservedData);
+            interpolatedMdValues = RigWellPathGeometryTools::interpolateMdFromTvd(mdValuesOfObbservedData, tvdValuesOfObbservedData, tvDepthValues);
+            CVF_ASSERT(interpolatedMdValues.size() == tvDepthValues.size());
+        }
+    }
+    return interpolatedMdValues;
+}

ApplicationCode/ProjectDataModel/RimWellLogRftCurve.h

@@ -96,6 +96,7 @@ private:
     std::vector<double>  xValues();
     std::vector<double>  tvDepthValues();
     std::vector<double>  measuredDepthValues();
+    std::vector<double>  interpolatedMeasuredDepthValuesFromObservedData(const std::vector<double>& tvDepthValues);
 
 private:
     caf::PdmPtrField<RimEclipseResultCase*>     m_eclipseResultCase;
@@ -109,6 +110,7 @@ private:
 
     std::map<size_t, size_t>                                                 m_idxInWellPathToIdxInRftFile;
     bool                                                                     m_isUsingPseudoLength;
+    bool                                                                     m_derivingMDFromObservedData;
     caf::PdmField<caf::AppEnum<RifEclipseRftAddress::RftWellLogChannelType>> m_wellLogChannelName;
 
 };

ApplicationCode/ReservoirDataModel/RigWellPathGeometryTools.cpp

@@ -79,6 +79,69 @@ std::vector<cvf::Vec3d> RigWellPathGeometryTools::calculateLineSegmentNormals(co
     return interpolateUndefinedNormals(up, pointNormals, vertices);
 }
 
+//--------------------------------------------------------------------------------------------------
+/// Lets you estimate MD values from an existing md/tvd relationship and a new set of TVD-values
+/// Requires the points to be ordered from the start/top of the well path to the end/bottom.
+//--------------------------------------------------------------------------------------------------
+std::vector<double> RigWellPathGeometryTools::interpolateMdFromTvd(const std::vector<double>& originalMdValues, const std::vector<double>& originalTvdValues, const std::vector<double>& tvdValuesToInterpolateFrom)
+{
+    CVF_ASSERT(!originalMdValues.empty());
+    if (originalMdValues.size() < 2u)
+    {
+        return {originalMdValues};
+    }
+
+    std::vector<double> interpolatedMdValues;
+    interpolatedMdValues.reserve(tvdValuesToInterpolateFrom.size());
+
+    std::vector<double>::const_iterator last_it = originalTvdValues.begin();
+    for (std::vector<double>::const_iterator it = tvdValuesToInterpolateFrom.begin(); it != tvdValuesToInterpolateFrom.end(); ++it)
+    {
+        double tvdValue = *it;
+        double sign = 0.0;
+        if (it != tvdValuesToInterpolateFrom.begin())
+        {
+            sign = *it - *(it - 1);
+        }
+        else
+        {
+            sign = *(it + 1) - *it;
+        }
+        if (std::fabs(sign) < 1.0e-8)
+        {
+            continue;
+        }
+            
+        sign /= std::fabs(sign);
+
+        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,
                                                                               const std::vector<cvf::Vec3d>& normals,
                                                                               const std::vector<cvf::Vec3d>& vertices)
@@ -160,3 +223,21 @@ cvf::Vec3d RigWellPathGeometryTools::estimateDominantDirectionInXYPlane(const st
 
     return directionSum.getNormalized();
 }
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+double RigWellPathGeometryTools::linearInterpolation(const std::vector<double>& xValues, const std::vector<double>& yValues, int valueIndex, double x)
+{
+    int N = (int) xValues.size() - 1;
+    int i = cvf::Math::clamp(valueIndex, 1, N);
+        
+    double x1 = xValues[i - 1];
+    double x2 = xValues[i];
+    double y1 = yValues[i - 1];
+    double y2 = yValues[i];
+
+    double M = (y2 - y1) / (x2 - x1);
+
+    return M * (x - x1) + y1;
+}

ApplicationCode/ReservoirDataModel/RigWellPathGeometryTools.h

@@ -40,9 +40,14 @@ public:
 public:
     static std::vector<cvf::Vec3d> calculateLineSegmentNormals(const std::vector<cvf::Vec3d>& vertices,
                                                                double                         angle);
+    static std::vector<double> interpolateMdFromTvd(const std::vector<double>& originalMdValues, const std::vector<double>& originalTvdValues, const std::vector<double>& tvdValuesToInterpolateFrom);
+
 private:
     static std::vector<cvf::Vec3d> interpolateUndefinedNormals(const cvf::Vec3d& planeNormal,
                                                                const std::vector<cvf::Vec3d>& normals,
                                                                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);
+
 };

ApplicationCode/UnitTests/CMakeLists_files.cmake

@@ -42,6 +42,7 @@ ${CMAKE_CURRENT_LIST_DIR}/RiaEclipseUnitTools-Test.cpp
 ${CMAKE_CURRENT_LIST_DIR}/RiaTextFileCompare-Test.cpp
 ${CMAKE_CURRENT_LIST_DIR}/RifCaseRealizationParametersReader-Test.cpp
 ${CMAKE_CURRENT_LIST_DIR}/RigWellLogExtractor-Test.cpp
+${CMAKE_CURRENT_LIST_DIR}/RigWellPathGeometryTools-Test.cpp
 ${CMAKE_CURRENT_LIST_DIR}/RifEclipseSummaryAddress-Test.cpp
 ${CMAKE_CURRENT_LIST_DIR}/RiaTimeHistoryCurveTools-Test.cpp
 ${CMAKE_CURRENT_LIST_DIR}/SolveSpaceSolver-Test.cpp

ApplicationCode/UnitTests/RigWellPathGeometryTools-Test.cpp

@@ -0,0 +1,108 @@
+/////////////////////////////////////////////////////////////////////////////////
+//
+//  Copyright (C) 2019-     Equinor ASA
+//
+//  ResInsight is free software: you can redistribute it and/or modify
+//  it under the terms of the GNU General Public License as published by
+//  the Free Software Foundation, either version 3 of the License, or
+//  (at your option) any later version.
+//
+//  ResInsight is distributed in the hope that it will be useful, but WITHOUT ANY
+//  WARRANTY; without even the implied warranty of MERCHANTABILITY or
+//  FITNESS FOR A PARTICULAR PURPOSE.
+//
+//  See the GNU General Public License at <http://www.gnu.org/licenses/gpl.html>
+//  for more details.
+//
+/////////////////////////////////////////////////////////////////////////////////
+
+#include "gtest/gtest.h"
+
+#include "RigWellPathGeometryTools.h"
+
+#include <vector>
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+TEST(RigWellPathGeometryTools, VerticalPath)
+{
+    std::vector<double> mdValues  = {100, 500, 1000};
+    std::vector<double> tvdValues = {100, 500, 1000};
+    std::vector<double> fullTVDValues = {100, 200, 300, 400, 500, 600, 700, 800, 900, 1000};
+    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(fullTVDValues[i], fullMDValues[i]);
+    }
+}
+
+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]);
+    }
+}
+
+double quadraticFunction(double x)
+{
+    return 0.0015 * std::pow(x, 2) - 0.25 * x + 100;
+}
+
+TEST(RigWellPathGeometryTools, QuadraticPath)
+{
+    std::vector<double> mdValues = {100, 300, 600, 1000};
+    std::vector<double> tvdValues;
+    for (double md : mdValues)
+    {
+        tvdValues.push_back(quadraticFunction(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(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]);
+    }
+}
+
+double cubicFunction(double x)
+{
+    return 0.000012 * std::pow(x, 3) - 0.0175 * std::pow(x, 2) + 7 * x;
+}
+
+TEST(RigWellPathGeometryTools, CubicPath)
+{
+    std::vector<double> mdValues = {100, 300, 600, 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)
+    {
+        EXPECT_DOUBLE_EQ(fullMDValues[i], estimatedFullMDValues[i]);
+    }
+}