diff --git a/ApplicationLibCode/ReservoirDataModel/RigGeoMechWellLogExtractor.cpp b/ApplicationLibCode/ReservoirDataModel/RigGeoMechWellLogExtractor.cpp index 86ded27774b..7226c306208 100644 --- a/ApplicationLibCode/ReservoirDataModel/RigGeoMechWellLogExtractor.cpp +++ b/ApplicationLibCode/ReservoirDataModel/RigGeoMechWellLogExtractor.cpp @@ -152,18 +152,12 @@ QString RigGeoMechWellLogExtractor::curveData( const RigFemResultAddress& resAdd wellBoreFGShale( RigWbsParameter::FG_Shale(), timeStepIndex, frameIndex, values ); values->front() = wbsCurveValuesAtMsl(); } - else if ( resAddr.fieldName == RiaResultNames::wbsSFGResult().toStdString() ) + else if ( resAddr.fieldName == RiaResultNames::wbsSFGResult().toStdString() || + resAddr.fieldName == RiaResultNames::wbsFGMkMinResult().toStdString() || + resAddr.fieldName == RiaResultNames::wbsFGMkExpResult().toStdString() ) { wellBoreWallCurveData( resAddr, timeStepIndex, frameIndex, values ); } - else if ( resAddr.fieldName == RiaResultNames::wbsFGMkMinResult().toStdString() ) - { - wellBoreFG_MatthewsKelly( RigWbsParameter::FG_MkMin(), timeStepIndex, frameIndex, values ); - } - else if ( resAddr.fieldName == RiaResultNames::wbsFGMkExpResult().toStdString() ) - { - wellBoreFG_MatthewsKelly( RigWbsParameter::FG_MkExp(), timeStepIndex, frameIndex, values ); - } else if ( resAddr.fieldName == RiaResultNames::wbsPPResult().toStdString() || resAddr.fieldName == RiaResultNames::wbsOBGResult().toStdString() || resAddr.fieldName == RiaResultNames::wbsSHResult().toStdString() ) @@ -628,34 +622,18 @@ void RigGeoMechWellLogExtractor::wellBoreWallCurveData( const RigFemResultAddres resAddr.fieldName == RiaResultNames::wbsFGMkMinResult().toStdString() || resAddr.fieldName == RiaResultNames::wbsFGMkExpResult().toStdString() ); - // The result addresses needed - RigFemResultAddress stressResAddr( RIG_ELEMENT_NODAL, "ST", "" ); - RigFemResultAddress porBarResAddr = RigFemAddressDefines::elementNodalPorBarAddress(); - RigFemPartResultsCollection* resultCollection = m_caseData->femPartResults(); auto mapFGResultToPP = []( const QString& fgResultName ) { - if ( fgResultName == RiaResultNames::wbsFGMkMinResult() ) return RigWbsParameter::PP_Min(); - if ( fgResultName == RiaResultNames::wbsFGMkExpResult() ) return RigWbsParameter::PP_Exp(); + // if ( fgResultName == RiaResultNames::wbsFGMkMinResult() ) return RigWbsParameter::PP_Min(); + // if ( fgResultName == RiaResultNames::wbsFGMkExpResult() ) return RigWbsParameter::PP_Exp(); return RigWbsParameter::PP_Reservoir(); }; - RigWbsParameter ppParameter = mapFGResultToPP( QString::fromStdString( resAddr.fieldName ) ); - - // Load results - std::vector vertexStressesFloat = resultCollection->tensors( stressResAddr, m_partId, timeStepIndex, frameIndex ); - if ( vertexStressesFloat.empty() ) return; - - std::vector vertexStresses; - vertexStresses.reserve( vertexStressesFloat.size() ); - for ( const caf::Ten3f& floatTensor : vertexStressesFloat ) - { - vertexStresses.push_back( caf::Ten3d( floatTensor ) ); - } - - std::vector interpolatedInterfaceStressBar = - interpolateInterfaceValues( stressResAddr, timeStepIndex, frameIndex, vertexStresses ); + RigWbsParameter ppParameter = mapFGResultToPP( QString::fromStdString( resAddr.fieldName ) ); + bool useGridStress = !( resAddr.fieldName == RiaResultNames::wbsFGMkMinResult().toStdString() || + resAddr.fieldName == RiaResultNames::wbsFGMkExpResult().toStdString() ); values->resize( intersections().size(), std::numeric_limits::infinity() ); @@ -669,15 +647,79 @@ void RigGeoMechWellLogExtractor::wellBoreWallCurveData( const RigFemResultAddres std::vector ucsAllSegments( intersections().size(), std::numeric_limits::infinity() ); calculateWbsParameterForAllSegments( RigWbsParameter::UCS(), timeStepIndex, frameIndex, &ucsAllSegments, false ); + std::vector> segmentStresses( intersections().size(), { caf::Ten3d::invalid(), false } ); + + if ( useGridStress ) + { + // The result addresses needed + RigFemResultAddress stressResAddr( RIG_ELEMENT_NODAL, "ST", "" ); + + // Load results + std::vector vertexStressesFloat = resultCollection->tensors( stressResAddr, m_partId, timeStepIndex, frameIndex ); + if ( vertexStressesFloat.empty() ) return; + + std::vector vertexStresses; + vertexStresses.reserve( vertexStressesFloat.size() ); + for ( const caf::Ten3f& floatTensor : vertexStressesFloat ) + { + vertexStresses.push_back( caf::Ten3d( floatTensor ) ); + } + + std::vector interpolatedInterfaceStressBar = + interpolateInterfaceValues( stressResAddr, timeStepIndex, frameIndex, vertexStresses ); + +#pragma omp parallel for + for ( int64_t intersectionIdx = 0; intersectionIdx < static_cast( intersections().size() ); ++intersectionIdx ) + { + caf::Ten3d segmentStress; + bool validSegmentStress = + averageIntersectionValuesToSegmentValue( intersectionIdx, interpolatedInterfaceStressBar, caf::Ten3d::invalid(), &segmentStress ); + segmentStresses[intersectionIdx] = { segmentStress, validSegmentStress }; + } + } + else + { + std::vector obgAllSegments( intersections().size(), std::numeric_limits::infinity() ); + calculateWbsParameterForAllSegments( RigWbsParameter::OBG(), timeStepIndex, frameIndex, &obgAllSegments, false ); + + auto mapFGResultToSH = []( const QString& fgResultName ) + { + if ( fgResultName == RiaResultNames::wbsFGMkMinResult() ) + return RiaResultNames::wbsSHMkMinResult(); + else + return RiaResultNames::wbsSHMkExpResult(); + }; + + std::vector SH; + QString SHMkResultName = mapFGResultToSH( QString::fromStdString( resAddr.fieldName ) ); + RigFemResultAddress SHMkAddr( RIG_WELLPATH_DERIVED, SHMkResultName.toStdString(), "" ); + + curveData( SHMkAddr, timeStepIndex, frameIndex, &SH ); + + CVF_ASSERT( SH.size() == intersections().size() ); + +#pragma omp parallel for + for ( int64_t intersectionIdx = 0; intersectionIdx < static_cast( intersections().size() ); ++intersectionIdx ) + { + double horizontalStress = obgAllSegments[intersectionIdx]; + double verticalStress = SH[intersectionIdx]; + double shear = 0.0; + caf::Ten3d segmentStress( horizontalStress, horizontalStress, verticalStress, shear, shear, shear ); + // Only for pp defined?? + bool validSegmentStress = true; + segmentStresses[intersectionIdx] = { segmentStress, validSegmentStress }; + } + } + + CAF_ASSERT( segmentStresses.size() == intersections().size() ); + #pragma omp parallel for for ( int64_t intersectionIdx = 0; intersectionIdx < static_cast( intersections().size() ); ++intersectionIdx ) { - // FG is for sands, SFG for shale. Sands has valid PP, shale does not. bool isFGregion = ppSources[intersectionIdx] == RigWbsParameter::GRID; - if ( resAddr.fieldName == RiaResultNames::wbsFGMkMinResult().toStdString() || - resAddr.fieldName == RiaResultNames::wbsFGMkExpResult().toStdString() ) + // TODO: hack + if ( !useGridStress ) { - // Assume only FG for entire well log for FG_MK_MIN/EXP. isFGregion = true; } @@ -692,15 +734,11 @@ void RigGeoMechWellLogExtractor::wellBoreWallCurveData( const RigFemResultAddres double poissonRatio = poissonAllSegments[intersectionIdx]; double ucsBar = ucsAllSegments[intersectionIdx]; - caf::Ten3d segmentStress; - bool validSegmentStress = - averageIntersectionValuesToSegmentValue( intersectionIdx, interpolatedInterfaceStressBar, caf::Ten3d::invalid(), &segmentStress ); + auto [segmentStress, validSegmentStress] = segmentStresses[intersectionIdx]; + cvf::Vec3d wellPathTangent = calculateWellPathTangent( intersectionIdx, TangentConstantWithinCell ); + caf::Ten3d wellPathStress = transformTensorToWellPathOrientation( wellPathTangent, segmentStress ); - cvf::Vec3d wellPathTangent = calculateWellPathTangent( intersectionIdx, TangentConstantWithinCell ); - caf::Ten3d wellPathStressFloat = transformTensorToWellPathOrientation( wellPathTangent, segmentStress ); - caf::Ten3d wellPathStressDouble( wellPathStressFloat ); - - RigGeoMechBoreHoleStressCalculator sigmaCalculator( wellPathStressDouble, porePressureBar, poissonRatio, ucsBar, 32 ); + RigGeoMechBoreHoleStressCalculator sigmaCalculator( wellPathStress, porePressureBar, poissonRatio, ucsBar, 32 ); double resultValue = std::numeric_limits::infinity(); if ( resAddr.fieldName == RiaResultNames::wbsFGResult().toStdString() || resAddr.fieldName == RiaResultNames::wbsFGMkMinResult().toStdString() || @@ -812,66 +850,6 @@ void RigGeoMechWellLogExtractor::wellBoreFGDerivedFromK0FG( const QString& } } -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -void RigGeoMechWellLogExtractor::wellBoreFG_MatthewsKelly( const RigWbsParameter& parameter, - int timeStepIndex, - int frameIndex, - std::vector* values ) -{ - values->resize( intersections().size(), std::numeric_limits::infinity() ); - - // Use FG_Shale source to avoid creating more options. - WbsParameterSource source = m_parameterSources.at( RigWbsParameter::FG_Shale() ); - if ( source == RigWbsParameter::DERIVED_FROM_K0FG ) - { - auto mapParameterToPPResult = []( const RigWbsParameter& parameter ) - { - if ( parameter.name() == RiaResultNames::wbsFGMkMinResult() ) return RiaResultNames::wbsPPMinResult(); - if ( parameter.name() == RiaResultNames::wbsFGMkExpResult() ) return RiaResultNames::wbsPPExpResult(); - return RiaResultNames::wbsPPResult(); - }; - - QString ppResultName = mapParameterToPPResult( parameter ); - bool onlyForPPReservoir = true; - wellBoreFGDerivedFromK0FG( ppResultName, timeStepIndex, frameIndex, values, onlyForPPReservoir ); - } - else - { - auto mapParameterToSHMkResult = []( const RigWbsParameter& parameter ) - { - if ( parameter.name() == RiaResultNames::wbsFGMkMinResult() ) return RiaResultNames::wbsSHMkMinResult(); - if ( parameter.name() == RiaResultNames::wbsFGMkExpResult() ) return RiaResultNames::wbsSHMkExpResult(); - return RiaResultNames::wbsSHMkResult(); - }; - - std::vector SH; - QString SHMkResultName = mapParameterToSHMkResult( parameter ); - RigFemResultAddress SHMkAddr( RIG_WELLPATH_DERIVED, SHMkResultName.toStdString(), "" ); - - curveData( SHMkAddr, timeStepIndex, frameIndex, &SH ); - - CVF_ASSERT( SH.size() == intersections().size() ); - - std::vector PP( intersections().size(), std::numeric_limits::infinity() ); - std::vector ppSources = - calculateWbsParameterForAllSegments( RigWbsParameter::PP_Reservoir(), timeStepIndex, frameIndex, &PP, false ); - - double multiplier = m_userDefinedValues.at( parameter ); - CVF_ASSERT( multiplier != std::numeric_limits::infinity() ); -#pragma omp parallel for - for ( int64_t intersectionIdx = 0; intersectionIdx < static_cast( intersections().size() ); ++intersectionIdx ) - { - if ( !isValid( ( *values )[intersectionIdx] ) && ppSources[intersectionIdx] == RigWbsParameter::GRID && - isValid( PP[intersectionIdx] ) && isValid( SH[intersectionIdx] ) ) - { - ( *values )[intersectionIdx] = SH[intersectionIdx] * multiplier; - } - } - } -} - //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- diff --git a/ApplicationLibCode/ReservoirDataModel/RigGeoMechWellLogExtractor.h b/ApplicationLibCode/ReservoirDataModel/RigGeoMechWellLogExtractor.h index e846bf440d8..5178dca748d 100644 --- a/ApplicationLibCode/ReservoirDataModel/RigGeoMechWellLogExtractor.h +++ b/ApplicationLibCode/ReservoirDataModel/RigGeoMechWellLogExtractor.h @@ -133,8 +133,6 @@ class RigGeoMechWellLogExtractor : public RigWellLogExtractor void wellBoreFGDerivedFromK0FG( const QString& ppResult, int timeStepIndex, int frameIndex, std::vector* values, bool onlyForPPReservoir ); - void wellBoreFG_MatthewsKelly( const RigWbsParameter& parameter, int timeStepIndex, int frameIndex, std::vector* values ); - template T interpolateGridResultValue( RigFemResultPosEnum resultPosType, const std::vector& gridResultValues, int64_t intersectionIdx ) const; size_t gridResultIndexFace( size_t elementIdx, cvf::StructGridInterface::FaceType cellFace, int faceLocalNodeIdx ) const;