Fault Reactivation: Use computed pore pressure for over and underburden.

This commit is contained in:
Kristian Bendiksen
2024-02-05 17:49:01 +01:00
parent f494813210
commit 5037473f5c
6 changed files with 78 additions and 48 deletions

View File

@@ -112,8 +112,12 @@ double RimFaultReactivationDataAccessorPorePressure::valueAtPosition( const cvf:
auto [values, intersections] =
RimFaultReactivationDataAccessorWellLogExtraction::extractValuesAndIntersections( *m_resultAccessor.p(), *extractor.p(), *wellPath );
auto [value, pos] =
RimFaultReactivationDataAccessorWellLogExtraction::calculatePorBar( intersections, values, position, m_defaultPorePressureGradient );
auto [value, pos] = RimFaultReactivationDataAccessorWellLogExtraction::calculatePorBar( model,
gridPart,
intersections,
values,
position,
m_defaultPorePressureGradient );
if ( pos.isUndefined() )
{
auto cellIdx = m_mainGrid->findReservoirCellIndexFromPoint( position );

View File

@@ -199,7 +199,7 @@ std::pair<double, cvf::Vec3d> RimFaultReactivationDataAccessorStressEclipse::cal
RimFaultReactivationDataAccessorWellLogExtraction::extractValuesAndIntersections( *m_resultAccessor.p(), *extractor.p(), *wellPath );
auto [value, extractionPos] =
RimFaultReactivationDataAccessorWellLogExtraction::calculatePorBar( intersections, values, position, m_gradient );
RimFaultReactivationDataAccessorWellLogExtraction::calculatePorBar( *m_model, gridPart, intersections, values, position, m_gradient );
if ( extractionPos.isUndefined() )
{
auto cellIdx = m_mainGrid->findReservoirCellIndexFromPoint( position );
@@ -248,7 +248,8 @@ std::vector<double>
double deltaDepth = previousDepth - currentDepth;
double density = previousDensity;
auto [isOk, elementSet] = findElementSetForPoint( *part, intersections[i], elementSets );
auto [isOk, elementSet] =
RimFaultReactivationDataAccessorWellLogExtraction::findElementSetForPoint( *part, intersections[i], elementSets );
if ( isOk )
{
// Unit: kg/m^3
@@ -265,37 +266,6 @@ std::vector<double>
return values;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::pair<bool, RimFaultReactivation::ElementSets> RimFaultReactivationDataAccessorStressEclipse::findElementSetForPoint(
const RigGriddedPart3d& part,
const cvf::Vec3d& point,
const std::map<RimFaultReactivation::ElementSets, std::vector<unsigned int>>& elementSets )
{
for ( auto [elementSet, elements] : elementSets )
{
for ( unsigned int elementIndex : elements )
{
// Find nodes from element
auto positions = part.elementCorners( elementIndex );
std::array<cvf::Vec3d, 8> coordinates;
for ( size_t i = 0; i < positions.size(); ++i )
{
coordinates[i] = positions[i];
}
if ( RigHexIntersectionTools::isPointInCell( point, coordinates.data() ) )
{
return { true, elementSet };
}
}
}
return { false, RimFaultReactivation::ElementSets::OverBurden };
};
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------

View File

@@ -80,11 +80,6 @@ private:
static void addOverburdenAndUnderburdenPoints( std::vector<cvf::Vec3d>& intersections, const std::vector<cvf::Vec3d>& wellPathPoints );
static std::pair<bool, RimFaultReactivation::ElementSets>
findElementSetForPoint( const RigGriddedPart3d& part,
const cvf::Vec3d& point,
const std::map<RimFaultReactivation::ElementSets, std::vector<unsigned int>>& elementSets );
RimEclipseCase* m_eclipseCase;
RigEclipseCaseData* m_caseData;
const RigMainGrid* m_mainGrid;

View File

@@ -229,8 +229,12 @@ std::pair<double, cvf::Vec3d> RimFaultReactivationDataAccessorStressGeoMech::cal
std::vector<double> values;
extractor->curveData( resAddr, timeStepIndex, frameIndex, &values );
auto [value, extractionPos] =
RimFaultReactivationDataAccessorWellLogExtraction::calculatePorBar( extractor->intersections(), values, position, gradient );
auto [value, extractionPos] = RimFaultReactivationDataAccessorWellLogExtraction::calculatePorBar( *m_model,
gridPart,
extractor->intersections(),
values,
position,
gradient );
if ( extractionPos.isUndefined() )
{

View File

@@ -63,11 +63,26 @@ RimFaultReactivationDataAccessorWellLogExtraction::~RimFaultReactivationDataAcce
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::pair<double, cvf::Vec3d> RimFaultReactivationDataAccessorWellLogExtraction::calculatePorBar( const std::vector<cvf::Vec3d>& intersections,
std::vector<double>& values,
const cvf::Vec3d& position,
double gradient )
std::pair<double, cvf::Vec3d> RimFaultReactivationDataAccessorWellLogExtraction::calculatePorBar( const RigFaultReactivationModel& model,
RimFaultReactivation::GridPart gridPart,
const std::vector<cvf::Vec3d>& intersections,
std::vector<double>& values,
const cvf::Vec3d& position,
double gradient )
{
auto part = model.grid( gridPart );
CAF_ASSERT( part );
auto elementSets = part->elementSets();
auto [isOk, elementSet] = RimFaultReactivationDataAccessorWellLogExtraction::findElementSetForPoint( *part, position, elementSets );
if ( isOk && ( elementSet == RimFaultReactivation::ElementSets::OverBurden || elementSet == RimFaultReactivation::ElementSets::UnderBurden ) )
{
auto calculatePorePressure = []( double depth, double gradient )
{ return RiaEclipseUnitTools::pascalToBar( gradient * 9.81 * depth * 1000.0 ); };
return { calculatePorePressure( std::abs( position.z() ), gradient ), position };
}
// Fill in missing values
fillInMissingValuesWithGradient( intersections, values, gradient );
auto [value, extractionPosition] = findValueAndPosition( intersections, values, position );
@@ -437,3 +452,34 @@ double RimFaultReactivationDataAccessorWellLogExtraction::computeMinimumDistance
return std::sqrt( minDistance );
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::pair<bool, RimFaultReactivation::ElementSets> RimFaultReactivationDataAccessorWellLogExtraction::findElementSetForPoint(
const RigGriddedPart3d& part,
const cvf::Vec3d& point,
const std::map<RimFaultReactivation::ElementSets, std::vector<unsigned int>>& elementSets )
{
for ( auto [elementSet, elements] : elementSets )
{
for ( unsigned int elementIndex : elements )
{
// Find nodes from element
auto positions = part.elementCorners( elementIndex );
std::array<cvf::Vec3d, 8> coordinates;
for ( size_t i = 0; i < positions.size(); ++i )
{
coordinates[i] = positions[i];
}
if ( RigHexIntersectionTools::isPointInCell( point, coordinates.data() ) )
{
return { true, elementSet };
}
}
}
return { false, RimFaultReactivation::ElementSets::OverBurden };
}

View File

@@ -18,6 +18,7 @@
#pragma once
#include "RigGriddedPart3d.h"
#include "RimFaultReactivationDataAccessor.h"
#include "RimFaultReactivationEnums.h"
@@ -29,6 +30,7 @@ class RigWellPath;
class RigEclipseWellLogExtractor;
class RigEclipseCaseData;
class RigResultAccessor;
class RigGriddedPart3d;
//==================================================================================================
///
@@ -40,8 +42,12 @@ public:
RimFaultReactivationDataAccessorWellLogExtraction();
~RimFaultReactivationDataAccessorWellLogExtraction();
static std::pair<double, cvf::Vec3d>
calculatePorBar( const std::vector<cvf::Vec3d>& intersections, std::vector<double>& values, const cvf::Vec3d& position, double gradient );
static std::pair<double, cvf::Vec3d> calculatePorBar( const RigFaultReactivationModel& model,
RimFaultReactivation::GridPart gridPart,
const std::vector<cvf::Vec3d>& intersections,
std::vector<double>& values,
const cvf::Vec3d& position,
double gradient );
static std::pair<double, cvf::Vec3d> calculateTemperature( const std::vector<cvf::Vec3d>& intersections,
std::vector<double>& values,
@@ -64,6 +70,11 @@ public:
const RigWellPath& wellPath );
static std::pair<int, int> findIntersectionsForTvd( const std::vector<cvf::Vec3d>& intersections, double tvd );
static std::pair<bool, RimFaultReactivation::ElementSets>
findElementSetForPoint( const RigGriddedPart3d& part,
const cvf::Vec3d& point,
const std::map<RimFaultReactivation::ElementSets, std::vector<unsigned int>>& elementSets );
protected:
static std::pair<int, int> findOverburdenAndUnderburdenIndex( const std::vector<double>& values );
static double computeValueWithGradient( const std::vector<cvf::Vec3d>& intersections,