Add reference well depth adjustment to geomech case (#9547)

- Fix bug in retrieving bottom TVD value from ref well
- Add result calculator for K-indices in geomech case
- Improve depth adjustment algorithm

ApplicationLibCode/GeoMech/GeoMechDataModel/CMakeLists.txt

@@ -102,6 +102,8 @@ add_library(
   RigFemPartResultCalculatorMudWeightWindow.cpp
   RigFemPartResultCalculatorShearSlipIndicator.h
   RigFemPartResultCalculatorShearSlipIndicator.cpp
+  RigFemPartResultCalculatorKIndices.h
+  RigFemPartResultCalculatorKIndices.cpp
   RimGeoMechGeometrySelectionItem.h
   RimGeoMechGeometrySelectionItem.cpp
 )

ApplicationLibCode/GeoMech/GeoMechDataModel/RigFemPartResultCalculatorKIndices.cpp

@@ -0,0 +1,109 @@
+/////////////////////////////////////////////////////////////////////////////////
+//
+//  Copyright (C) 2022-     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 "RigFemPartResultCalculatorKIndices.h"
+
+#include "RigFemPart.h"
+#include "RigFemPartCollection.h"
+#include "RigFemPartGrid.h"
+#include "RigFemPartResultsCollection.h"
+#include "RigFemResultAddress.h"
+#include "RigFemScalarResultFrames.h"
+#include "RigFormationNames.h"
+
+#include "cafProgressInfo.h"
+
+#include <QString>
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+RigFemPartResultCalculatorKIndices::RigFemPartResultCalculatorKIndices( RigFemPartResultsCollection& collection )
+    : RigFemPartResultCalculator( collection )
+{
+}
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+RigFemPartResultCalculatorKIndices::~RigFemPartResultCalculatorKIndices()
+{
+}
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+bool RigFemPartResultCalculatorKIndices::isMatching( const RigFemResultAddress& resVarAddr ) const
+{
+    return ( resVarAddr.fieldName == "INDEX" && resVarAddr.componentName == "INDEX_K" );
+}
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+RigFemScalarResultFrames* RigFemPartResultCalculatorKIndices::calculate( int                        partIndex,
+                                                                         const RigFemResultAddress& resVarAddr )
+{
+    caf::ProgressInfo frameCountProgress( 2, "" );
+    frameCountProgress.setProgressDescription(
+        "Calculating " + QString::fromStdString( resVarAddr.fieldName + ": " + resVarAddr.componentName ) );
+
+    RigFemScalarResultFrames* resFrames = m_resultCollection->createScalarResult( partIndex, resVarAddr );
+    resFrames->enableAsSingleFrameResult();
+
+    const RigFemPart*   femPart      = m_resultCollection->parts()->part( partIndex );
+    std::vector<float>& dstFrameData = resFrames->frameData( 0 );
+
+    const size_t valCount = femPart->elementNodeResultCount();
+    dstFrameData.resize( valCount, std::numeric_limits<float>::infinity() );
+
+    const RigFormationNames* activeFormNames = m_resultCollection->activeFormationNames();
+
+    frameCountProgress.incrementProgress();
+
+    if ( activeFormNames )
+    {
+        // Has to be done before the parallel loop because the first call allocates.
+        const RigFemPartGrid* structGrid = femPart->getOrCreateStructGrid();
+
+        const int elementCount = femPart->elementCount();
+
+        // Using max() as std::numeric_limits<size_t>::infinity() returns 0
+        constexpr size_t maxValue = std::numeric_limits<size_t>::max();
+
+#pragma omp parallel for
+        for ( int elmIdx = 0; elmIdx < elementCount; ++elmIdx )
+        {
+            RigElementType elmType      = femPart->elementType( elmIdx );
+            int            elmNodeCount = RigFemTypes::elementNodeCount( elmType );
+
+            size_t i, j, k = maxValue;
+            bool   validIndex = structGrid->ijkFromCellIndex( elmIdx, &i, &j, &k );
+            if ( validIndex )
+            {
+                for ( int elmNodIdx = 0; elmNodIdx < elmNodeCount; ++elmNodIdx )
+                {
+                    size_t elmNodResIdx        = femPart->elementNodeResultIdx( elmIdx, elmNodIdx );
+                    dstFrameData[elmNodResIdx] = k != maxValue ? static_cast<float>( k ) : HUGE_VALF;
+                }
+            }
+        }
+    }
+
+    return resFrames;
+}

ApplicationLibCode/GeoMech/GeoMechDataModel/RigFemPartResultCalculatorKIndices.h

@@ -0,0 +1,37 @@
+/////////////////////////////////////////////////////////////////////////////////
+//
+//  Copyright (C) 2022-     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.
+//
+/////////////////////////////////////////////////////////////////////////////////
+
+#pragma once
+
+#include "RigFemPartResultCalculator.h"
+
+class RigFemPartResultsCollection;
+class RigFemScalarResultFrames;
+class RigFemResultAddress;
+
+//==================================================================================================
+///
+//==================================================================================================
+class RigFemPartResultCalculatorKIndices : public RigFemPartResultCalculator
+{
+public:
+    explicit RigFemPartResultCalculatorKIndices( RigFemPartResultsCollection& collection );
+    ~RigFemPartResultCalculatorKIndices() override;
+    bool                      isMatching( const RigFemResultAddress& resVarAddr ) const override;
+    RigFemScalarResultFrames* calculate( int partIndex, const RigFemResultAddress& resVarAddr ) override;
+};

ApplicationLibCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.cpp

@@ -37,6 +37,7 @@
 #include "RigFemPartResultCalculatorFormationIndices.h"
 #include "RigFemPartResultCalculatorGamma.h"
 #include "RigFemPartResultCalculatorInitialPorosity.h"
+#include "RigFemPartResultCalculatorKIndices.h"
 #include "RigFemPartResultCalculatorMudWeightWindow.h"
 #include "RigFemPartResultCalculatorNE.h"
 #include "RigFemPartResultCalculatorNodalGradients.h"
@@ -199,6 +200,8 @@ RigFemPartResultsCollection::RigFemPartResultsCollection( RifGeoMechReaderInterf
         std::unique_ptr<RigFemPartResultCalculator>( new RigFemPartResultCalculatorShearSlipIndicator( *this ) ) );
     m_resultCalculators.push_back(
         std::unique_ptr<RigFemPartResultCalculator>( new RigFemPartResultCalculatorFormationIndices( *this ) ) );
+    m_resultCalculators.push_back(
+        std::unique_ptr<RigFemPartResultCalculator>( new RigFemPartResultCalculatorKIndices( *this ) ) );
 }
 
 //--------------------------------------------------------------------------------------------------
@@ -684,6 +687,8 @@ std::map<std::string, std::vector<std::string>>
                 fieldCompNames["PE"].push_back( "PE2" );
                 fieldCompNames["PE"].push_back( "PE3" );
             }
+
+            fieldCompNames["INDEX"].push_back( "INDEX_K" );
         }
         else if ( resPos == RIG_INTEGRATION_POINT )
         {

ApplicationLibCode/ProjectDataModel/WellLog/RimWellLogExtractionCurve.cpp

@@ -552,12 +552,6 @@ RimWellLogExtractionCurve::WellLogExtractionCurveData
     // Reference well adjustment does not support simulated wells
     if ( m_trajectoryType == WELL_PATH && wellExtractor.notNull() && refWellExtractor.notNull() )
     {
-        // ************************************************
-        //
-        // Adjust measured dept values according to refWell
-        //
-        // ************************************************
-
         RigEclipseResultAddress indexKResAdr( RiaDefines::ResultCatType::STATIC_NATIVE,
                                               RiaResultNames::indexKResultName() );
         eclipseCase->eclipseCaseData()->results( RiaDefines::PorosityModelType::MATRIX_MODEL )->ensureKnownResultLoaded( indexKResAdr );
@@ -601,8 +595,6 @@ RimWellLogExtractionCurve::WellLogExtractionCurveData
                                                 bool            performDataSmoothing,
                                                 double          smoothingThreshold )
 {
-    // TODO: Add depth adjustements for reference well
-
     WellLogExtractionCurveData           curveData;
     RimWellLogPlotCollection*            wellLogCollection = RimMainPlotCollection::current()->wellLogPlotCollection();
     cvf::ref<RigGeoMechWellLogExtractor> wellExtractor = wellLogCollection->findOrCreateExtractor( m_wellPath, geomCase );
@@ -641,23 +633,51 @@ RimWellLogExtractionCurve::WellLogExtractionCurveData
         }
     }
 
+    const std::string fieldName = m_geomResultDefinition->resultAddress().fieldName;
+    if ( fieldName == RiaResultNames::wbsAzimuthResult().toStdString() ||
+         fieldName == RiaResultNames::wbsInclinationResult().toStdString() )
+    {
+        // Do not adjust depth values of Azimuth and Inclination as they are
+        // dependent of the well geometry
+        return curveData;
+    }
+
     if ( wellExtractor.notNull() && refWellExtractor.notNull() )
     {
-        // Add reference well depth adjustments
+        // NOTE: Perform performCurveDataSmoothing() for refWell, as done for well?
 
-        // ************************************************
+        RigFemResultAddress indexKResAdr( RigFemResultPosEnum::RIG_ELEMENT_NODAL, "INDEX", "INDEX_K" );
-        //
+
-        // Adjust measured dept values according to refWell
+        const size_t        frameIdx                   = 0;
-        //
+        std::vector<double> refWellMeasuredDepthValues = refWellExtractor->cellIntersectionMDs();
-        // ************************************************
+        std::vector<double> refWellTvDepthValues       = refWellExtractor->cellIntersectionTVDs();
+        std::vector<double> wellIndexKValues;
+        std::vector<double> refWellIndexKValues;
+        if ( indexKResAdr.isValid() )
+        {
+            wellExtractor->curveData( indexKResAdr, frameIdx, &wellIndexKValues );
+            refWellExtractor->curveData( indexKResAdr, frameIdx, &refWellIndexKValues );
+        }
+        if ( !wellIndexKValues.empty() && !refWellIndexKValues.empty() )
+        {
+            adjustWellDepthValuesToReferenceWell( curveData.measuredDepthValues,
+                                                  curveData.tvDepthValues,
+                                                  wellIndexKValues,
+                                                  refWellMeasuredDepthValues,
+                                                  refWellTvDepthValues,
+                                                  refWellIndexKValues );
+        }
     }
 
     return curveData;
 }
 
 //--------------------------------------------------------------------------------------------------
-// Utility function to adjust well depth values according to reference well - match k-index
+/// Utility function to adjust well depth values according to reference well - match
-//     enter/exit values and linearize between enter/exit of k-index
+/// enter and exit values for k-layer and linearize for depth values between enter/exit of k-layer.
+///
+/// Performs depth adjustment from first common k-layer to last common k-layer - i.e common min and
+/// max index-k, as long as the k-layers are asymptotically increasing.
 //--------------------------------------------------------------------------------------------------
 void RimWellLogExtractionCurve::adjustWellDepthValuesToReferenceWell( std::vector<double>&       rMeasuredDepthValues,
                                                                       std::vector<double>&       rTvDepthValues,
@@ -666,12 +686,6 @@ void RimWellLogExtractionCurve::adjustWellDepthValuesToReferenceWell( std::vecto
                                                                       const std::vector<double>& refWellTvDepthValues,
                                                                       const std::vector<double>& refWellIndexKValues )
 {
-    // Description:
-    // - Adjust values up to largest common k-index value
-    // Assumptions:
-    // - Both wells have min k-index equal to 1.
-    // - k-index values are continously increasing values between top and bottom of well
-
     CAF_ASSERT( rMeasuredDepthValues.size() == rTvDepthValues.size() &&
                 "Number of depth values must be equal for well!" );
     CAF_ASSERT( rMeasuredDepthValues.size() == indexKValues.size() &&
@@ -682,13 +696,15 @@ void RimWellLogExtractionCurve::adjustWellDepthValuesToReferenceWell( std::vecto
                 "Number of index K values must be number of depth values for reference well!" );
     CAF_ASSERT( *std::min( indexKValues.cbegin(), indexKValues.cend() ) ==
                     *std::min( refWellIndexKValues.cbegin(), refWellIndexKValues.cend() ) &&
-                "Both index-K value vectors must have same start index-K layer" );
+                "Both index-K value vectors must containt common min index-K layer" );
 
-    // Find common min and max k-index value for range depth values adjustment
+    // Find common min and max k-index value for range of depth values to adjust
-    const double minLayerK = std::max( *std::min_element( refWellIndexKValues.cbegin(), refWellIndexKValues.cend() ),
+    const auto minLayerK =
-                                       *std::min_element( indexKValues.cbegin(), indexKValues.cend() ) );
+        static_cast<int>( std::max( *std::min_element( refWellIndexKValues.cbegin(), refWellIndexKValues.cend() ),
-    const double maxLayerK = std::min( *std::max_element( refWellIndexKValues.cbegin(), refWellIndexKValues.cend() ),
+                                    *std::min_element( indexKValues.cbegin(), indexKValues.cend() ) ) );
-                                       *std::max_element( indexKValues.cbegin(), indexKValues.cend() ) );
+    const auto maxLayerK =
+        static_cast<int>( std::min( *std::max_element( refWellIndexKValues.cbegin(), refWellIndexKValues.cend() ),
+                                    *std::max_element( indexKValues.cbegin(), indexKValues.cend() ) ) );
     if ( minLayerK > maxLayerK )
     {
         RiaLogging::error(
@@ -696,32 +712,39 @@ void RimWellLogExtractionCurve::adjustWellDepthValuesToReferenceWell( std::vecto
         return;
     }
 
-    RigWellLogIndexDepthOffset refWellLogIndexDepthOffset;
+    // Only allow asymptotically increasing k-layers - break at first decreasing k-layer value
-    for ( int kLayer = static_cast<int>( minLayerK ); kLayer <= static_cast<int>( maxLayerK ); kLayer++ )
+    auto createKLayerAndIndexMap = []( const std::vector<double>& indexKValues, int minLayerK, int maxLayerK ) {
+        int                                prevKLayer          = -1;
+        std::map<int, std::vector<size_t>> kLayerAndIndexesMap = {};
+        for ( size_t i = 0; i < indexKValues.size(); ++i )
+        {
+            const auto kLayer = static_cast<int>( indexKValues[i] );
+            if ( kLayer < minLayerK ) continue;
+            if ( kLayer < prevKLayer || kLayer > maxLayerK ) break;
+
+            kLayerAndIndexesMap[kLayer].push_back( i );
+            prevKLayer = kLayer;
+        }
+        return kLayerAndIndexesMap;
+    };
+
+    RigWellLogIndexDepthOffset         refWellLogIndexDepthOffset;
+    std::map<int, std::vector<size_t>> refWellKLayerAndIndexesMap =
+        createKLayerAndIndexMap( refWellIndexKValues, minLayerK, maxLayerK );
+    for ( const auto& [kLayer, indexes] : refWellKLayerAndIndexesMap )
     {
-        const auto kLayerTopIter =
+        if ( indexes.empty() ) continue;
-            std::find( refWellIndexKValues.cbegin(), refWellIndexKValues.cend(), static_cast<double>( kLayer ) );
+        const auto indexTop    = indexes.front();
-        const auto kLayerBottomRIter =
+        const auto indexBottom = indexes.back();
-            std::find( refWellIndexKValues.crbegin(), refWellIndexKValues.crend(), static_cast<double>( kLayer ) );
+        refWellLogIndexDepthOffset.setIndexOffsetDepth( kLayer,
-        const auto indexTop    = std::distance( refWellIndexKValues.cbegin(), kLayerTopIter );
+                                                        refWellMeasuredDepthValues[indexTop],
-        const auto indexBottom = refWellIndexKValues.size() - 1 -
+                                                        refWellMeasuredDepthValues[indexBottom],
-                                 std::distance( refWellIndexKValues.crbegin(), kLayerBottomRIter );
+                                                        refWellTvDepthValues[indexTop],
-
+                                                        refWellTvDepthValues[indexBottom] );
-        const auto topMd     = refWellMeasuredDepthValues[indexTop];
-        const auto bottomMd  = refWellMeasuredDepthValues[indexBottom];
-        const auto topTvd    = refWellTvDepthValues[indexTop];
-        const auto bottomTvd = refWellTvDepthValues[indexBottom];
-
-        refWellLogIndexDepthOffset.setIndexOffsetDepth( kLayer, topMd, bottomMd, topTvd, bottomTvd );
-    }
-
-    std::map<int, std::vector<size_t>> wellKLayerAndIndexesMap = {};
-    for ( size_t i = 0; i < indexKValues.size(); i++ )
-    {
-        const int kLayer = static_cast<int>( indexKValues[i] );
-        wellKLayerAndIndexesMap[kLayer].push_back( i );
     }
 
+    std::map<int, std::vector<size_t>> wellKLayerAndIndexesMap =
+        createKLayerAndIndexMap( indexKValues, minLayerK, maxLayerK );
     for ( const auto& [kLayer, indexes] : wellKLayerAndIndexesMap )
     {
         const auto firstIdx = indexes.front();
@@ -731,7 +754,7 @@ void RimWellLogExtractionCurve::adjustWellDepthValuesToReferenceWell( std::vecto
             rMeasuredDepthValues[firstIdx] = refWellLogIndexDepthOffset.getTopMd( kLayer );
             rMeasuredDepthValues[lastIdx]  = refWellLogIndexDepthOffset.getBottomMd( kLayer );
             rTvDepthValues[firstIdx]       = refWellLogIndexDepthOffset.getTopTvd( kLayer );
-            rTvDepthValues[lastIdx]        = refWellLogIndexDepthOffset.getBottomMd( kLayer );
+            rTvDepthValues[lastIdx]        = refWellLogIndexDepthOffset.getBottomTvd( kLayer );
         }
         else if ( indexes.size() > 2 && refWellLogIndexDepthOffset.hasIndex( kLayer ) )
         {