From ea496072b87e66a54ba7db0dc5bec19925010e3a Mon Sep 17 00:00:00 2001 From: Kristian Bendiksen Date: Thu, 21 Oct 2021 09:55:33 +0200 Subject: [PATCH] #8102 Ensemble Well Log Plot: Fix crash when using TVD --- .../WellLog/RimEnsembleWellLogCurveSet.cpp | 23 +++-- .../WellLog/RimEnsembleWellLogCurveSet.h | 4 +- .../WellLog/RimEnsembleWellLogStatistics.cpp | 90 +++++++++++++------ .../WellLog/RimEnsembleWellLogStatistics.h | 4 +- .../RimEnsembleWellLogStatisticsCurve.cpp | 29 ++++-- .../WellLog/RimWellLogFileCurve.cpp | 4 +- .../RigWellLogIndexDepthOffset.cpp | 47 ++++++++-- .../RigWellLogIndexDepthOffset.h | 11 ++- 8 files changed, 150 insertions(+), 62 deletions(-) diff --git a/ApplicationLibCode/ProjectDataModel/WellLog/RimEnsembleWellLogCurveSet.cpp b/ApplicationLibCode/ProjectDataModel/WellLog/RimEnsembleWellLogCurveSet.cpp index b3048b8db8..40c66b5f13 100644 --- a/ApplicationLibCode/ProjectDataModel/WellLog/RimEnsembleWellLogCurveSet.cpp +++ b/ApplicationLibCode/ProjectDataModel/WellLog/RimEnsembleWellLogCurveSet.cpp @@ -855,21 +855,24 @@ void RimEnsembleWellLogCurveSet::setLogScaleFromSelectedResult( const QString re //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- -void RimEnsembleWellLogCurveSet::updateStatistics() +bool RimEnsembleWellLogCurveSet::updateStatistics() { - updateStatistics( std::vector() ); + return updateStatistics( std::vector() ); } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- -void RimEnsembleWellLogCurveSet::updateStatistics( const std::vector& sumCases ) +bool RimEnsembleWellLogCurveSet::updateStatistics( const std::vector& sumCases ) { RimEnsembleWellLogs* ensembleWellLogs = m_ensembleWellLogs; QString wellLogChannelName = m_wellLogChannelName(); - if ( !isCurvesVisible() || m_disableStatisticCurves || !ensembleWellLogs || wellLogChannelName == "None" ) return; - + if ( !isCurvesVisible() || m_disableStatisticCurves || !ensembleWellLogs || wellLogChannelName == "None" ) + { + m_ensembleWellLogStatistics->clearData(); + return false; + } // Calculate std::vector wellLogFiles = sumCases; if ( wellLogFiles.empty() ) @@ -881,6 +884,7 @@ void RimEnsembleWellLogCurveSet::updateStatistics( const std::vectorcalculate( wellLogFiles, wellLogChannelName, m_depthEqualization() ); + return true; } //-------------------------------------------------------------------------------------------------- @@ -888,7 +892,9 @@ void RimEnsembleWellLogCurveSet::updateStatistics( const std::vector& sumCases ) { - updateStatistics( sumCases ); + deleteStatisticsCurves(); + + if ( !updateStatistics( sumCases ) ) return; RimWellLogPlot* plot = nullptr; firstAncestorOrThisOfType( plot ); @@ -919,7 +925,6 @@ void RimEnsembleWellLogCurveSet::updateStatisticsCurves( const std::vector for ( int kLayer : offsets->sortedIndexes() ) { RigWellPathFormation wellPathFormation; - wellPathFormation.mdTop = offsets->getTopDepth( kLayer ); - wellPathFormation.mdBase = offsets->getBottomDepth( kLayer ); + wellPathFormation.mdTop = offsets->getTopMd( kLayer ); + wellPathFormation.mdBase = offsets->getBottomMd( kLayer ); wellPathFormation.formationName = formationNames->formationNameFromKLayerIdx( kLayer - 1 ); wellPathFormationItems.push_back( wellPathFormation ); } diff --git a/ApplicationLibCode/ProjectDataModel/WellLog/RimEnsembleWellLogCurveSet.h b/ApplicationLibCode/ProjectDataModel/WellLog/RimEnsembleWellLogCurveSet.h index a78cef5dc6..adfff2d148 100644 --- a/ApplicationLibCode/ProjectDataModel/WellLog/RimEnsembleWellLogCurveSet.h +++ b/ApplicationLibCode/ProjectDataModel/WellLog/RimEnsembleWellLogCurveSet.h @@ -119,7 +119,7 @@ public: RimEnsembleWellLogStatistics::DepthEqualization depthEqualization() const; void setDepthEqualization( RimEnsembleWellLogStatistics::DepthEqualization depthEqualization ); - void updateStatistics(); + bool updateStatistics(); void setEnsembleWellLogs( RimEnsembleWellLogs* ensembleWellLogs ); void setWellLogChannelName( const QString& wellLogChannelName ); @@ -131,7 +131,7 @@ public: private: void updateEnsembleCurves( const std::vector& curves ); void updateStatisticsCurves( const std::vector& curves ); - void updateStatistics( const std::vector& sumCases ); + bool updateStatistics( const std::vector& sumCases ); caf::PdmFieldHandle* userDescriptionField() override; caf::PdmFieldHandle* objectToggleField() override; diff --git a/ApplicationLibCode/ProjectDataModel/WellLog/RimEnsembleWellLogStatistics.cpp b/ApplicationLibCode/ProjectDataModel/WellLog/RimEnsembleWellLogStatistics.cpp index f28cb04898..212a7ed5a3 100644 --- a/ApplicationLibCode/ProjectDataModel/WellLog/RimEnsembleWellLogStatistics.cpp +++ b/ApplicationLibCode/ProjectDataModel/WellLog/RimEnsembleWellLogStatistics.cpp @@ -87,6 +87,7 @@ void RimEnsembleWellLogStatistics::calculate( const std::vector const QString& wellLogChannelName ) { RiaCurveMerger curveMerger; + RiaCurveMerger tvdCurveMerger; RiaWeightedMeanCalculator dataSetSizeCalc; @@ -111,12 +112,14 @@ void RimEnsembleWellLogStatistics::calculate( const std::vector } m_logChannelUnitString = logChannelUnitString; - std::vector depths = fileData->depthValues(); - std::vector values = fileData->values( wellLogChannelName ); - if ( !depths.empty() && !values.empty() ) + std::vector depths = fileData->depthValues(); + std::vector tvdMslValues = fileData->tvdMslValues(); + std::vector values = fileData->values( wellLogChannelName ); + if ( !depths.empty() && !values.empty() && !tvdMslValues.empty() ) { dataSetSizeCalc.addValueAndWeight( depths.size(), 1.0 ); curveMerger.addCurveData( depths, values ); + tvdCurveMerger.addCurveData( depths, tvdMslValues ); } } else @@ -125,6 +128,7 @@ void RimEnsembleWellLogStatistics::calculate( const std::vector } } curveMerger.computeInterpolatedValues( true ); + tvdCurveMerger.computeInterpolatedValues( true ); clearData(); @@ -144,6 +148,7 @@ void RimEnsembleWellLogStatistics::calculate( const std::vector const std::vector& curveValues = curveMerger.interpolatedYValuesForAllXValues( curveIdx ); valuesAtDepth.push_back( curveValues[depthIdx] ); } + double p10, p50, p90, mean; RigStatisticsMath::calculateStatisticsCurves( valuesAtDepth, &p10, @@ -152,6 +157,30 @@ void RimEnsembleWellLogStatistics::calculate( const std::vector &mean, RigStatisticsMath::PercentileStyle::SWITCHED ); + // TVD is the mean TVD at a given MD + std::vector tvdsAtDepth; + tvdsAtDepth.reserve( tvdCurveMerger.curveCount() ); + + for ( size_t curveIdx = 0; curveIdx < tvdCurveMerger.curveCount(); ++curveIdx ) + { + const std::vector& curveValues = tvdCurveMerger.interpolatedYValuesForAllXValues( curveIdx ); + tvdsAtDepth.push_back( curveValues[depthIdx] ); + } + + double sumTvds = 0.0; + int numTvds = 0; + for ( auto tvd : tvdsAtDepth ) + { + if ( !std::isinf( tvd ) ) + { + sumTvds += tvd; + numTvds++; + } + } + + double meanTvd = sumTvds / numTvds; + + m_tvDepths.push_back( meanTvd ); m_measuredDepths.push_back( allDepths[depthIdx] ); m_p10Data.push_back( p10 ); m_p50Data.push_back( p50 ); @@ -217,8 +246,8 @@ void RimEnsembleWellLogStatistics::calculateByKLayer( const std::vector kIndexes = offsets->sortedIndexes(); for ( auto kIndex : kIndexes ) { - double topMean = -1.0; - double bottomMean = -2.3; + double topMean = 0.0; + double bottomMean = 0.0; // Top first { std::vector valuesAtDepth = topValues[kIndex]; @@ -229,7 +258,8 @@ void RimEnsembleWellLogStatistics::calculateByKLayer( const std::vectorgetTopDepth( kIndex ) ); + m_measuredDepths.push_back( offsets->getTopMd( kIndex ) ); + m_tvDepths.push_back( offsets->getTopTvd( kIndex ) ); m_p10Data.push_back( p10 ); m_p50Data.push_back( p50 ); m_p90Data.push_back( p90 ); @@ -248,7 +278,8 @@ void RimEnsembleWellLogStatistics::calculateByKLayer( const std::vectorgetBottomDepth( kIndex ) ); + m_measuredDepths.push_back( offsets->getBottomMd( kIndex ) ); + m_tvDepths.push_back( offsets->getBottomTvd( kIndex ) ); m_p10Data.push_back( p10 ); m_p50Data.push_back( p50 ); m_p90Data.push_back( p90 ); @@ -259,8 +290,8 @@ void RimEnsembleWellLogStatistics::calculateByKLayer( const std::vectorgetTopDepth( kIndex ) ) - .arg( offsets->getBottomDepth( kIndex ) ) + .arg( offsets->getTopMd( kIndex ) ) + .arg( offsets->getBottomMd( kIndex ) ) .arg( topMean ) .arg( bottomMean ) ); } @@ -272,11 +303,13 @@ void RimEnsembleWellLogStatistics::calculateByKLayer( const std::vector RimEnsembleWellLogStatistics::calculateIndexDepthOffset( const std::vector& wellLogFiles ) { - std::map sumTopDepths; - std::map numTopDepths; + std::map sumTopMds; + std::map sumTopTvds; + std::map numTopMds; - std::map sumBottomDepths; - std::map numBottomDepths; + std::map sumBottomMds; + std::map sumBottomTvds; + std::map numBottomMds; int minLayerK = std::numeric_limits::max(); int maxLayerK = -std::numeric_limits::max(); @@ -291,11 +324,12 @@ std::shared_ptr RigWellLogFile* fileData = wellLogFile->wellLogFileData(); std::vector depths = fileData->depthValues(); + std::vector tvdDepths = fileData->tvdMslValues(); std::vector kIndexValues = fileData->values( RiaResultNames::indexKResultName() ); std::set seenTopIndexes; std::set seenBottomIndexes; - if ( !depths.empty() && !kIndexValues.empty() ) + if ( !depths.empty() && !tvdDepths.empty() && !kIndexValues.empty() ) { // Find top indexes for ( size_t i = 0; i < kIndexValues.size(); i++ ) @@ -307,8 +341,9 @@ std::shared_ptr // This is depth of the top of the index since the file is // sorted by increasing depth. seenTopIndexes.insert( kLayer ); - sumTopDepths[kLayer] += depths[i]; - numTopDepths[kLayer] += 1; + sumTopMds[kLayer] += depths[i]; + sumTopTvds[kLayer] += tvdDepths[i]; + numTopMds[kLayer] += 1; minLayerK = std::min( minLayerK, kLayer ); maxLayerK = std::max( maxLayerK, kLayer ); } @@ -324,8 +359,9 @@ std::shared_ptr // This is depth of the bottom of the index since the file is // sorted by increasing depth. seenBottomIndexes.insert( kLayer ); - sumBottomDepths[kLayer] += depths[i]; - numBottomDepths[kLayer] += 1; + sumBottomMds[kLayer] += depths[i]; + sumBottomTvds[kLayer] += tvdDepths[i]; + numBottomMds[kLayer] += 1; } } } @@ -346,17 +382,19 @@ std::shared_ptr std::shared_ptr offset = std::make_shared(); for ( int kLayer = minLayerK; kLayer <= maxLayerK; kLayer++ ) { - if ( numTopDepths[kLayer] > 0 && numBottomDepths[kLayer] > 0 ) + if ( numTopMds[kLayer] > 0 && numBottomMds[kLayer] > 0 ) { - double topDepth = sumTopDepths[kLayer] / numTopDepths[kLayer]; - double bottomDepth = sumBottomDepths[kLayer] / numBottomDepths[kLayer]; + double topMd = sumTopMds[kLayer] / numTopMds[kLayer]; + double bottomMd = sumBottomMds[kLayer] / numBottomMds[kLayer]; + double topTvd = sumTopTvds[kLayer] / numBottomMds[kLayer]; + double bottomTvd = sumBottomTvds[kLayer] / numBottomMds[kLayer]; RiaLogging::debug( QString( "K: %1 mean depth range: %2 - %3 Samples: %4 - %5" ) .arg( kLayer ) - .arg( topDepth ) - .arg( bottomDepth ) - .arg( numTopDepths[kLayer] ) - .arg( numBottomDepths[kLayer] ) ); - offset->setIndexOffsetDepth( kLayer, topDepth, bottomDepth ); + .arg( topMd ) + .arg( bottomMd ) + .arg( numTopMds[kLayer] ) + .arg( numBottomMds[kLayer] ) ); + offset->setIndexOffsetDepth( kLayer, topMd, bottomMd, topTvd, bottomTvd ); } } diff --git a/ApplicationLibCode/ProjectDataModel/WellLog/RimEnsembleWellLogStatistics.h b/ApplicationLibCode/ProjectDataModel/WellLog/RimEnsembleWellLogStatistics.h index 68dc295183..0b647afb98 100644 --- a/ApplicationLibCode/ProjectDataModel/WellLog/RimEnsembleWellLogStatistics.h +++ b/ApplicationLibCode/ProjectDataModel/WellLog/RimEnsembleWellLogStatistics.h @@ -72,12 +72,12 @@ public: static std::shared_ptr calculateIndexDepthOffset( const std::vector& wellLogFiles ); + void clearData(); + private: void calculate( const std::vector& sumCases, const QString& wellLogChannelName ); void calculateByKLayer( const std::vector& sumCases, const QString& wellLogChannelName ); - void clearData(); - QString m_logChannelUnitString; RiaDefines::DepthUnitType m_depthUnit; std::vector m_measuredDepths; diff --git a/ApplicationLibCode/ProjectDataModel/WellLog/RimEnsembleWellLogStatisticsCurve.cpp b/ApplicationLibCode/ProjectDataModel/WellLog/RimEnsembleWellLogStatisticsCurve.cpp index ae5569e292..ca1b996876 100644 --- a/ApplicationLibCode/ProjectDataModel/WellLog/RimEnsembleWellLogStatisticsCurve.cpp +++ b/ApplicationLibCode/ProjectDataModel/WellLog/RimEnsembleWellLogStatisticsCurve.cpp @@ -108,36 +108,49 @@ void RimEnsembleWellLogStatisticsCurve::performDataExtraction( bool* isUsingPseu { values = ensembleWellLogStatistics->mean(); measuredDepthValues = ensembleWellLogStatistics->measuredDepths(); + tvDepthValues = ensembleWellLogStatistics->tvDepths(); } else if ( m_statisticsType == RimEnsembleWellLogStatistics::StatisticsType::P10 ) { values = ensembleWellLogStatistics->p10(); measuredDepthValues = ensembleWellLogStatistics->measuredDepths(); + tvDepthValues = ensembleWellLogStatistics->tvDepths(); } else if ( m_statisticsType == RimEnsembleWellLogStatistics::StatisticsType::P50 ) { values = ensembleWellLogStatistics->p50(); measuredDepthValues = ensembleWellLogStatistics->measuredDepths(); + tvDepthValues = ensembleWellLogStatistics->tvDepths(); } else if ( m_statisticsType == RimEnsembleWellLogStatistics::StatisticsType::P90 ) { values = ensembleWellLogStatistics->p90(); measuredDepthValues = ensembleWellLogStatistics->measuredDepths(); + tvDepthValues = ensembleWellLogStatistics->tvDepths(); } - bool performDataSmoothing = false; if ( !values.empty() && !measuredDepthValues.empty() && measuredDepthValues.size() == values.size() ) { if ( m_ensembleWellLogCurveSet->depthEqualization() == RimEnsembleWellLogStatistics::DepthEqualization::NONE ) + { + std::vector dummyValues( values.begin(), values.end() ); addDatapointsForBottomOfSegment( measuredDepthValues, values ); + // Pass copy of values vector to avoid adding bottom segments twice + addDatapointsForBottomOfSegment( tvDepthValues, dummyValues ); + } - this->setValuesAndDepths( values, - measuredDepthValues, - RiaDefines::DepthTypeEnum::MEASURED_DEPTH, - rkbDiff, - depthUnit, - !performDataSmoothing, - xUnits ); + std::map> validDepths; + if ( values.size() == measuredDepthValues.size() ) + { + validDepths.insert( std::make_pair( RiaDefines::DepthTypeEnum::MEASURED_DEPTH, measuredDepthValues ) ); + } + if ( values.size() == tvDepthValues.size() ) + { + validDepths.insert( std::make_pair( RiaDefines::DepthTypeEnum::TRUE_VERTICAL_DEPTH, tvDepthValues ) ); + validDepths.insert( std::make_pair( RiaDefines::DepthTypeEnum::TRUE_VERTICAL_DEPTH_RKB, tvDepthValues ) ); + } + + this->setValuesAndDepths( values, validDepths, rkbDiff, depthUnit, false ); } } } diff --git a/ApplicationLibCode/ProjectDataModel/WellLog/RimWellLogFileCurve.cpp b/ApplicationLibCode/ProjectDataModel/WellLog/RimWellLogFileCurve.cpp index e066c596cf..fecf2f9492 100644 --- a/ApplicationLibCode/ProjectDataModel/WellLog/RimWellLogFileCurve.cpp +++ b/ApplicationLibCode/ProjectDataModel/WellLog/RimWellLogFileCurve.cpp @@ -242,11 +242,11 @@ std::pair, std::vector> if ( firstIndex != kIndexValues.size() && lastIndex != kIndexValues.size() ) { // Add top - measuredDepthValuesAdjusted.push_back( m_indexDepthOffsets->getTopDepth( kLayer ) ); + measuredDepthValuesAdjusted.push_back( m_indexDepthOffsets->getTopMd( kLayer ) ); valuesAdjusted.push_back( values[firstIndex] ); // Add bottom of layer - measuredDepthValuesAdjusted.push_back( m_indexDepthOffsets->getBottomDepth( kLayer ) ); + measuredDepthValuesAdjusted.push_back( m_indexDepthOffsets->getBottomMd( kLayer ) ); valuesAdjusted.push_back( values[lastIndex] ); } } diff --git a/ApplicationLibCode/ReservoirDataModel/RigWellLogIndexDepthOffset.cpp b/ApplicationLibCode/ReservoirDataModel/RigWellLogIndexDepthOffset.cpp index d3f45dc623..3a37b52043 100644 --- a/ApplicationLibCode/ReservoirDataModel/RigWellLogIndexDepthOffset.cpp +++ b/ApplicationLibCode/ReservoirDataModel/RigWellLogIndexDepthOffset.cpp @@ -23,18 +23,19 @@ //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- -void RigWellLogIndexDepthOffset::setIndexOffsetDepth( int kIndex, double topDepth, double bottomDepth ) +void RigWellLogIndexDepthOffset::setIndexOffsetDepth( int kIndex, double topMd, double bottomMd, double topTvd, double bottomTvd ) { - m_depthOffsets[kIndex] = std::pair( topDepth, bottomDepth ); + m_mdOffsets[kIndex] = std::pair( topMd, bottomMd ); + m_tvdOffsets[kIndex] = std::pair( topTvd, bottomTvd ); } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- -double RigWellLogIndexDepthOffset::getTopDepth( int kIndex ) const +double RigWellLogIndexDepthOffset::getTopMd( int kIndex ) const { - auto hit = m_depthOffsets.find( kIndex ); - if ( hit != m_depthOffsets.end() ) + auto hit = m_mdOffsets.find( kIndex ); + if ( hit != m_mdOffsets.end() ) { return hit->second.first; } @@ -45,10 +46,38 @@ double RigWellLogIndexDepthOffset::getTopDepth( int kIndex ) const //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- -double RigWellLogIndexDepthOffset::getBottomDepth( int kIndex ) const +double RigWellLogIndexDepthOffset::getBottomMd( int kIndex ) const { - auto hit = m_depthOffsets.find( kIndex ); - if ( hit != m_depthOffsets.end() ) + auto hit = m_mdOffsets.find( kIndex ); + if ( hit != m_mdOffsets.end() ) + { + return hit->second.second; + } + + return std::numeric_limits::infinity(); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +double RigWellLogIndexDepthOffset::getTopTvd( int kIndex ) const +{ + auto hit = m_tvdOffsets.find( kIndex ); + if ( hit != m_tvdOffsets.end() ) + { + return hit->second.first; + } + + return std::numeric_limits::infinity(); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +double RigWellLogIndexDepthOffset::getBottomTvd( int kIndex ) const +{ + auto hit = m_tvdOffsets.find( kIndex ); + if ( hit != m_tvdOffsets.end() ) { return hit->second.second; } @@ -62,7 +91,7 @@ double RigWellLogIndexDepthOffset::getBottomDepth( int kIndex ) const std::vector RigWellLogIndexDepthOffset::sortedIndexes() const { std::vector indexes; - for ( auto m : m_depthOffsets ) + for ( auto m : m_mdOffsets ) { indexes.push_back( m.first ); } diff --git a/ApplicationLibCode/ReservoirDataModel/RigWellLogIndexDepthOffset.h b/ApplicationLibCode/ReservoirDataModel/RigWellLogIndexDepthOffset.h index cd5bf2c8d1..3bf5cb1852 100644 --- a/ApplicationLibCode/ReservoirDataModel/RigWellLogIndexDepthOffset.h +++ b/ApplicationLibCode/ReservoirDataModel/RigWellLogIndexDepthOffset.h @@ -28,11 +28,14 @@ public: RigWellLogIndexDepthOffset() = default; ~RigWellLogIndexDepthOffset() override = default; - void setIndexOffsetDepth( int kIndex, double topDepth, double bottomDepth ); - double getTopDepth( int kIndex ) const; - double getBottomDepth( int kIndex ) const; + void setIndexOffsetDepth( int kIndex, double topMd, double bottomMd, double topTvd, double bottomTvd ); + double getTopMd( int kIndex ) const; + double getBottomMd( int kIndex ) const; + double getTopTvd( int kIndex ) const; + double getBottomTvd( int kIndex ) const; std::vector sortedIndexes() const; private: - std::map> m_depthOffsets; + std::map> m_mdOffsets; + std::map> m_tvdOffsets; };