Fault Reactivation: Improve temperature data access.

New approach: find the top reservoir and compute gradient of temperature change
from sea bed. Values from Eclipse grid is used when in reservoir.
This commit is contained in:
Kristian Bendiksen 2024-02-14 14:02:05 +01:00
parent 85753568b6
commit d8b1e643e0
4 changed files with 86 additions and 39 deletions

View File

@ -17,6 +17,7 @@
/////////////////////////////////////////////////////////////////////////////////
#include "RimFaultReactivationDataAccessorTemperature.h"
#include "RigFaultReactivationModel.h"
#include "RimFaultReactivationEnums.h"
#include "RiaDefines.h"
@ -78,6 +79,48 @@ void RimFaultReactivationDataAccessorTemperature::updateResultAccessor()
RimFaultReactivationDataAccessorWellLogExtraction::createEclipseWellPathExtractors( *m_model, *m_caseData, m_seabedDepth );
m_wellPaths = wellPaths;
m_extractors = extractors;
m_gradient = computeGradient();
}
//--------------------------------------------------------------------------------------------------
/// Find the top encounter with reservoir (of the two well paths), and create gradient from that point
//--------------------------------------------------------------------------------------------------
double RimFaultReactivationDataAccessorTemperature::computeGradient() const
{
double gradient = std::numeric_limits<double>::infinity();
double minDepth = -std::numeric_limits<double>::max();
for ( auto gridPart : m_model->allGridParts() )
{
auto extractor = m_extractors.find( gridPart )->second;
auto wellPath = m_wellPaths.find( gridPart )->second;
auto [values, intersections] =
RimFaultReactivationDataAccessorWellLogExtraction::extractValuesAndIntersections( *m_resultAccessor.p(), *extractor.p(), *wellPath );
int lastOverburdenIndex = RimFaultReactivationDataAccessorWellLogExtraction::findLastOverburdenIndex( values );
if ( lastOverburdenIndex != -1 )
{
double depth = intersections[lastOverburdenIndex].z();
double value = values[lastOverburdenIndex];
if ( !std::isinf( value ) )
{
double currentGradient =
RimFaultReactivationDataAccessorWellLogExtraction::computeGradient( intersections[0].z(),
m_seabedTemperature,
intersections[lastOverburdenIndex].z(),
values[lastOverburdenIndex] );
if ( !std::isinf( value ) && !std::isnan( currentGradient ) && depth > minDepth )
{
gradient = currentGradient;
minDepth = depth;
}
}
}
}
return gradient;
}
//--------------------------------------------------------------------------------------------------
@ -100,6 +143,13 @@ double RimFaultReactivationDataAccessorTemperature::valueAtPosition( const cvf::
{
if ( ( m_mainGrid != nullptr ) && m_resultAccessor.notNull() )
{
auto cellIdx = m_mainGrid->findReservoirCellIndexFromPoint( position );
if ( cellIdx != cvf::UNDEFINED_SIZE_T )
{
double tempFromEclipse = m_resultAccessor->cellScalar( cellIdx );
if ( !std::isinf( tempFromEclipse ) ) return tempFromEclipse;
}
CAF_ASSERT( m_extractors.find( gridPart ) != m_extractors.end() );
auto extractor = m_extractors.find( gridPart )->second;
@ -110,16 +160,7 @@ double RimFaultReactivationDataAccessorTemperature::valueAtPosition( const cvf::
RimFaultReactivationDataAccessorWellLogExtraction::extractValuesAndIntersections( *m_resultAccessor.p(), *extractor.p(), *wellPath );
auto [value, pos] =
RimFaultReactivationDataAccessorWellLogExtraction::calculateTemperature( intersections, values, position, m_seabedTemperature );
if ( pos.isUndefined() )
{
auto cellIdx = m_mainGrid->findReservoirCellIndexFromPoint( position );
if ( cellIdx != cvf::UNDEFINED_SIZE_T )
{
double tempFromEclipse = m_resultAccessor->cellScalar( cellIdx );
if ( !std::isinf( tempFromEclipse ) ) return tempFromEclipse;
}
}
RimFaultReactivationDataAccessorWellLogExtraction::calculateTemperature( intersections, position, m_seabedTemperature, m_gradient );
return value;
}

View File

@ -51,13 +51,15 @@ public:
size_t elementIndex = std::numeric_limits<size_t>::max() ) const override;
private:
void updateResultAccessor() override;
void updateResultAccessor() override;
double computeGradient() const;
RimEclipseCase* m_eclipseCase;
RigEclipseCaseData* m_caseData;
const RigMainGrid* m_mainGrid;
double m_seabedTemperature;
double m_seabedDepth;
double m_gradient;
cvf::ref<RigResultAccessor> m_resultAccessor;

View File

@ -96,23 +96,11 @@ std::pair<double, cvf::Vec3d> RimFaultReactivationDataAccessorWellLogExtraction:
//--------------------------------------------------------------------------------------------------
std::pair<double, cvf::Vec3d>
RimFaultReactivationDataAccessorWellLogExtraction::calculateTemperature( const std::vector<cvf::Vec3d>& intersections,
std::vector<double>& values,
const cvf::Vec3d& position,
double seabedTemperature )
double seabedTemperature,
double gradient )
{
// Fill in missing values
fillInMissingValuesWithTopValue( intersections, values, seabedTemperature );
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 };
}
return { calculateTemperature( seabedTemperature, intersections[0].z(), std::abs( position.z() ), gradient ), position };
}
//--------------------------------------------------------------------------------------------------
@ -182,16 +170,6 @@ std::pair<int, int> RimFaultReactivationDataAccessorWellLogExtraction::findInter
//--------------------------------------------------------------------------------------------------
std::pair<int, int> RimFaultReactivationDataAccessorWellLogExtraction::findOverburdenAndUnderburdenIndex( const std::vector<double>& values )
{
auto findLastOverburdenIndex = []( const std::vector<double>& values )
{
for ( size_t i = 0; i < values.size(); i++ )
{
if ( !std::isinf( values[i] ) ) return static_cast<int>( i );
}
return -1;
};
auto findFirstUnderburdenIndex = []( const std::vector<double>& values )
{
for ( size_t i = values.size() - 1; i > 0; i-- )
@ -207,6 +185,19 @@ std::pair<int, int> RimFaultReactivationDataAccessorWellLogExtraction::findOverb
return { lastOverburdenIndex, firstUnderburdenIndex };
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
int RimFaultReactivationDataAccessorWellLogExtraction::findLastOverburdenIndex( const std::vector<double>& values )
{
for ( size_t i = 0; i < values.size(); i++ )
{
if ( !std::isinf( values[i] ) ) return static_cast<int>( i );
}
return -1;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
@ -481,3 +472,12 @@ double RimFaultReactivationDataAccessorWellLogExtraction::calculatePorePressure(
{
return RiaEclipseUnitTools::pascalToBar( gradient * 9.81 * depth * 1000.0 );
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double RimFaultReactivationDataAccessorWellLogExtraction::calculateTemperature( double topValue, double topDepth, double depth, double gradient )
{
double tvdDiff = topDepth - depth;
return tvdDiff * gradient + topValue;
}

View File

@ -51,9 +51,9 @@ public:
double gradient );
static std::pair<double, cvf::Vec3d> calculateTemperature( const std::vector<cvf::Vec3d>& intersections,
std::vector<double>& values,
const cvf::Vec3d& position,
double seabedTemperature );
double seabedTemperature,
double gradient );
static std::pair<std::map<RimFaultReactivation::GridPart, cvf::ref<RigWellPath>>,
std::map<RimFaultReactivation::GridPart, cvf::ref<RigEclipseWellLogExtractor>>>
createEclipseWellPathExtractors( const RigFaultReactivationModel& model, RigEclipseCaseData& eclipseCaseData, double seabedDepth );
@ -76,6 +76,10 @@ public:
const cvf::Vec3d& point,
const std::map<RimFaultReactivation::ElementSets, std::vector<unsigned int>>& elementSets );
static int findLastOverburdenIndex( const std::vector<double>& values );
static double computeGradient( double depth1, double value1, double depth2, double value2 );
protected:
static std::pair<int, int> findOverburdenAndUnderburdenIndex( const std::vector<double>& values );
static double computeValueWithGradient( const std::vector<cvf::Vec3d>& intersections,
@ -89,7 +93,6 @@ protected:
static std::pair<double, cvf::Vec3d>
findValueAndPosition( const std::vector<cvf::Vec3d>& intersections, const std::vector<double>& values, const cvf::Vec3d& position );
static double computeGradient( double depth1, double value1, double depth2, double value2 );
static std::vector<double> extractDepthValues( const std::vector<cvf::Vec3d>& intersections );
static void insertUnderburdenValues( const std::vector<cvf::Vec3d>& intersections,
@ -103,4 +106,5 @@ protected:
static double computeMinimumDistance( const cvf::Vec3d& position, const std::vector<cvf::Vec3d>& positions );
static double calculatePorePressure( double depth, double gradient );
static double calculateTemperature( double topValue, double topDepth, double depth, double gradient );
};