From 89f90ee9a9fa4dea96881c330a570807b00cd6e6 Mon Sep 17 00:00:00 2001 From: Kristian Bendiksen Date: Wed, 31 Jan 2024 21:20:08 +0100 Subject: [PATCH] Fault reactivation: use values from GeoMech grid for stress. No longer snapping to the fake well path along the border between the parts. --- ...ltReactivationDataAccessorPorePressure.cpp | 10 +++ ...tReactivationDataAccessorStressGeoMech.cpp | 75 ++++++++++++------- ...ultReactivationDataAccessorStressGeoMech.h | 15 ++-- ...ctivationDataAccessorWellLogExtraction.cpp | 12 ++- 4 files changed, 76 insertions(+), 36 deletions(-) diff --git a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorPorePressure.cpp b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorPorePressure.cpp index 487c937ea0..c569224881 100644 --- a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorPorePressure.cpp +++ b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorPorePressure.cpp @@ -114,6 +114,16 @@ double RimFaultReactivationDataAccessorPorePressure::valueAtPosition( const cvf: auto [value, pos] = RimFaultReactivationDataAccessorWellLogExtraction::calculatePorBar( intersections, values, position, m_defaultPorePressureGradient ); + if ( pos.isUndefined() ) + { + auto cellIdx = m_mainGrid->findReservoirCellIndexFromPoint( position ); + if ( cellIdx != cvf::UNDEFINED_SIZE_T ) + { + double valueFromEclipse = m_resultAccessor->cellScalar( cellIdx ); + if ( !std::isinf( valueFromEclipse ) ) return RiaEclipseUnitTools::barToPascal( valueFromEclipse ); + } + } + return RiaEclipseUnitTools::barToPascal( value ); } diff --git a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStressGeoMech.cpp b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStressGeoMech.cpp index a1c000f5a6..a19cfa4c52 100644 --- a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStressGeoMech.cpp +++ b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStressGeoMech.cpp @@ -88,6 +88,7 @@ void RimFaultReactivationDataAccessorStressGeoMech::updateResultAccessor() m_s33Frames = loadFrameLambda( femParts, getResultAddress( "ST", "S33" ), timeStepIndex ); m_s11Frames = loadFrameLambda( femParts, getResultAddress( "ST", "S11" ), timeStepIndex ); m_s22Frames = loadFrameLambda( femParts, getResultAddress( "ST", "S22" ), timeStepIndex ); + m_porBarFrames = loadFrameLambda( femParts, RigFemAddressDefines::nodalPorBarAddress(), timeStepIndex ); auto [faultTopPosition, faultBottomPosition] = m_model->faultTopBottom(); auto faultNormal = m_model->modelNormal() ^ cvf::Vec3d::Z_AXIS; @@ -96,31 +97,25 @@ void RimFaultReactivationDataAccessorStressGeoMech::updateResultAccessor() double distanceFromFault = 1.0; auto [topDepth, bottomDepth] = m_model->depthTopBottom(); - RigFemPartCollection* geoMechPartCollection = m_geoMechCaseData->femParts(); - std::string errorName = "fault reactivation data access"; - + for ( auto gridPart : m_model->allGridParts() ) { + double sign = m_model->normalPointsAt() == gridPart ? -1.0 : 1.0; std::vector wellPoints = RimFaultReactivationDataAccessorWellLogExtraction::generateWellPoints( faultTopPosition, faultBottomPosition, m_seabedDepth, bottomDepth, - faultNormal * distanceFromFault ); - m_faceAWellPath = new RigWellPath( wellPoints, RimFaultReactivationDataAccessorWellLogExtraction::generateMds( wellPoints ) ); - m_partIndexA = geoMechPartCollection->getPartIndexFromPoint( wellPoints[1] ); - m_extractorA = new RigGeoMechWellLogExtractor( m_geoMechCaseData, partIndex, m_faceAWellPath.p(), errorName ); - } + sign * faultNormal * distanceFromFault ); - { - std::vector wellPoints = - RimFaultReactivationDataAccessorWellLogExtraction::generateWellPoints( faultTopPosition, - faultBottomPosition, - m_seabedDepth, - bottomDepth, - -faultNormal * distanceFromFault ); - m_faceBWellPath = new RigWellPath( wellPoints, RimFaultReactivationDataAccessorWellLogExtraction::generateMds( wellPoints ) ); - m_partIndexB = geoMechPartCollection->getPartIndexFromPoint( wellPoints[1] ); - m_extractorB = new RigGeoMechWellLogExtractor( m_geoMechCaseData, partIndex, m_faceBWellPath.p(), errorName ); + cvf::ref wellPath = + new RigWellPath( wellPoints, RimFaultReactivationDataAccessorWellLogExtraction::generateMds( wellPoints ) ); + m_wellPaths[gridPart] = wellPath; + + std::string errorName = "fault reactivation data access"; + cvf::ref extractor = + new RigGeoMechWellLogExtractor( m_geoMechCaseData, partIndex, wellPath.p(), errorName ); + + m_extractors[gridPart] = extractor; } } @@ -138,7 +133,7 @@ RigFemResultAddress RimFaultReactivationDataAccessorStressGeoMech::getResultAddr //-------------------------------------------------------------------------------------------------- bool RimFaultReactivationDataAccessorStressGeoMech::isDataAvailable() const { - return m_s11Frames && m_s22Frames && m_s33Frames && m_femPart; + return m_s11Frames && m_s22Frames && m_s33Frames && m_porBarFrames && m_femPart; } //-------------------------------------------------------------------------------------------------- @@ -183,7 +178,7 @@ std::pair RimFaultReactivationDataAccessorStressGeoMech::cal { int timeStepIndex = 0; int frameIndex = m_s33Frames->frameCount( timeStepIndex ) - 1; - return calculatePorBar( position, m_gradient, timeStepIndex, frameIndex ); + return calculatePorBar( position, m_gradient, gridPart, timeStepIndex, frameIndex ); } //-------------------------------------------------------------------------------------------------- @@ -215,14 +210,15 @@ double RimFaultReactivationDataAccessorStressGeoMech::interpolatedResultValue( R //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- -std::pair RimFaultReactivationDataAccessorStressGeoMech::calculatePorBar( const cvf::Vec3d& position, - double gradient, - int timeStepIndex, - int frameIndex ) const +std::pair RimFaultReactivationDataAccessorStressGeoMech::calculatePorBar( const cvf::Vec3d& position, + double gradient, + RimFaultReactivation::GridPart gridPart, + int timeStepIndex, + int frameIndex ) const { - RigFemPartCollection* partCollection = m_geoMechCaseData->femParts(); - cvf::ref extractor = m_partIndexA == partCollection->getPartIndexFromPoint( position ) ? m_extractorA - : m_extractorB; + CAF_ASSERT( m_extractors.find( gridPart ) != m_extractors.end() ); + auto extractor = m_extractors.find( gridPart )->second; + if ( !extractor->valid() ) { RiaLogging::error( "Invalid extractor when extracting PorBar" ); @@ -233,5 +229,28 @@ std::pair RimFaultReactivationDataAccessorStressGeoMech::cal std::vector values; extractor->curveData( resAddr, timeStepIndex, frameIndex, &values ); - return RimFaultReactivationDataAccessorWellLogExtraction::calculatePorBar( extractor->intersections(), values, position, gradient ); + auto [value, extractionPos] = + RimFaultReactivationDataAccessorWellLogExtraction::calculatePorBar( extractor->intersections(), values, position, gradient ); + + if ( extractionPos.isUndefined() ) + { + // If extraction position is not defined the position is not close to the border between the two parts. + // This means it should be safe to use POR-BAR from the model. + const std::vector& frameData = m_porBarFrames->frameData( timeStepIndex, frameIndex ); + + // Use data from geo mech grid if defined (only position is reservoir). + RimWellIADataAccess iaDataAccess( m_geoMechCase ); + double gridValue = iaDataAccess.interpolatedResultValue( m_femPart, frameData, RIG_NODAL, position ); + if ( !std::isinf( gridValue ) ) + { + return { gridValue, position }; + } + + // Use calculated value when POR-BAR is inf (outside of reservoir). + return { value, position }; + } + else + { + return { value, extractionPos }; + } } diff --git a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStressGeoMech.h b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStressGeoMech.h index 5405bfac31..87d581bb55 100644 --- a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStressGeoMech.h +++ b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStressGeoMech.h @@ -77,7 +77,11 @@ private: const cvf::Vec3d& position, const std::vector& scalarResults ) const; - std::pair calculatePorBar( const cvf::Vec3d& position, double gradient, int timeStepIndex, int frameIndex ) const; + std::pair calculatePorBar( const cvf::Vec3d& position, + double gradient, + RimFaultReactivation::GridPart gridPart, + int timeStepIndex, + int frameIndex ) const; RigFemScalarResultFrames* dataFrames( StressType stressType ) const; @@ -86,12 +90,9 @@ private: RigFemScalarResultFrames* m_s11Frames; RigFemScalarResultFrames* m_s22Frames; RigFemScalarResultFrames* m_s33Frames; + RigFemScalarResultFrames* m_porBarFrames; const RigFemPart* m_femPart; - cvf::ref m_faceAWellPath; - cvf::ref m_faceBWellPath; - int m_partIndexA; - int m_partIndexB; - cvf::ref m_extractorA; - cvf::ref m_extractorB; + std::map> m_wellPaths; + std::map> m_extractors; }; diff --git a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorWellLogExtraction.cpp b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorWellLogExtraction.cpp index 39a752e101..c5afb0513f 100644 --- a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorWellLogExtraction.cpp +++ b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorWellLogExtraction.cpp @@ -70,7 +70,17 @@ std::pair RimFaultReactivationDataAccessorWellLogExtraction: { // Fill in missing values fillInMissingValuesWithGradient( intersections, values, gradient ); - return findValueAndPosition( intersections, values, position ); + auto [value, extractionPosition] = findValueAndPosition( intersections, values, position ); + + double minDistance = computeMinimumDistance( position, intersections ); + if ( minDistance < 1.0 ) + { + return { value, extractionPosition }; + } + else + { + return { value, cvf::Vec3d::UNDEFINED }; + } } //--------------------------------------------------------------------------------------------------