From f16c60f394914989918d40d8d24de1f6c68c9ad2 Mon Sep 17 00:00:00 2001 From: Kristian Bendiksen Date: Fri, 8 Dec 2023 16:00:20 +0100 Subject: [PATCH] WIP: try to extract porbar --- .../Faults/RimFaultReactivationDataAccess.cpp | 2 +- .../Faults/RimFaultReactivationDataAccessor.h | 5 +- ...imFaultReactivationDataAccessorGeoMech.cpp | 4 +- .../RimFaultReactivationDataAccessorGeoMech.h | 5 +- ...ltReactivationDataAccessorPorePressure.cpp | 3 +- ...aultReactivationDataAccessorPorePressure.h | 5 +- ...RimFaultReactivationDataAccessorStress.cpp | 117 +++++++++++++++++- .../RimFaultReactivationDataAccessorStress.h | 23 +++- ...ultReactivationDataAccessorTemperature.cpp | 3 +- ...FaultReactivationDataAccessorTemperature.h | 5 +- ...FaultReactivationDataAccessorVoidRatio.cpp | 3 +- ...imFaultReactivationDataAccessorVoidRatio.h | 5 +- 12 files changed, 154 insertions(+), 26 deletions(-) diff --git a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccess.cpp b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccess.cpp index 725eab421fa..f1c73fae51f 100644 --- a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccess.cpp +++ b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccess.cpp @@ -153,7 +153,7 @@ std::vector RimFaultReactivationDataAccess::extractModelData( const RigF double bottomDepth = computeAverageDepth( corners, { 4, 5, 6, 7 } ); cvf::Vec3d position = RigCaseToCaseCellMapperTools::calculateCellCenter( corners.data() ); - double value = accessor->valueAtPosition( position, model, gridPart, topDepth, bottomDepth ); + double value = accessor->valueAtPosition( position, model, gridPart, topDepth, bottomDepth, elementIndex ); values.push_back( value ); } } diff --git a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessor.h b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessor.h index 3d3d3504fc0..ddc49bdd433 100644 --- a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessor.h +++ b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessor.h @@ -45,8 +45,9 @@ class RimFaultReactivationDataAccessor virtual double valueAtPosition( const cvf::Vec3d& position, const RigFaultReactivationModel& model, RimFaultReactivation::GridPart gridPart, - double topDepth = std::numeric_limits::infinity(), - double bottomDepth = std::numeric_limits::infinity() ) const = 0; + double topDepth = std::numeric_limits::infinity(), + double bottomDepth = std::numeric_limits::infinity(), + size_t elementIndex = std::numeric_limits::max() ) const = 0; protected: virtual void updateResultAccessor() = 0; diff --git a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorGeoMech.cpp b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorGeoMech.cpp index 7da47f84c5b..4c81dfd05ba 100644 --- a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorGeoMech.cpp +++ b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorGeoMech.cpp @@ -101,8 +101,8 @@ double RimFaultReactivationDataAccessorGeoMech::valueAtPosition( const cvf::Vec3 const RigFaultReactivationModel& model, RimFaultReactivation::GridPart gridPart, double topDepth, - double bottomDepth ) const - + double bottomDepth, + size_t elementIndex ) const { if ( !m_resultFrames ) return std::numeric_limits::infinity(); diff --git a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorGeoMech.h b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorGeoMech.h index e89241d1758..5e7f9bdb068 100644 --- a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorGeoMech.h +++ b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorGeoMech.h @@ -43,8 +43,9 @@ class RimFaultReactivationDataAccessorGeoMech : public RimFaultReactivationDataA double valueAtPosition( const cvf::Vec3d& position, const RigFaultReactivationModel& model, RimFaultReactivation::GridPart gridPart, - double topDepth = std::numeric_limits::infinity(), - double bottomDepth = std::numeric_limits::infinity() ) const override; + double topDepth = std::numeric_limits::infinity(), + double bottomDepth = std::numeric_limits::infinity(), + size_t elementIndex = std::numeric_limits::max() ) const override; private: void updateResultAccessor() override; diff --git a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorPorePressure.cpp b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorPorePressure.cpp index 9d32921eb21..e9a386b1e90 100644 --- a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorPorePressure.cpp +++ b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorPorePressure.cpp @@ -90,7 +90,8 @@ double RimFaultReactivationDataAccessorPorePressure::valueAtPosition( const cvf: const RigFaultReactivationModel& model, RimFaultReactivation::GridPart gridPart, double topDepth, - double bottomDepth ) const + double bottomDepth, + size_t elementIndex ) const { if ( ( m_mainGrid != nullptr ) && m_resultAccessor.notNull() ) { diff --git a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorPorePressure.h b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorPorePressure.h index 6ae024a7a01..647537a7cb4 100644 --- a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorPorePressure.h +++ b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorPorePressure.h @@ -44,8 +44,9 @@ class RimFaultReactivationDataAccessorPorePressure : public RimFaultReactivation double valueAtPosition( const cvf::Vec3d& position, const RigFaultReactivationModel& model, RimFaultReactivation::GridPart gridPart, - double topDepth = std::numeric_limits::infinity(), - double bottomDepth = std::numeric_limits::infinity() ) const override; + double topDepth = std::numeric_limits::infinity(), + double bottomDepth = std::numeric_limits::infinity(), + size_t elementIndex = std::numeric_limits::max() ) const override; private: void updateResultAccessor() override; diff --git a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStress.cpp b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStress.cpp index e14c1d7bee7..5bb4b798151 100644 --- a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStress.cpp +++ b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStress.cpp @@ -20,15 +20,18 @@ #include "RiaEclipseUnitTools.h" +#include "RigFaultReactivationModel.h" #include "RigFemAddressDefines.h" #include "RigFemPartCollection.h" #include "RigFemPartResultsCollection.h" #include "RigFemResultAddress.h" #include "RigFemScalarResultFrames.h" #include "RigGeoMechCaseData.h" +#include "RigGriddedPart3d.h" #include "RigResultAccessorFactory.h" #include "RimFaultReactivationEnums.h" +#include "RimFracture.h" #include "RimGeoMechCase.h" #include "RimWellIADataAccess.h" @@ -101,7 +104,8 @@ double RimFaultReactivationDataAccessorStress::valueAtPosition( const cvf::Vec3d const RigFaultReactivationModel& model, RimFaultReactivation::GridPart gridPart, double topDepth, - double bottomDepth ) const + double bottomDepth, + size_t elementIndex ) const { if ( !m_porFrames || !m_s11Frames || !m_s22Frames || !m_s33Frames || !m_femPart ) return std::numeric_limits::infinity(); @@ -123,16 +127,18 @@ double RimFaultReactivationDataAccessorStress::valueAtPosition( const cvf::Vec3d const std::vector& s33Data = m_s33Frames->frameData( timeStepIndex, frameIndex ); const std::vector& porData = m_porFrames->frameData( timeStepIndex, frameIndex ); + const RigGriddedPart3d* grid = model.grid( gridPart ); + if ( m_property == RimFaultReactivation::Property::StressTop ) { double s33 = interpolatedResultValue( iaDataAccess, m_femPart, topPosition, s33Data ); - double porBar = interpolatedResultValue( iaDataAccess, m_femPart, topPosition, porData ); + double porBar = getPorBar( iaDataAccess, m_femPart, topPosition, porData, *grid, elementIndex ); return RiaEclipseUnitTools::barToPascal( s33 - porBar ); } else if ( m_property == RimFaultReactivation::Property::StressBottom ) { double s33 = interpolatedResultValue( iaDataAccess, m_femPart, bottomPosition, s33Data ); - double porBar = interpolatedResultValue( iaDataAccess, m_femPart, bottomPosition, porData ); + double porBar = getPorBar( iaDataAccess, m_femPart, bottomPosition, porData, *grid, elementIndex ); return RiaEclipseUnitTools::barToPascal( s33 - porBar ); } else if ( m_property == RimFaultReactivation::Property::DepthTop ) @@ -147,14 +153,14 @@ double RimFaultReactivationDataAccessorStress::valueAtPosition( const cvf::Vec3d { double s11 = interpolatedResultValue( iaDataAccess, m_femPart, position, s11Data ); double s33 = interpolatedResultValue( iaDataAccess, m_femPart, position, s33Data ); - double porBar = interpolatedResultValue( iaDataAccess, m_femPart, position, porData ); + double porBar = getPorBar( iaDataAccess, m_femPart, position, porData, *grid, elementIndex ); return ( s11 - porBar ) / ( s33 - porBar ); } else if ( m_property == RimFaultReactivation::Property::LateralStressComponentY ) { double s22 = interpolatedResultValue( iaDataAccess, m_femPart, position, s22Data ); double s33 = interpolatedResultValue( iaDataAccess, m_femPart, position, s33Data ); - double porBar = interpolatedResultValue( iaDataAccess, m_femPart, position, porData ); + double porBar = getPorBar( iaDataAccess, m_femPart, position, porData, *grid, elementIndex ); return ( s22 - porBar ) / ( s33 - porBar ); } } @@ -172,3 +178,104 @@ double RimFaultReactivationDataAccessorStress::interpolatedResultValue( RimWellI { return iaDataAccess.interpolatedResultValue( femPart, scalarResults, RIG_ELEMENT_NODAL, position ); } + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +double RimFaultReactivationDataAccessorStress::getPorBar( RimWellIADataAccess& iaDataAccess, + const RigFemPart* femPart, + const cvf::Vec3d& position, + const std::vector& porData, + const RigGriddedPart3d& grid, + size_t elementIndex ) const +{ + // Find element set for the position + const std::map>& elementSets = grid.elementSets(); + auto [isOk, elementSet] = findElementSetContainingElement( elementSets, elementIndex ); + + if ( !isOk ) + { + printf( "No data for pos: [%f, %f, %f] %zu.\n", position.x(), position.y(), position.z(), elementIndex ); + return std::numeric_limits::infinity(); + } + + // + std::vector> layers = grid.layers( elementSet ); + + if ( elementSet == RimFaultReactivation::ElementSets::OverBurden ) + { + // Find top reservoir pore pressure + std::vector> reservoirLayers = grid.layers( RimFaultReactivation::ElementSets::Reservoir ); + if ( reservoirLayers.empty() ) + { + printf( "MISSING RESERVOIR LAYER FOR OVERBURDEN POS: %f\n", position.z() ); + return std::numeric_limits::infinity(); + } + + cvf::Vec3d reservoirTopPos( position.x(), position.y(), reservoirLayers[0].first ); + double value = interpolatedResultValue( iaDataAccess, femPart, reservoirTopPos, porData ); + printf( "Overburden: depth: %f. data depth: %f (distance: %f) Value: %g\n", + position.z(), + reservoirTopPos.z(), + position.z() - reservoirTopPos.z(), + value ); + return value; + } + else if ( elementSet == RimFaultReactivation::ElementSets::Reservoir ) + { + double value = interpolatedResultValue( iaDataAccess, femPart, position, porData ); + printf( "Reservoir: depth: %f. Value: %g\n", position.z(), value ); + return value; + } + else if ( elementSet == RimFaultReactivation::ElementSets::IntraReservoir ) + { + double value = interpolatedResultValue( iaDataAccess, femPart, position, porData ); + printf( "Intra-reservoir: depth: %f. Value: %g\n", position.z(), value ); + + return value; + } + else if ( elementSet == RimFaultReactivation::ElementSets::UnderBurden ) + { + // Find top reservoir pore pressure + std::vector> reservoirLayers = grid.layers( RimFaultReactivation::ElementSets::Reservoir ); + if ( reservoirLayers.empty() ) + { + printf( "MISSING RESERVOIR LAYER FOR UNDERBURDEN POS: %f\n", position.z() ); + return std::numeric_limits::infinity(); + } + + cvf::Vec3d reservoirPos( position.x(), position.y(), reservoirLayers.back().second ); + double value = interpolatedResultValue( iaDataAccess, femPart, reservoirPos, porData ); + printf( "Underburden: depth: %f. data depth: %f (distance: %f) Value: %g\n", + position.z(), + reservoirPos.z(), + position.z() - reservoirPos.z(), + value ); + return value; + } + + return interpolatedResultValue( iaDataAccess, m_femPart, position, porData ); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +std::pair RimFaultReactivationDataAccessorStress::findElementSetContainingElement( + const std::map>& elementSets, + unsigned int elmIdx ) +{ + // Iterate through the map to find the ElementSets for the given unsigned int + for ( auto [element, elementSet] : elementSets ) + { + // Check if the targetValue is in the current vector + auto it = std::find( elementSet.begin(), elementSet.end(), elmIdx ); + + // If found, print the ElementSets + if ( it != elementSet.end() ) + { + return { true, element }; + } + } + + return { false, RimFaultReactivation::ElementSets::OverBurden }; +} diff --git a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStress.h b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStress.h index 4fdcfa4bd22..df8bb184123 100644 --- a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStress.h +++ b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStress.h @@ -25,11 +25,12 @@ #include +class RigFemPart; class RigFemPartResultsCollection; -class RimGeoMechCase; -class RigGeoMechCaseData; class RigFemScalarResultFrames; -class RigFemPart; +class RigGeoMechCaseData; +class RigGriddedPart3d; +class RimGeoMechCase; class RimWellIADataAccess; //================================================================================================== @@ -47,8 +48,9 @@ class RimFaultReactivationDataAccessorStress : public RimFaultReactivationDataAc double valueAtPosition( const cvf::Vec3d& position, const RigFaultReactivationModel& model, RimFaultReactivation::GridPart gridPart, - double topDepth = std::numeric_limits::infinity(), - double bottomDepth = std::numeric_limits::infinity() ) const override; + double topDepth = std::numeric_limits::infinity(), + double bottomDepth = std::numeric_limits::infinity(), + size_t elementIndex = std::numeric_limits::max() ) const override; private: void updateResultAccessor() override; @@ -60,6 +62,17 @@ class RimFaultReactivationDataAccessorStress : public RimFaultReactivationDataAc const cvf::Vec3d& position, const std::vector& scalarResults ) const; + double getPorBar( RimWellIADataAccess& iaDataAccess, + const RigFemPart* femPart, + const cvf::Vec3d& position, + const std::vector& porData, + const RigGriddedPart3d& grid, + size_t elementIndex ) const; + + static std::pair + findElementSetContainingElement( const std::map>& elementSets, + unsigned int elmIdx ); + RimGeoMechCase* m_geoMechCase; RimFaultReactivation::Property m_property; RigGeoMechCaseData* m_geoMechCaseData; diff --git a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorTemperature.cpp b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorTemperature.cpp index 10607338c2f..52c5eef8bfe 100644 --- a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorTemperature.cpp +++ b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorTemperature.cpp @@ -87,7 +87,8 @@ double RimFaultReactivationDataAccessorTemperature::valueAtPosition( const cvf:: const RigFaultReactivationModel& model, RimFaultReactivation::GridPart gridPart, double topDepth, - double bottomDepth ) const + double bottomDepth, + size_t elementIndex ) const { if ( ( m_mainGrid != nullptr ) && m_resultAccessor.notNull() ) { diff --git a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorTemperature.h b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorTemperature.h index 8cb50c065cc..3d310152d6f 100644 --- a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorTemperature.h +++ b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorTemperature.h @@ -44,8 +44,9 @@ class RimFaultReactivationDataAccessorTemperature : public RimFaultReactivationD double valueAtPosition( const cvf::Vec3d& position, const RigFaultReactivationModel& model, RimFaultReactivation::GridPart gridPart, - double topDepth = std::numeric_limits::infinity(), - double bottomDepth = std::numeric_limits::infinity() ) const override; + double topDepth = std::numeric_limits::infinity(), + double bottomDepth = std::numeric_limits::infinity(), + size_t elementIndex = std::numeric_limits::max() ) const override; private: void updateResultAccessor() override; diff --git a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorVoidRatio.cpp b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorVoidRatio.cpp index c1dc6ccfb3e..96fe0c6f3be 100644 --- a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorVoidRatio.cpp +++ b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorVoidRatio.cpp @@ -88,7 +88,8 @@ double RimFaultReactivationDataAccessorVoidRatio::valueAtPosition( const cvf::Ve const RigFaultReactivationModel& model, RimFaultReactivation::GridPart gridPart, double topDepth, - double bottomDepth ) const + double bottomDepth, + size_t elementIndex ) const { if ( ( m_mainGrid != nullptr ) && m_resultAccessor.notNull() ) { diff --git a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorVoidRatio.h b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorVoidRatio.h index 069b3b344e0..519a5727da4 100644 --- a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorVoidRatio.h +++ b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorVoidRatio.h @@ -44,8 +44,9 @@ class RimFaultReactivationDataAccessorVoidRatio : public RimFaultReactivationDat double valueAtPosition( const cvf::Vec3d& position, const RigFaultReactivationModel& model, RimFaultReactivation::GridPart gridPart, - double topDepth = std::numeric_limits::infinity(), - double bottomDepth = std::numeric_limits::infinity() ) const override; + double topDepth = std::numeric_limits::infinity(), + double bottomDepth = std::numeric_limits::infinity(), + size_t elementIndex = std::numeric_limits::max() ) const override; private: void updateResultAccessor() override;