From adf52ed3018c8963f136680899652bfd077dd190 Mon Sep 17 00:00:00 2001 From: Kristian Bendiksen Date: Fri, 5 Jan 2024 14:59:05 +0100 Subject: [PATCH] #10827 Fault reactivation: extract temperatures using well log. --- .../RifFaultReactivationModelExporter.cpp | 3 +- .../Faults/RimFaultReactivationDataAccess.cpp | 17 +- .../Faults/RimFaultReactivationDataAccess.h | 6 +- ...ltReactivationDataAccessorPorePressure.cpp | 11 +- ...aultReactivationDataAccessorPorePressure.h | 3 +- ...RimFaultReactivationDataAccessorStress.cpp | 6 +- .../RimFaultReactivationDataAccessorStress.h | 6 +- ...ultReactivationDataAccessorTemperature.cpp | 48 +++-- ...FaultReactivationDataAccessorTemperature.h | 16 +- ...ctivationDataAccessorWellLogExtraction.cpp | 170 +++++++++++++----- ...eactivationDataAccessorWellLogExtraction.h | 40 ++++- .../Faults/RimFaultReactivationModel.cpp | 14 ++ .../Faults/RimFaultReactivationModel.h | 2 + 13 files changed, 256 insertions(+), 86 deletions(-) diff --git a/ApplicationLibCode/FileInterface/RifFaultReactivationModelExporter.cpp b/ApplicationLibCode/FileInterface/RifFaultReactivationModelExporter.cpp index a6d08042c0..f9b7163378 100644 --- a/ApplicationLibCode/FileInterface/RifFaultReactivationModelExporter.cpp +++ b/ApplicationLibCode/FileInterface/RifFaultReactivationModelExporter.cpp @@ -850,7 +850,8 @@ std::shared_ptr if ( eCase == nullptr ) return nullptr; // extract data for each timestep - auto dataAccess = std::make_shared( eCase, rimModel.geoMechCase(), rimModel.selectedTimeStepIndexes() ); + auto dataAccess = + std::make_shared( rimModel, eCase, rimModel.geoMechCase(), rimModel.selectedTimeStepIndexes() ); dataAccess->extractModelData( *rimModel.model() ); return dataAccess; } diff --git a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccess.cpp b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccess.cpp index 5ff6135df9..44275676aa 100644 --- a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccess.cpp +++ b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccess.cpp @@ -39,19 +39,23 @@ #include "RimFaultReactivationDataAccessorTemperature.h" #include "RimFaultReactivationDataAccessorVoidRatio.h" #include "RimFaultReactivationEnums.h" +#include "RimFaultReactivationModel.h" //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- -RimFaultReactivationDataAccess::RimFaultReactivationDataAccess( RimEclipseCase* thecase, - RimGeoMechCase* geoMechCase, - const std::vector& timeSteps ) +RimFaultReactivationDataAccess::RimFaultReactivationDataAccess( const RimFaultReactivationModel& model, + RimEclipseCase* thecase, + RimGeoMechCase* geoMechCase, + const std::vector& timeSteps ) : m_timeSteps( timeSteps ) { double porePressureGradient = 1.0; - m_accessors.push_back( std::make_shared( thecase, porePressureGradient ) ); + double topTemperature = model.seabedTemperature(); + double seabedDepth = -model.seaBedDepth(); + m_accessors.push_back( std::make_shared( thecase, porePressureGradient, seabedDepth ) ); m_accessors.push_back( std::make_shared( thecase, 0.0001 ) ); - m_accessors.push_back( std::make_shared( thecase ) ); + m_accessors.push_back( std::make_shared( thecase, topTemperature, seabedDepth ) ); if ( geoMechCase ) { std::vector properties = { RimFaultReactivation::Property::YoungsModulus, @@ -71,7 +75,8 @@ RimFaultReactivationDataAccess::RimFaultReactivationDataAccess( RimEclipseCase* RimFaultReactivation::Property::LateralStressComponentY }; for ( auto property : stressProperties ) { - m_accessors.push_back( std::make_shared( geoMechCase, property, porePressureGradient ) ); + m_accessors.push_back( + std::make_shared( geoMechCase, property, porePressureGradient, seabedDepth ) ); } } } diff --git a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccess.h b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccess.h index 07f5b74014..0ede2cf555 100644 --- a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccess.h +++ b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccess.h @@ -32,6 +32,7 @@ class RimEclipseCase; class RimGeoMechCase; class RimFaultReactivationDataAccessor; class RigFaultReactivationModel; +class RimFaultReactivationModel; //================================================================================================== /// @@ -40,7 +41,10 @@ class RigFaultReactivationModel; class RimFaultReactivationDataAccess { public: - RimFaultReactivationDataAccess( RimEclipseCase* eclipseCase, RimGeoMechCase* geoMechCase, const std::vector& timeSteps ); + RimFaultReactivationDataAccess( const RimFaultReactivationModel& model, + RimEclipseCase* eclipseCase, + RimGeoMechCase* geoMechCase, + const std::vector& timeSteps ); ~RimFaultReactivationDataAccess(); void extractModelData( const RigFaultReactivationModel& model ); diff --git a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorPorePressure.cpp b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorPorePressure.cpp index 2cd7edc73e..487c937ea0 100644 --- a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorPorePressure.cpp +++ b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorPorePressure.cpp @@ -42,9 +42,11 @@ /// //-------------------------------------------------------------------------------------------------- RimFaultReactivationDataAccessorPorePressure::RimFaultReactivationDataAccessorPorePressure( RimEclipseCase* eclipseCase, - double porePressureGradient ) + double porePressureGradient, + double seabedDepth ) : m_eclipseCase( eclipseCase ) , m_defaultPorePressureGradient( porePressureGradient ) + , m_seabedDepth( seabedDepth ) , m_caseData( nullptr ) , m_mainGrid( nullptr ) { @@ -75,9 +77,10 @@ void RimFaultReactivationDataAccessorPorePressure::updateResultAccessor() m_resultAccessor = RigResultAccessorFactory::createFromResultAddress( m_caseData, 0, RiaDefines::PorosityModelType::MATRIX_MODEL, m_timeStep, resVarAddress ); - auto [wellPaths, extractors] = RimFaultReactivationDataAccessorWellLogExtraction::createEclipseWellPathExtractors( *m_model, *m_caseData ); - m_wellPaths = wellPaths; - m_extractors = extractors; + auto [wellPaths, extractors] = + RimFaultReactivationDataAccessorWellLogExtraction::createEclipseWellPathExtractors( *m_model, *m_caseData, m_seabedDepth ); + m_wellPaths = wellPaths; + m_extractors = extractors; } //-------------------------------------------------------------------------------------------------- diff --git a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorPorePressure.h b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorPorePressure.h index 43e0e25d0e..e8e759bbe7 100644 --- a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorPorePressure.h +++ b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorPorePressure.h @@ -37,7 +37,7 @@ class RigWellPath; class RimFaultReactivationDataAccessorPorePressure : public RimFaultReactivationDataAccessor { public: - RimFaultReactivationDataAccessorPorePressure( RimEclipseCase* eclipseCase, double porePressureGradient ); + RimFaultReactivationDataAccessorPorePressure( RimEclipseCase* eclipseCase, double porePressureGradient, double seabedDepth ); ~RimFaultReactivationDataAccessorPorePressure(); bool isMatching( RimFaultReactivation::Property property ) const override; @@ -58,6 +58,7 @@ private: RigEclipseCaseData* m_caseData; const RigMainGrid* m_mainGrid; double m_defaultPorePressureGradient; + double m_seabedDepth; cvf::ref m_resultAccessor; std::map> m_wellPaths; diff --git a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStress.cpp b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStress.cpp index 8bcca8b66c..9dc4b5afc2 100644 --- a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStress.cpp +++ b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStress.cpp @@ -48,10 +48,12 @@ //-------------------------------------------------------------------------------------------------- RimFaultReactivationDataAccessorStress::RimFaultReactivationDataAccessorStress( RimGeoMechCase* geoMechCase, RimFaultReactivation::Property property, - double gradient ) + double gradient, + double seabedDepth ) : m_geoMechCase( geoMechCase ) , m_property( property ) , m_gradient( gradient ) + , m_seabedDepth( seabedDepth ) { m_geoMechCaseData = geoMechCase->geoMechData(); } @@ -99,6 +101,7 @@ void RimFaultReactivationDataAccessorStress::updateResultAccessor() std::vector wellPoints = RimFaultReactivationDataAccessorWellLogExtraction::generateWellPoints( faultTopPosition, faultBottomPosition, + m_seabedDepth, faultNormal * distanceFromFault ); m_faceAWellPath = new RigWellPath( wellPoints, RimFaultReactivationDataAccessorWellLogExtraction::generateMds( wellPoints ) ); m_partIndexA = geoMechPartCollection->getPartIndexFromPoint( wellPoints[1] ); @@ -109,6 +112,7 @@ void RimFaultReactivationDataAccessorStress::updateResultAccessor() std::vector wellPoints = RimFaultReactivationDataAccessorWellLogExtraction::generateWellPoints( faultTopPosition, faultBottomPosition, + m_seabedDepth, -faultNormal * distanceFromFault ); m_faceBWellPath = new RigWellPath( wellPoints, RimFaultReactivationDataAccessorWellLogExtraction::generateMds( wellPoints ) ); m_partIndexB = geoMechPartCollection->getPartIndexFromPoint( wellPoints[1] ); diff --git a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStress.h b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStress.h index 056048ef98..d14c8147a6 100644 --- a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStress.h +++ b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStress.h @@ -49,7 +49,10 @@ class RigGeoMechWellLogExtractor; class RimFaultReactivationDataAccessorStress : public RimFaultReactivationDataAccessor { public: - RimFaultReactivationDataAccessorStress( RimGeoMechCase* geoMechCase, RimFaultReactivation::Property property, double gradient ); + RimFaultReactivationDataAccessorStress( RimGeoMechCase* geoMechCase, + RimFaultReactivation::Property property, + double gradient, + double seabedDepth ); ~RimFaultReactivationDataAccessorStress(); bool isMatching( RimFaultReactivation::Property property ) const override; @@ -76,6 +79,7 @@ private: RimGeoMechCase* m_geoMechCase; RimFaultReactivation::Property m_property; double m_gradient; + double m_seabedDepth; RigGeoMechCaseData* m_geoMechCaseData; RigFemScalarResultFrames* m_s11Frames; RigFemScalarResultFrames* m_s22Frames; diff --git a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorTemperature.cpp b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorTemperature.cpp index 52c5eef8bf..3ce6c0f4b5 100644 --- a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorTemperature.cpp +++ b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorTemperature.cpp @@ -25,10 +25,13 @@ #include "RigCaseCellResultsData.h" #include "RigEclipseCaseData.h" #include "RigEclipseResultAddress.h" +#include "RigEclipseWellLogExtractor.h" #include "RigMainGrid.h" #include "RigResultAccessorFactory.h" +#include "RigWellPath.h" #include "RimEclipseCase.h" +#include "RimFaultReactivationDataAccessorWellLogExtraction.h" #include #include @@ -36,10 +39,14 @@ //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- -RimFaultReactivationDataAccessorTemperature::RimFaultReactivationDataAccessorTemperature( RimEclipseCase* eclipseCase ) +RimFaultReactivationDataAccessorTemperature::RimFaultReactivationDataAccessorTemperature( RimEclipseCase* eclipseCase, + double seabedTemperature, + double seabedDepth ) : m_eclipseCase( eclipseCase ) , m_caseData( nullptr ) , m_mainGrid( nullptr ) + , m_seabedTemperature( seabedTemperature ) + , m_seabedDepth( seabedDepth ) { if ( m_eclipseCase ) { @@ -60,16 +67,17 @@ RimFaultReactivationDataAccessorTemperature::~RimFaultReactivationDataAccessorTe //-------------------------------------------------------------------------------------------------- void RimFaultReactivationDataAccessorTemperature::updateResultAccessor() { - if ( m_caseData ) - { - RigEclipseResultAddress resVarAddress( RiaDefines::ResultCatType::DYNAMIC_NATIVE, "TEMP" ); - m_eclipseCase->results( RiaDefines::PorosityModelType::MATRIX_MODEL )->ensureKnownResultLoaded( resVarAddress ); - m_resultAccessor = RigResultAccessorFactory::createFromResultAddress( m_caseData, - 0, - RiaDefines::PorosityModelType::MATRIX_MODEL, - m_timeStep, - resVarAddress ); - } + if ( !m_caseData ) return; + + RigEclipseResultAddress resVarAddress( RiaDefines::ResultCatType::DYNAMIC_NATIVE, "SOIL" ); + m_eclipseCase->results( RiaDefines::PorosityModelType::MATRIX_MODEL )->ensureKnownResultLoaded( resVarAddress ); + m_resultAccessor = + RigResultAccessorFactory::createFromResultAddress( m_caseData, 0, RiaDefines::PorosityModelType::MATRIX_MODEL, m_timeStep, resVarAddress ); + + auto [wellPaths, extractors] = + RimFaultReactivationDataAccessorWellLogExtraction::createEclipseWellPathExtractors( *m_model, *m_caseData, m_seabedDepth ); + m_wellPaths = wellPaths; + m_extractors = extractors; } //-------------------------------------------------------------------------------------------------- @@ -92,11 +100,19 @@ double RimFaultReactivationDataAccessorTemperature::valueAtPosition( const cvf:: { if ( ( m_mainGrid != nullptr ) && m_resultAccessor.notNull() ) { - auto cellIdx = m_mainGrid->findReservoirCellIndexFromPoint( position ); - if ( cellIdx != cvf::UNDEFINED_SIZE_T ) - { - return m_resultAccessor->cellScalar( cellIdx ); - } + CAF_ASSERT( m_extractors.find( gridPart ) != m_extractors.end() ); + auto extractor = m_extractors.find( gridPart )->second; + + CAF_ASSERT( m_wellPaths.find( gridPart ) != m_wellPaths.end() ); + auto wellPath = m_wellPaths.find( gridPart )->second; + + auto [values, intersections] = + RimFaultReactivationDataAccessorWellLogExtraction::extractValuesAndIntersections( *m_resultAccessor.p(), *extractor.p(), *wellPath ); + + auto [value, pos] = + RimFaultReactivationDataAccessorWellLogExtraction::calculateTemperature( intersections, values, position, m_seabedTemperature ); + + return value; } return std::numeric_limits::infinity(); diff --git a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorTemperature.h b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorTemperature.h index 3d310152d6..058d23cdba 100644 --- a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorTemperature.h +++ b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorTemperature.h @@ -28,6 +28,8 @@ class RimEclipseCase; class RigEclipseCaseData; class RigMainGrid; class RigResultAccessor; +class RigWellPath; +class RigEclipseWellLogExtractor; //================================================================================================== /// @@ -36,7 +38,7 @@ class RigResultAccessor; class RimFaultReactivationDataAccessorTemperature : public RimFaultReactivationDataAccessor { public: - RimFaultReactivationDataAccessorTemperature( RimEclipseCase* eclipseCase ); + RimFaultReactivationDataAccessorTemperature( RimEclipseCase* eclipseCase, double seabedTemperature, double seabedDepth ); ~RimFaultReactivationDataAccessorTemperature(); bool isMatching( RimFaultReactivation::Property property ) const override; @@ -51,8 +53,14 @@ public: private: void updateResultAccessor() override; - RimEclipseCase* m_eclipseCase; - RigEclipseCaseData* m_caseData; - const RigMainGrid* m_mainGrid; + RimEclipseCase* m_eclipseCase; + RigEclipseCaseData* m_caseData; + const RigMainGrid* m_mainGrid; + double m_seabedTemperature; + double m_seabedDepth; + cvf::ref m_resultAccessor; + + std::map> m_wellPaths; + std::map> m_extractors; }; diff --git a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorWellLogExtraction.cpp b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorWellLogExtraction.cpp index d83aadbb97..1e4d088aa5 100644 --- a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorWellLogExtraction.cpp +++ b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorWellLogExtraction.cpp @@ -63,13 +63,37 @@ RimFaultReactivationDataAccessorWellLogExtraction::~RimFaultReactivationDataAcce /// //-------------------------------------------------------------------------------------------------- std::pair RimFaultReactivationDataAccessorWellLogExtraction::calculatePorBar( const std::vector& intersections, - std::vector& values, - const cvf::Vec3d& position, - double gradient ) + std::vector& values, + const cvf::Vec3d& position, + double gradient ) { // Fill in missing values - fillInMissingValues( intersections, values, gradient ); + fillInMissingValuesWithGradient( intersections, values, gradient ); + return findValueAndPosition( intersections, values, position ); +} +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +std::pair + RimFaultReactivationDataAccessorWellLogExtraction::calculateTemperature( const std::vector& intersections, + std::vector& values, + const cvf::Vec3d& position, + double seabedTemperature ) +{ + // Fill in missing values + fillInMissingValuesWithTopValue( intersections, values, seabedTemperature ); + return findValueAndPosition( intersections, values, position ); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +std::pair + RimFaultReactivationDataAccessorWellLogExtraction::findValueAndPosition( const std::vector& intersections, + const std::vector& values, + const cvf::Vec3d& position ) +{ // Linear interpolation between two points auto lerp = []( const cvf::Vec3d& start, const cvf::Vec3d& end, double t ) { return start + t * ( end - start ); }; @@ -99,7 +123,7 @@ std::pair RimFaultReactivationDataAccessorWellLogExtraction: /// //-------------------------------------------------------------------------------------------------- std::pair RimFaultReactivationDataAccessorWellLogExtraction::findIntersectionsForTvd( const std::vector& intersections, - double tvd ) + double tvd ) { int topIdx = -1; int bottomIdx = -1; @@ -109,7 +133,7 @@ std::pair RimFaultReactivationDataAccessorWellLogExtraction::findInter { auto top = intersections[i - 1]; auto bottom = intersections[i]; - if ( top.z() > tvd && bottom.z() < tvd ) + if ( top.z() >= tvd && bottom.z() < tvd ) { topIdx = static_cast( i ) - 1; bottomIdx = static_cast( i ); @@ -153,11 +177,11 @@ std::pair RimFaultReactivationDataAccessorWellLogExtraction::findOverb //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- -double RimFaultReactivationDataAccessorWellLogExtraction::computePorBarWithGradient( const std::vector& intersections, - const std::vector& values, - int i1, - int i2, - double gradient ) +double RimFaultReactivationDataAccessorWellLogExtraction::computeValueWithGradient( const std::vector& intersections, + const std::vector& values, + int i1, + int i2, + double gradient ) { double tvdDiff = intersections[i2].z() - intersections[i1].z(); return tvdDiff * gradient + values[i2]; @@ -166,50 +190,51 @@ double RimFaultReactivationDataAccessorWellLogExtraction::computePorBarWithGradi //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- -void RimFaultReactivationDataAccessorWellLogExtraction::fillInMissingValues( const std::vector& intersections, - std::vector& values, - double gradient ) +void RimFaultReactivationDataAccessorWellLogExtraction::fillInMissingValuesWithGradient( const std::vector& intersections, + std::vector& values, + double gradient ) { CAF_ASSERT( intersections.size() == values.size() ); auto calculatePorePressure = []( double depth, double gradient ) { return RiaEclipseUnitTools::pascalToBar( gradient * 9.81 * depth * 1000.0 ); }; - auto computeGradient = []( double depth1, double value1, double depth2, double value2 ) - { return ( value2 - value1 ) / ( depth2 - depth1 ); }; - auto [lastOverburdenIndex, firstUnderburdenIndex] = findOverburdenAndUnderburdenIndex( values ); // Fill in overburden values using gradient double topPorePressure = calculatePorePressure( std::abs( intersections[0].z() ), gradient ); - double overburdenGradient = - computeGradient( intersections[0].z(), topPorePressure, intersections[lastOverburdenIndex].z(), values[lastOverburdenIndex] ); - - for ( int i = 0; i < lastOverburdenIndex; i++ ) - { - values[i] = computePorBarWithGradient( intersections, values, i, lastOverburdenIndex, -overburdenGradient ); - } + insertOverburdenValues( intersections, values, lastOverburdenIndex, topPorePressure ); // Fill in underburden values using gradient - int lastElementIndex = static_cast( values.size() ) - 1; - double bottomPorePressure = calculatePorePressure( std::abs( intersections[lastElementIndex].z() ), gradient ); - double underburdenGradient = computeGradient( intersections[firstUnderburdenIndex].z(), - values[firstUnderburdenIndex], - intersections[lastElementIndex].z(), - bottomPorePressure ); - - for ( int i = lastElementIndex; i >= firstUnderburdenIndex; i-- ) - { - values[i] = computePorBarWithGradient( intersections, values, i, firstUnderburdenIndex, -underburdenGradient ); - } + int lastElementIndex = static_cast( values.size() ) - 1; + double bottomPorePressure = calculatePorePressure( std::abs( intersections[lastElementIndex].z() ), gradient ); + insertUnderburdenValues( intersections, values, firstUnderburdenIndex, bottomPorePressure ); // Interpolate the missing values (should only be intra-reservoir by now) - std::vector intersectionsZ; - for ( auto i : intersections ) - { - intersectionsZ.push_back( i.z() ); - } + std::vector intersectionsZ = extractDepthValues( intersections ); + RiaInterpolationTools::interpolateMissingValues( intersectionsZ, values ); +} +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +void RimFaultReactivationDataAccessorWellLogExtraction::fillInMissingValuesWithTopValue( const std::vector& intersections, + std::vector& values, + double topValue ) +{ + CAF_ASSERT( intersections.size() == values.size() ); + + auto [lastOverburdenIndex, firstUnderburdenIndex] = findOverburdenAndUnderburdenIndex( values ); + + // Fill in overburden values using gradient + insertOverburdenValues( intersections, values, lastOverburdenIndex, topValue ); + + // Fill in underburden values + double bottomValue = values[firstUnderburdenIndex - 1]; + insertUnderburdenValues( intersections, values, firstUnderburdenIndex, bottomValue ); + + // Interpolate the missing values (should only be intra-reservoir by now) + std::vector intersectionsZ = extractDepthValues( intersections ); RiaInterpolationTools::interpolateMissingValues( intersectionsZ, values ); } @@ -218,10 +243,11 @@ void RimFaultReactivationDataAccessorWellLogExtraction::fillInMissingValues( con //-------------------------------------------------------------------------------------------------- std::vector RimFaultReactivationDataAccessorWellLogExtraction::generateWellPoints( const cvf::Vec3d& faultTopPosition, const cvf::Vec3d& faultBottomPosition, + double seabedDepth, const cvf::Vec3d& offset ) { cvf::Vec3d faultTop = faultTopPosition + offset; - cvf::Vec3d seabed( faultTop.x(), faultTop.y(), 0.0 ); + cvf::Vec3d seabed( faultTop.x(), faultTop.y(), seabedDepth ); cvf::Vec3d faultBottom = faultBottomPosition + offset; cvf::Vec3d underburdenBottom( faultBottom.x(), faultBottom.y(), -10000.0 ); return { seabed, faultTop, faultBottom, underburdenBottom }; @@ -250,7 +276,8 @@ std::vector RimFaultReactivationDataAccessorWellLogExtraction::generateM //-------------------------------------------------------------------------------------------------- std::pair>, std::map>> RimFaultReactivationDataAccessorWellLogExtraction::createEclipseWellPathExtractors( const RigFaultReactivationModel& model, - RigEclipseCaseData& eclipseCaseData ) + RigEclipseCaseData& eclipseCaseData, + double seabedDepth ) { auto [faultTopPosition, faultBottomPosition] = model.faultTopBottom(); auto faultNormal = model.faultNormal(); @@ -265,6 +292,7 @@ std::pair>, std:: std::vector wellPoints = RimFaultReactivationDataAccessorWellLogExtraction::generateWellPoints( faultTopPosition, faultBottomPosition, + seabedDepth, sign * faultNormal * distanceFromFault ); cvf::ref wellPath = new RigWellPath( wellPoints, RimFaultReactivationDataAccessorWellLogExtraction::generateMds( wellPoints ) ); @@ -302,3 +330,61 @@ std::pair, std::vector> return { values, intersections }; } + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +double RimFaultReactivationDataAccessorWellLogExtraction::computeGradient( double depth1, double value1, double depth2, double value2 ) +{ + return ( value2 - value1 ) / ( depth2 - depth1 ); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +std::vector RimFaultReactivationDataAccessorWellLogExtraction::extractDepthValues( const std::vector& intersections ) +{ + std::vector intersectionsZ; + for ( auto i : intersections ) + { + intersectionsZ.push_back( i.z() ); + } + return intersectionsZ; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +void RimFaultReactivationDataAccessorWellLogExtraction::insertUnderburdenValues( const std::vector& intersections, + std::vector& values, + int firstUnderburdenIndex, + double bottomValue ) +{ + int lastElementIndex = static_cast( values.size() ) - 1; + double underburdenGradient = computeGradient( intersections[firstUnderburdenIndex].z(), + values[firstUnderburdenIndex], + intersections[lastElementIndex].z(), + bottomValue ); + + for ( int i = lastElementIndex; i >= firstUnderburdenIndex; i-- ) + { + values[i] = computeValueWithGradient( intersections, values, i, firstUnderburdenIndex, -underburdenGradient ); + } +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +void RimFaultReactivationDataAccessorWellLogExtraction::insertOverburdenValues( const std::vector& intersections, + std::vector& values, + int lastOverburdenIndex, + double topValue ) +{ + double overburdenGradient = + computeGradient( intersections[0].z(), topValue, intersections[lastOverburdenIndex].z(), values[lastOverburdenIndex] ); + + for ( int i = 0; i < lastOverburdenIndex; i++ ) + { + values[i] = computeValueWithGradient( intersections, values, i, lastOverburdenIndex, -overburdenGradient ); + } +} diff --git a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorWellLogExtraction.h b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorWellLogExtraction.h index a7611ae6d0..94f518204f 100644 --- a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorWellLogExtraction.h +++ b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorWellLogExtraction.h @@ -43,12 +43,18 @@ public: static std::pair calculatePorBar( const std::vector& intersections, std::vector& values, const cvf::Vec3d& position, double gradient ); + static std::pair calculateTemperature( const std::vector& intersections, + std::vector& values, + const cvf::Vec3d& position, + double seabedTemperature ); static std::pair>, std::map>> - createEclipseWellPathExtractors( const RigFaultReactivationModel& model, RigEclipseCaseData& eclipseCaseData ); + createEclipseWellPathExtractors( const RigFaultReactivationModel& model, RigEclipseCaseData& eclipseCaseData, double seabedDepth ); - static std::vector - generateWellPoints( const cvf::Vec3d& faultTopPosition, const cvf::Vec3d& faultBottomPosition, const cvf::Vec3d& offset ); + static std::vector generateWellPoints( const cvf::Vec3d& faultTopPosition, + const cvf::Vec3d& faultBottomPosition, + double seabedDepth, + const cvf::Vec3d& offset ); static std::vector generateMds( const std::vector& points ); @@ -59,10 +65,26 @@ public: protected: static std::pair findIntersectionsForTvd( const std::vector& intersections, double tvd ); static std::pair findOverburdenAndUnderburdenIndex( const std::vector& values ); - static double computePorBarWithGradient( const std::vector& intersections, - const std::vector& values, - int i1, - int i2, - double gradient ); - static void fillInMissingValues( const std::vector& intersections, std::vector& values, double gradient ); + static double computeValueWithGradient( const std::vector& intersections, + const std::vector& values, + int i1, + int i2, + double gradient ); + static void fillInMissingValuesWithGradient( const std::vector& intersections, std::vector& values, double gradient ); + static void fillInMissingValuesWithTopValue( const std::vector& intersections, std::vector& values, double topValue ); + + static std::pair + findValueAndPosition( const std::vector& intersections, const std::vector& values, const cvf::Vec3d& position ); + + static double computeGradient( double depth1, double value1, double depth2, double value2 ); + static std::vector extractDepthValues( const std::vector& intersections ); + + static void insertUnderburdenValues( const std::vector& intersections, + std::vector& values, + int firstUnderburdenIndex, + double bottomValue ); + static void insertOverburdenValues( const std::vector& intersections, + std::vector& values, + int lastOverburdenIndex, + double topValue ); }; diff --git a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationModel.cpp b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationModel.cpp index b71487b435..8b02ef0baa 100644 --- a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationModel.cpp +++ b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationModel.cpp @@ -127,6 +127,7 @@ RimFaultReactivationModel::RimFaultReactivationModel() CAF_PDM_InitField( &m_waterDensity, "WaterDensity", 1030.0, "Water Density [kg/m3]" ); CAF_PDM_InitField( &m_frictionAngleDeg, "FrictionAngle", 20.0, "Friction Angle [degree]" ); + CAF_PDM_InitField( &m_seabedTemperature, "SeabedTemperature", 5.0, "Seabed Temperature [C]" ); CAF_PDM_InitFieldNoDefault( &m_targets, "Targets", "Targets" ); m_targets.uiCapability()->setUiEditorTypeName( caf::PdmUiTableViewEditor::uiEditorTypeName() ); @@ -480,6 +481,11 @@ void RimFaultReactivationModel::defineUiOrdering( QString uiConfigName, caf::Pdm propertiesGrp->add( &m_useGridElasticProperties ); propertiesGrp->add( &m_useGridStress ); propertiesGrp->add( &m_waterDensity ); + + propertiesGrp->add( &m_seabedTemperature ); + + bool useTemperatureFromGrid = m_useGridTemperature(); + m_seabedTemperature.uiCapability()->setUiReadOnly( !useTemperatureFromGrid ); } propertiesGrp->add( &m_frictionAngleDeg ); @@ -790,3 +796,11 @@ double RimFaultReactivationModel::waterDensity() const { return m_waterDensity; } + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +double RimFaultReactivationModel::seabedTemperature() const +{ + return m_seabedTemperature; +} diff --git a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationModel.h b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationModel.h index 37df7b424b..4785b55a40 100644 --- a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationModel.h +++ b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationModel.h @@ -128,6 +128,7 @@ public: double seaBedDepth() const; double waterDensity() const; double frictionAngleDeg() const; + double seabedTemperature() const; RimEclipseCase* eclipseCase() const; RimGeoMechCase* geoMechCase() const; @@ -182,6 +183,7 @@ private: caf::PdmField m_waterDensity; caf::PdmField m_frictionAngleDeg; + caf::PdmField m_seabedTemperature; caf::PdmField m_startCellIndex; caf::PdmField m_startCellFace;