diff --git a/ApplicationCode/ProjectDataModel/Flow/RimWellRftPlot.cpp b/ApplicationCode/ProjectDataModel/Flow/RimWellRftPlot.cpp index 1d6f795075..4be3a3efbd 100644 --- a/ApplicationCode/ProjectDataModel/Flow/RimWellRftPlot.cpp +++ b/ApplicationCode/ProjectDataModel/Flow/RimWellRftPlot.cpp @@ -450,6 +450,7 @@ void RimWellRftPlot::updateCurvesInPlot(const std::set 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 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); @@ -731,7 +733,7 @@ QList RimWellRftPlot::calculateValueOptions(const caf::P item.setLevel(1); options.push_back(item); } - const std::vector observedFmuRftCases = RimWellPlotTools::observedFmuRftDataForWell(m_wellPathNameOrSimWellName); + const std::vector observedFmuRftCases = RimWellPlotTools::observedFmuRftDataForWell(m_wellPathNameOrSimWellName); if (!observedFmuRftCases.empty()) { options.push_back(caf::PdmOptionItemInfo::createHeader( @@ -1124,3 +1126,19 @@ void RimWellRftPlot::defineCurveColorsAndSymbols(const std::setrftReader()->availableTimeSteps(wellPathName).count(timeStep)) + { + return observedData; + } + } + return nullptr; +} diff --git a/ApplicationCode/ProjectDataModel/Flow/RimWellRftPlot.h b/ApplicationCode/ProjectDataModel/Flow/RimWellRftPlot.h index 591cf9673f..c3acfed6f6 100644 --- a/ApplicationCode/ProjectDataModel/Flow/RimWellRftPlot.h +++ b/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 selectedCurveDefs() const; std::set curveDefsFromCurves() const; @@ -161,5 +163,4 @@ private: std::map m_dataSourceColors; std::map m_timeStepSymbols; bool m_isOnLoad; - }; diff --git a/ApplicationCode/ProjectDataModel/RimObservedFmuRftData.cpp b/ApplicationCode/ProjectDataModel/RimObservedFmuRftData.cpp index 9d507469da..0c20313679 100644 --- a/ApplicationCode/ProjectDataModel/RimObservedFmuRftData.cpp +++ b/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 allWells = wells(); for (const QString& well : allWells) { - if (well == simWellName) + if (well == wellPathName) { return true; } diff --git a/ApplicationCode/ProjectDataModel/RimObservedFmuRftData.h b/ApplicationCode/ProjectDataModel/RimObservedFmuRftData.h index 8aafd09484..ecdb891bb0 100644 --- a/ApplicationCode/ProjectDataModel/RimObservedFmuRftData.h +++ b/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 wells() const; + std::vector> derivedWellPathTvdMd(const QString& wellPathName) const; + private: cvf::ref m_fmuRftReader; diff --git a/ApplicationCode/ProjectDataModel/RimWellLogRftCurve.cpp b/ApplicationCode/ProjectDataModel/RimWellLogRftCurve.cpp index c4bef49bfb..03858fa329 100644 --- a/ApplicationCode/ProjectDataModel/RimWellLogRftCurve.cpp +++ b/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(m_curveData->xPlotValues().size())); m_isUsingPseudoLength = false; + m_derivingMDFromObservedData = false; } - if (m_isUsingPseudoLength) + if (m_isUsingPseudoLength || m_derivingMDFromObservedData) { RimWellLogTrack* wellLogTrack; firstAncestorOrThisOfType(wellLogTrack); @@ -443,7 +455,10 @@ void RimWellLogRftCurve::onLoadDataAndUpdate(bool updateParentPlot) RiuWellLogTrack* viewer = wellLogTrack->viewer(); if (viewer) { - viewer->setDepthTitle("PL/" + wellLogPlot->depthPlotTitle()); + if (m_derivingMDFromObservedData) + viewer->setDepthTitle("OBS/" + wellLogPlot->depthPlotTitle()); + else + viewer->setDepthTitle("PL/" + wellLogPlot->depthPlotTitle()); } } @@ -835,7 +850,7 @@ std::vector RimWellLogRftCurve::tvDepthValues() //-------------------------------------------------------------------------------------------------- std::vector RimWellLogRftCurve::measuredDepthValues() { - if (m_observedFmuRftData) + if (m_observedFmuRftData && !m_ensemble && !m_summaryCase) { RifReaderRftInterface* reader = rftReader(); std::vector values; @@ -881,3 +896,29 @@ std::vector RimWellLogRftCurve::measuredDepthValues() return measuredDepthForCellsWhichHasRftData; } + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +std::vector RimWellLogRftCurve::interpolatedMeasuredDepthValuesFromObservedData(const std::vector& tvDepthValues) +{ + std::vector interpolatedMdValues; + if (m_observedFmuRftData) + { + RifReaderRftInterface* reader = m_observedFmuRftData->rftReader(); + if (reader) + { + std::vector tvdValuesOfObbservedData; + std::vector 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; +} diff --git a/ApplicationCode/ProjectDataModel/RimWellLogRftCurve.h b/ApplicationCode/ProjectDataModel/RimWellLogRftCurve.h index 3252aa3557..d135fe4235 100644 --- a/ApplicationCode/ProjectDataModel/RimWellLogRftCurve.h +++ b/ApplicationCode/ProjectDataModel/RimWellLogRftCurve.h @@ -96,6 +96,7 @@ private: std::vector xValues(); std::vector tvDepthValues(); std::vector measuredDepthValues(); + std::vector interpolatedMeasuredDepthValuesFromObservedData(const std::vector& tvDepthValues); private: caf::PdmPtrField m_eclipseResultCase; @@ -109,6 +110,7 @@ private: std::map m_idxInWellPathToIdxInRftFile; bool m_isUsingPseudoLength; + bool m_derivingMDFromObservedData; caf::PdmField> m_wellLogChannelName; }; diff --git a/ApplicationCode/ReservoirDataModel/RigWellPathGeometryTools.cpp b/ApplicationCode/ReservoirDataModel/RigWellPathGeometryTools.cpp index 188ec95d92..0093af617e 100644 --- a/ApplicationCode/ReservoirDataModel/RigWellPathGeometryTools.cpp +++ b/ApplicationCode/ReservoirDataModel/RigWellPathGeometryTools.cpp @@ -79,6 +79,69 @@ std::vector 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 RigWellPathGeometryTools::interpolateMdFromTvd(const std::vector& originalMdValues, const std::vector& originalTvdValues, const std::vector& tvdValuesToInterpolateFrom) +{ + CVF_ASSERT(!originalMdValues.empty()); + if (originalMdValues.size() < 2u) + { + return {originalMdValues}; + } + + std::vector interpolatedMdValues; + interpolatedMdValues.reserve(tvdValuesToInterpolateFrom.size()); + + std::vector::const_iterator last_it = originalTvdValues.begin(); + for (std::vector::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(current_it - originalTvdValues.begin()); + double mdValue = linearInterpolation(originalTvdValues, originalMdValues, valueIndex, tvdValue); + interpolatedMdValues.push_back(mdValue); + last_it = current_it; + } + return interpolatedMdValues; +} + std::vector RigWellPathGeometryTools::interpolateUndefinedNormals(const cvf::Vec3d& planeNormal, const std::vector& normals, const std::vector& vertices) @@ -160,3 +223,21 @@ cvf::Vec3d RigWellPathGeometryTools::estimateDominantDirectionInXYPlane(const st return directionSum.getNormalized(); } + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +double RigWellPathGeometryTools::linearInterpolation(const std::vector& xValues, const std::vector& 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; +} diff --git a/ApplicationCode/ReservoirDataModel/RigWellPathGeometryTools.h b/ApplicationCode/ReservoirDataModel/RigWellPathGeometryTools.h index 838ba966ec..b2be6dccc3 100644 --- a/ApplicationCode/ReservoirDataModel/RigWellPathGeometryTools.h +++ b/ApplicationCode/ReservoirDataModel/RigWellPathGeometryTools.h @@ -40,9 +40,14 @@ public: public: static std::vector calculateLineSegmentNormals(const std::vector& vertices, double angle); + static std::vector interpolateMdFromTvd(const std::vector& originalMdValues, const std::vector& originalTvdValues, const std::vector& tvdValuesToInterpolateFrom); + private: static std::vector interpolateUndefinedNormals(const cvf::Vec3d& planeNormal, const std::vector& normals, const std::vector& vertices); static cvf::Vec3d estimateDominantDirectionInXYPlane(const std::vector& vertices); + + static double linearInterpolation(const std::vector& xValues, const std::vector& yValues, int valueIndex, double x); + }; diff --git a/ApplicationCode/UnitTests/CMakeLists_files.cmake b/ApplicationCode/UnitTests/CMakeLists_files.cmake index 6b60930eee..ba3413d9a2 100644 --- a/ApplicationCode/UnitTests/CMakeLists_files.cmake +++ b/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 diff --git a/ApplicationCode/UnitTests/RigWellPathGeometryTools-Test.cpp b/ApplicationCode/UnitTests/RigWellPathGeometryTools-Test.cpp new file mode 100644 index 0000000000..084d8864ce --- /dev/null +++ b/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 +// for more details. +// +///////////////////////////////////////////////////////////////////////////////// + +#include "gtest/gtest.h" + +#include "RigWellPathGeometryTools.h" + +#include + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +TEST(RigWellPathGeometryTools, VerticalPath) +{ + std::vector mdValues = {100, 500, 1000}; + std::vector tvdValues = {100, 500, 1000}; + std::vector fullTVDValues = {100, 200, 300, 400, 500, 600, 700, 800, 900, 1000}; + std::vector 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 mdValues = {100, 500, 1000}; + std::vector tvdValues = {50, 250, 500}; + std::vector fullTVDValues = {50, 100, 150, 200, 250, 300, 350, 400, 450, 500}; + std::vector 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 mdValues = {100, 300, 600, 1000}; + std::vector tvdValues; + for (double md : mdValues) + { + tvdValues.push_back(quadraticFunction(md)); + } + std::vector fullMDValues = {100, 200, 300, 400, 500, 600, 700, 800, 900, 1000}; + std::vector fullTvdValues; + for (double md : fullMDValues) + { + fullTvdValues.push_back(quadraticFunction(md)); + } + std::vector 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 mdValues = {100, 300, 600, 1000}; + std::vector tvdValues; + for (double md : mdValues) + { + tvdValues.push_back(cubicFunction(md)); + } + std::vector fullMDValues = {100, 200, 300, 400, 500, 600, 700, 800, 900, 1000}; + std::vector fullTvdValues; + for (double md : fullMDValues) + { + fullTvdValues.push_back(cubicFunction(md)); + } + std::vector 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]); + } +} \ No newline at end of file