Skip to content

Commit

Permalink
#10826 Fault Reactivation: improve pore pressure eclipse extraction.
Browse files Browse the repository at this point in the history
  • Loading branch information
@kriben
kriben committed Jan 3, 2024
1 parent 544e697 commit a64737a
Show file tree
Hide file tree
Showing 9 changed files with 386 additions and 215 deletions.
2 changes: 2 additions & 0 deletions ApplicationLibCode/ProjectDataModel/Faults/CMakeLists_files.cmake
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,7 @@ set(SOURCE_GROUP_HEADER_FILES
${CMAKE_CURRENT_LIST_DIR}/RimFaultReactivationDataAccessorTemperature.h
${CMAKE_CURRENT_LIST_DIR}/RimFaultReactivationDataAccessorGeoMech.h
${CMAKE_CURRENT_LIST_DIR}/RimFaultReactivationDataAccessorStress.h
${CMAKE_CURRENT_LIST_DIR}/RimFaultReactivationDataAccessorWellLogExtraction.h
${CMAKE_CURRENT_LIST_DIR}/RimFaultReactivationEnums.h
)

Expand All @@ -27,6 +28,7 @@ set(SOURCE_GROUP_SOURCE_FILES
${CMAKE_CURRENT_LIST_DIR}/RimFaultReactivationDataAccessorTemperature.cpp
${CMAKE_CURRENT_LIST_DIR}/RimFaultReactivationDataAccessorGeoMech.cpp
${CMAKE_CURRENT_LIST_DIR}/RimFaultReactivationDataAccessorStress.cpp
${CMAKE_CURRENT_LIST_DIR}/RimFaultReactivationDataAccessorWellLogExtraction.cpp
)

list(APPEND CODE_HEADER_FILES ${SOURCE_GROUP_HEADER_FILES})
Expand Down
62 changes: 42 additions & 20 deletions ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorPorePressure.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,21 +17,26 @@
/////////////////////////////////////////////////////////////////////////////////

#include "RimFaultReactivationDataAccessorPorePressure.h"
#include "RimEclipseCase.h"
#include "RimFaciesProperties.h"
#include "RimFaultReactivationDataAccessorWellLogExtraction.h"
#include "RimFaultReactivationEnums.h"

#include "RiaDefines.h"
#include "RiaEclipseUnitTools.h"
#include "RiaPorosityModel.h"

#include "RigCaseCellResultsData.h"
#include "RigEclipseCaseData.h"
#include "RigEclipseResultAddress.h"
#include "RigEclipseWellLogExtractor.h"
#include "RigFaultReactivationModel.h"
#include "RigMainGrid.h"
#include "RigResultAccessorFactory.h"

#include "RimEclipseCase.h"
#include "RigWellPath.h"

#include <cmath>
#include <limits>

//--------------------------------------------------------------------------------------------------
///
Expand Down Expand Up @@ -72,6 +77,24 @@ void RimFaultReactivationDataAccessorPorePressure::updateResultAccessor()
RiaDefines::PorosityModelType::MATRIX_MODEL,
m_timeStep,
resVarAddress );

auto [faultTopPosition, faultBottomPosition] = m_model->faultTopBottom();
auto faultNormal = m_model->faultNormal();
double distanceFromFault = 1.0;

std::string errorName = "fault reactivation data access";

{
std::vector<cvf::Vec3d> wellPoints = generateWellPoints( faultTopPosition, faultBottomPosition, faultNormal * distanceFromFault );
m_faceWellPathA = new RigWellPath( wellPoints, generateMds( wellPoints ) );
m_extractorA = new RigEclipseWellLogExtractor( m_caseData, m_faceWellPathA.p(), errorName );
}

{
std::vector<cvf::Vec3d> wellPoints = generateWellPoints( faultTopPosition, faultBottomPosition, -faultNormal * distanceFromFault );
m_faceWellPathB = new RigWellPath( wellPoints, generateMds( wellPoints ) );
m_extractorB = new RigEclipseWellLogExtractor( m_caseData, m_faceWellPathB.p(), errorName );
}
}
}

Expand All @@ -95,24 +118,23 @@ double RimFaultReactivationDataAccessorPorePressure::valueAtPosition( const cvf:
{
if ( ( m_mainGrid != nullptr ) && m_resultAccessor.notNull() )
{
auto cellIdx = m_mainGrid->findReservoirCellIndexFromPoint( position );
if ( cellIdx != cvf::UNDEFINED_SIZE_T )
{
double value = m_resultAccessor->cellScalar( cellIdx );
if ( !std::isinf( value ) )
{
return 100000.0 * value; // return in pascal, not bar
}
}
}
// TODO: choose correct extractor
auto extractor = m_extractorA;
std::vector<double> values;
extractor->curveData( m_resultAccessor.p(), &values );

return calculatePorePressure( std::abs( position.z() ), m_defaultPorePressureGradient );
}
auto intersections = extractor->intersections();

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double RimFaultReactivationDataAccessorPorePressure::calculatePorePressure( double depth, double gradient )
{
return gradient * 9.81 * depth * 1000.0;
intersections.insert( intersections.begin(), m_faceWellPathA->wellPathPoints().front() );
values.insert( values.begin(), std::numeric_limits<double>::infinity() );

intersections.push_back( m_faceWellPathA->wellPathPoints().back() );
values.push_back( std::numeric_limits<double>::infinity() );

auto [value, pos] =
RimFaultReactivationDataAccessorWellLogExtraction::getPorBar( intersections, values, position, m_defaultPorePressureGradient );
return RiaEclipseUnitTools::barToPascal( value );
}

return std::numeric_limits<double>::infinity();
}
7 changes: 7 additions & 0 deletions ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorPorePressure.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,6 +28,8 @@ class RimEclipseCase;
class RigEclipseCaseData;
class RigMainGrid;
class RigResultAccessor;
class RigEclipseWellLogExtractor;
class RigWellPath;

//==================================================================================================
///
Expand Down Expand Up @@ -58,4 +60,9 @@ class RimFaultReactivationDataAccessorPorePressure : public RimFaultReactivation
const RigMainGrid* m_mainGrid;
double m_defaultPorePressureGradient;
cvf::ref<RigResultAccessor> m_resultAccessor;

cvf::ref<RigWellPath> m_faceWellPathA;
cvf::ref<RigWellPath> m_faceWellPathB;
cvf::ref<RigEclipseWellLogExtractor> m_extractorA;
cvf::ref<RigEclipseWellLogExtractor> m_extractorB;
};
178 changes: 2 additions & 176 deletions ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStress.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -34,6 +34,7 @@
#include "RigResultAccessorFactory.h"
#include "RigWellPath.h"

#include "RimFaultReactivationDataAccessorWellLogExtraction.h"
#include "RimFaultReactivationEnums.h"
#include "RimFracture.h"
#include "RimGeoMechCase.h"
Expand Down Expand Up @@ -235,180 +236,5 @@ std::pair<double, cvf::Vec3d> RimFaultReactivationDataAccessorStress::getPorBar(
std::vector<double> values;
extractor->curveData( resAddr, timeStepIndex, frameIndex, &values );

// Fill in missing values
auto intersections = extractor->intersections();
fillInMissingValues( intersections, values, gradient );

// Linear interpolation between two points
auto lerp = []( const cvf::Vec3d& start, const cvf::Vec3d& end, double t ) { return start + t * ( end - start ); };

auto [topIdx, bottomIdx] = findIntersectionsForTvd( intersections, position.z() );
if ( topIdx != -1 && bottomIdx != -1 )
{
double topValue = values[topIdx];
double bottomValue = values[bottomIdx];
if ( !std::isinf( topValue ) && !std::isinf( bottomValue ) )
{
// Interpolate value from the two closest points.
std::vector<double> xs = { intersections[bottomIdx].z(), intersections[topIdx].z() };
std::vector<double> ys = { values[bottomIdx], values[topIdx] };
double porBar = RiaInterpolationTools::linear( xs, ys, position.z() );

// Interpolate position from depth
double fraction = RiaInterpolationTools::linear( xs, { 0.0, 1.0 }, position.z() );
cvf::Vec3d extractionPosition = lerp( intersections[bottomIdx], intersections[topIdx], fraction );
return { porBar, extractionPosition };
}
}

return { std::numeric_limits<double>::infinity(), cvf::Vec3d::UNDEFINED };
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::pair<int, int> RimFaultReactivationDataAccessorStress::findIntersectionsForTvd( const std::vector<cvf::Vec3d>& intersections, double tvd )
{
int topIdx = -1;
int bottomIdx = -1;
if ( intersections.size() >= 2 )
{
for ( size_t i = 1; i < intersections.size(); i++ )
{
auto top = intersections[i - 1];
auto bottom = intersections[i];
if ( top.z() > tvd && bottom.z() < tvd )
{
topIdx = static_cast<int>( i ) - 1;
bottomIdx = static_cast<int>( i );
break;
}
}
}
return std::make_pair( topIdx, bottomIdx );
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::pair<int, int> RimFaultReactivationDataAccessorStress::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-- )
{
if ( !std::isinf( values[i] ) ) return static_cast<int>( i );
}

return -1;
};

int lastOverburdenIndex = findLastOverburdenIndex( values );
int firstUnderburdenIndex = findFirstUnderburdenIndex( values );
return { lastOverburdenIndex, firstUnderburdenIndex };
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double RimFaultReactivationDataAccessorStress::computePorBarWithGradient( const std::vector<cvf::Vec3d>& intersections,
const std::vector<double>& values,
int i1,
int i2,
double gradient )
{
double tvdDiff = intersections[i2].z() - intersections[i1].z();
return tvdDiff * gradient + values[i2];
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RimFaultReactivationDataAccessorStress::fillInMissingValues( const std::vector<cvf::Vec3d>& intersections,
std::vector<double>& 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 );
}

// Fill in underburden values using gradient
int lastElementIndex = static_cast<int>( 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 );
}

// Interpolate the missing values (should only be intra-reservoir by now)
std::vector<double> intersectionsZ;
for ( auto i : intersections )
{
intersectionsZ.push_back( i.z() );
}

RiaInterpolationTools::interpolateMissingValues( intersectionsZ, values );
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::vector<double> RimFaultReactivationDataAccessorStress::generateMds( const std::vector<cvf::Vec3d>& points )
{
CAF_ASSERT( points.size() >= 2 );

// Assume first at zero, all other points relative to that.
std::vector<double> mds = { 0.0 };
double sum = 0.0;
for ( size_t i = 1; i < points.size(); i++ )
{
sum += points[i - 1].pointDistance( points[i] );
mds.push_back( sum );
}
return mds;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::vector<cvf::Vec3d> RimFaultReactivationDataAccessorStress::generateWellPoints( const cvf::Vec3d& faultTopPosition,
const cvf::Vec3d& faultBottomPosition,
const cvf::Vec3d& offset )
{
cvf::Vec3d faultTop = faultTopPosition + offset;
cvf::Vec3d seabed( faultTop.x(), faultTop.y(), 0.0 );
cvf::Vec3d faultBottom = faultBottomPosition + offset;
cvf::Vec3d underburdenBottom( faultBottom.x(), faultBottom.y(), -10000.0 );
return { seabed, faultTop, faultBottom, underburdenBottom };
return RimFaultReactivationDataAccessorWellLogExtraction::getPorBar( extractor->intersections(), values, position, gradient );
}
17 changes: 0 additions & 17 deletions ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStress.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -78,23 +78,6 @@ class RimFaultReactivationDataAccessorStress : public RimFaultReactivationDataAc
int timeStepIndex,
int frameIndex ) const;

static std::pair<bool, RimFaultReactivation::ElementSets>
findElementSetContainingElement( const std::map<RimFaultReactivation::ElementSets, std::vector<unsigned int>>& elementSets,
unsigned int elmIdx );

static std::pair<int, int> findIntersectionsForTvd( const std::vector<cvf::Vec3d>& intersections, double tvd );
static std::pair<int, int> findOverburdenAndUnderburdenIndex( const std::vector<double>& values );
static double computePorBarWithGradient( const std::vector<cvf::Vec3d>& intersections,
const std::vector<double>& values,
int i1,
int i2,
double gradient );
static void fillInMissingValues( const std::vector<cvf::Vec3d>& intersections, std::vector<double>& values, double gradient );

static std::vector<double> generateMds( const std::vector<cvf::Vec3d>& points );
static std::vector<cvf::Vec3d>
generateWellPoints( const cvf::Vec3d& faultTopPosition, const cvf::Vec3d& faultBottomPosition, const cvf::Vec3d& offset );

RimGeoMechCase* m_geoMechCase;
RimFaultReactivation::Property m_property;
double m_gradient;
Expand Down
Loading

0 comments on commit a64737a

Please sign in to comment.