diff --git a/ApplicationLibCode/GeoMech/GeoMechDataModel/CMakeLists.txt b/ApplicationLibCode/GeoMech/GeoMechDataModel/CMakeLists.txt index 4619c4a8d7..db22d5201a 100644 --- a/ApplicationLibCode/GeoMech/GeoMechDataModel/CMakeLists.txt +++ b/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 ) diff --git a/ApplicationLibCode/GeoMech/GeoMechDataModel/RigFemPartResultCalculatorKIndices.cpp b/ApplicationLibCode/GeoMech/GeoMechDataModel/RigFemPartResultCalculatorKIndices.cpp new file mode 100644 index 0000000000..25c00127f9 --- /dev/null +++ b/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 +// 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 + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +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& dstFrameData = resFrames->frameData( 0 ); + + const size_t valCount = femPart->elementNodeResultCount(); + dstFrameData.resize( valCount, std::numeric_limits::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::infinity() returns 0 + constexpr size_t maxValue = std::numeric_limits::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( k ) : HUGE_VALF; + } + } + } + } + + return resFrames; +} diff --git a/ApplicationLibCode/GeoMech/GeoMechDataModel/RigFemPartResultCalculatorKIndices.h b/ApplicationLibCode/GeoMech/GeoMechDataModel/RigFemPartResultCalculatorKIndices.h new file mode 100644 index 0000000000..a73397d97f --- /dev/null +++ b/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 +// 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; +}; diff --git a/ApplicationLibCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.cpp b/ApplicationLibCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.cpp index 2a155229ef..4bfc28c6ca 100644 --- a/ApplicationLibCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.cpp +++ b/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( new RigFemPartResultCalculatorShearSlipIndicator( *this ) ) ); m_resultCalculators.push_back( std::unique_ptr( new RigFemPartResultCalculatorFormationIndices( *this ) ) ); + m_resultCalculators.push_back( + std::unique_ptr( new RigFemPartResultCalculatorKIndices( *this ) ) ); } //-------------------------------------------------------------------------------------------------- @@ -684,6 +687,8 @@ std::map> fieldCompNames["PE"].push_back( "PE2" ); fieldCompNames["PE"].push_back( "PE3" ); } + + fieldCompNames["INDEX"].push_back( "INDEX_K" ); } else if ( resPos == RIG_INTEGRATION_POINT ) { diff --git a/ApplicationLibCode/ProjectDataModel/WellLog/RimWellLogExtractionCurve.cpp b/ApplicationLibCode/ProjectDataModel/WellLog/RimWellLogExtractionCurve.cpp index 4e078d21ab..5d03bcfb2b 100644 --- a/ApplicationLibCode/ProjectDataModel/WellLog/RimWellLogExtractionCurve.cpp +++ b/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 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? - // ************************************************ - // - // Adjust measured dept values according to refWell - // - // ************************************************ + RigFemResultAddress indexKResAdr( RigFemResultPosEnum::RIG_ELEMENT_NODAL, "INDEX", "INDEX_K" ); + + const size_t frameIdx = 0; + std::vector refWellMeasuredDepthValues = refWellExtractor->cellIntersectionMDs(); + std::vector refWellTvDepthValues = refWellExtractor->cellIntersectionTVDs(); + std::vector wellIndexKValues; + std::vector 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 -// enter/exit values and linearize between enter/exit of k-index +/// Utility function to adjust well depth values according to reference well - match +/// 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& rMeasuredDepthValues, std::vector& rTvDepthValues, @@ -666,12 +686,6 @@ void RimWellLogExtractionCurve::adjustWellDepthValuesToReferenceWell( std::vecto const std::vector& refWellTvDepthValues, const std::vector& 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 - const double minLayerK = std::max( *std::min_element( refWellIndexKValues.cbegin(), refWellIndexKValues.cend() ), - *std::min_element( indexKValues.cbegin(), indexKValues.cend() ) ); - const double maxLayerK = std::min( *std::max_element( refWellIndexKValues.cbegin(), refWellIndexKValues.cend() ), - *std::max_element( indexKValues.cbegin(), indexKValues.cend() ) ); + // Find common min and max k-index value for range of depth values to adjust + const auto minLayerK = + static_cast( std::max( *std::min_element( refWellIndexKValues.cbegin(), refWellIndexKValues.cend() ), + *std::min_element( indexKValues.cbegin(), indexKValues.cend() ) ) ); + const auto maxLayerK = + static_cast( 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; - for ( int kLayer = static_cast( minLayerK ); kLayer <= static_cast( maxLayerK ); kLayer++ ) + // Only allow asymptotically increasing k-layers - break at first decreasing k-layer value + auto createKLayerAndIndexMap = []( const std::vector& indexKValues, int minLayerK, int maxLayerK ) { + int prevKLayer = -1; + std::map> kLayerAndIndexesMap = {}; + for ( size_t i = 0; i < indexKValues.size(); ++i ) + { + const auto kLayer = static_cast( 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> refWellKLayerAndIndexesMap = + createKLayerAndIndexMap( refWellIndexKValues, minLayerK, maxLayerK ); + for ( const auto& [kLayer, indexes] : refWellKLayerAndIndexesMap ) { - const auto kLayerTopIter = - std::find( refWellIndexKValues.cbegin(), refWellIndexKValues.cend(), static_cast( kLayer ) ); - const auto kLayerBottomRIter = - std::find( refWellIndexKValues.crbegin(), refWellIndexKValues.crend(), static_cast( kLayer ) ); - const auto indexTop = std::distance( refWellIndexKValues.cbegin(), kLayerTopIter ); - const auto indexBottom = refWellIndexKValues.size() - 1 - - std::distance( refWellIndexKValues.crbegin(), kLayerBottomRIter ); - - 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> wellKLayerAndIndexesMap = {}; - for ( size_t i = 0; i < indexKValues.size(); i++ ) - { - const int kLayer = static_cast( indexKValues[i] ); - wellKLayerAndIndexesMap[kLayer].push_back( i ); + if ( indexes.empty() ) continue; + const auto indexTop = indexes.front(); + const auto indexBottom = indexes.back(); + refWellLogIndexDepthOffset.setIndexOffsetDepth( kLayer, + refWellMeasuredDepthValues[indexTop], + refWellMeasuredDepthValues[indexBottom], + refWellTvDepthValues[indexTop], + refWellTvDepthValues[indexBottom] ); } + std::map> 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 ) ) {