Fault reactivation: use values from GeoMech grid for stress.

No longer snapping to the fake well path along the border between the parts.
This commit is contained in:
Kristian Bendiksen
2024-01-31 21:20:08 +01:00
parent 9fea0853fd
commit 89f90ee9a9
4 changed files with 76 additions and 36 deletions

View File

@@ -114,6 +114,16 @@ double RimFaultReactivationDataAccessorPorePressure::valueAtPosition( const cvf:
auto [value, pos] =
RimFaultReactivationDataAccessorWellLogExtraction::calculatePorBar( intersections, values, position, m_defaultPorePressureGradient );
if ( pos.isUndefined() )
{
auto cellIdx = m_mainGrid->findReservoirCellIndexFromPoint( position );
if ( cellIdx != cvf::UNDEFINED_SIZE_T )
{
double valueFromEclipse = m_resultAccessor->cellScalar( cellIdx );
if ( !std::isinf( valueFromEclipse ) ) return RiaEclipseUnitTools::barToPascal( valueFromEclipse );
}
}
return RiaEclipseUnitTools::barToPascal( value );
}

View File

@@ -88,6 +88,7 @@ void RimFaultReactivationDataAccessorStressGeoMech::updateResultAccessor()
m_s33Frames = loadFrameLambda( femParts, getResultAddress( "ST", "S33" ), timeStepIndex );
m_s11Frames = loadFrameLambda( femParts, getResultAddress( "ST", "S11" ), timeStepIndex );
m_s22Frames = loadFrameLambda( femParts, getResultAddress( "ST", "S22" ), timeStepIndex );
m_porBarFrames = loadFrameLambda( femParts, RigFemAddressDefines::nodalPorBarAddress(), timeStepIndex );
auto [faultTopPosition, faultBottomPosition] = m_model->faultTopBottom();
auto faultNormal = m_model->modelNormal() ^ cvf::Vec3d::Z_AXIS;
@@ -96,31 +97,25 @@ void RimFaultReactivationDataAccessorStressGeoMech::updateResultAccessor()
double distanceFromFault = 1.0;
auto [topDepth, bottomDepth] = m_model->depthTopBottom();
RigFemPartCollection* geoMechPartCollection = m_geoMechCaseData->femParts();
std::string errorName = "fault reactivation data access";
for ( auto gridPart : m_model->allGridParts() )
{
double sign = m_model->normalPointsAt() == gridPart ? -1.0 : 1.0;
std::vector<cvf::Vec3d> wellPoints =
RimFaultReactivationDataAccessorWellLogExtraction::generateWellPoints( faultTopPosition,
faultBottomPosition,
m_seabedDepth,
bottomDepth,
faultNormal * distanceFromFault );
m_faceAWellPath = new RigWellPath( wellPoints, RimFaultReactivationDataAccessorWellLogExtraction::generateMds( wellPoints ) );
m_partIndexA = geoMechPartCollection->getPartIndexFromPoint( wellPoints[1] );
m_extractorA = new RigGeoMechWellLogExtractor( m_geoMechCaseData, partIndex, m_faceAWellPath.p(), errorName );
}
sign * faultNormal * distanceFromFault );
{
std::vector<cvf::Vec3d> wellPoints =
RimFaultReactivationDataAccessorWellLogExtraction::generateWellPoints( faultTopPosition,
faultBottomPosition,
m_seabedDepth,
bottomDepth,
-faultNormal * distanceFromFault );
m_faceBWellPath = new RigWellPath( wellPoints, RimFaultReactivationDataAccessorWellLogExtraction::generateMds( wellPoints ) );
m_partIndexB = geoMechPartCollection->getPartIndexFromPoint( wellPoints[1] );
m_extractorB = new RigGeoMechWellLogExtractor( m_geoMechCaseData, partIndex, m_faceBWellPath.p(), errorName );
cvf::ref<RigWellPath> wellPath =
new RigWellPath( wellPoints, RimFaultReactivationDataAccessorWellLogExtraction::generateMds( wellPoints ) );
m_wellPaths[gridPart] = wellPath;
std::string errorName = "fault reactivation data access";
cvf::ref<RigGeoMechWellLogExtractor> extractor =
new RigGeoMechWellLogExtractor( m_geoMechCaseData, partIndex, wellPath.p(), errorName );
m_extractors[gridPart] = extractor;
}
}
@@ -138,7 +133,7 @@ RigFemResultAddress RimFaultReactivationDataAccessorStressGeoMech::getResultAddr
//--------------------------------------------------------------------------------------------------
bool RimFaultReactivationDataAccessorStressGeoMech::isDataAvailable() const
{
return m_s11Frames && m_s22Frames && m_s33Frames && m_femPart;
return m_s11Frames && m_s22Frames && m_s33Frames && m_porBarFrames && m_femPart;
}
//--------------------------------------------------------------------------------------------------
@@ -183,7 +178,7 @@ std::pair<double, cvf::Vec3d> RimFaultReactivationDataAccessorStressGeoMech::cal
{
int timeStepIndex = 0;
int frameIndex = m_s33Frames->frameCount( timeStepIndex ) - 1;
return calculatePorBar( position, m_gradient, timeStepIndex, frameIndex );
return calculatePorBar( position, m_gradient, gridPart, timeStepIndex, frameIndex );
}
//--------------------------------------------------------------------------------------------------
@@ -215,14 +210,15 @@ double RimFaultReactivationDataAccessorStressGeoMech::interpolatedResultValue( R
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::pair<double, cvf::Vec3d> RimFaultReactivationDataAccessorStressGeoMech::calculatePorBar( const cvf::Vec3d& position,
double gradient,
int timeStepIndex,
int frameIndex ) const
std::pair<double, cvf::Vec3d> RimFaultReactivationDataAccessorStressGeoMech::calculatePorBar( const cvf::Vec3d& position,
double gradient,
RimFaultReactivation::GridPart gridPart,
int timeStepIndex,
int frameIndex ) const
{
RigFemPartCollection* partCollection = m_geoMechCaseData->femParts();
cvf::ref<RigGeoMechWellLogExtractor> extractor = m_partIndexA == partCollection->getPartIndexFromPoint( position ) ? m_extractorA
: m_extractorB;
CAF_ASSERT( m_extractors.find( gridPart ) != m_extractors.end() );
auto extractor = m_extractors.find( gridPart )->second;
if ( !extractor->valid() )
{
RiaLogging::error( "Invalid extractor when extracting PorBar" );
@@ -233,5 +229,28 @@ std::pair<double, cvf::Vec3d> RimFaultReactivationDataAccessorStressGeoMech::cal
std::vector<double> values;
extractor->curveData( resAddr, timeStepIndex, frameIndex, &values );
return RimFaultReactivationDataAccessorWellLogExtraction::calculatePorBar( extractor->intersections(), values, position, gradient );
auto [value, extractionPos] =
RimFaultReactivationDataAccessorWellLogExtraction::calculatePorBar( extractor->intersections(), values, position, gradient );
if ( extractionPos.isUndefined() )
{
// If extraction position is not defined the position is not close to the border between the two parts.
// This means it should be safe to use POR-BAR from the model.
const std::vector<float>& frameData = m_porBarFrames->frameData( timeStepIndex, frameIndex );
// Use data from geo mech grid if defined (only position is reservoir).
RimWellIADataAccess iaDataAccess( m_geoMechCase );
double gridValue = iaDataAccess.interpolatedResultValue( m_femPart, frameData, RIG_NODAL, position );
if ( !std::isinf( gridValue ) )
{
return { gridValue, position };
}
// Use calculated value when POR-BAR is inf (outside of reservoir).
return { value, position };
}
else
{
return { value, extractionPos };
}
}

View File

@@ -77,7 +77,11 @@ private:
const cvf::Vec3d& position,
const std::vector<float>& scalarResults ) const;
std::pair<double, cvf::Vec3d> calculatePorBar( const cvf::Vec3d& position, double gradient, int timeStepIndex, int frameIndex ) const;
std::pair<double, cvf::Vec3d> calculatePorBar( const cvf::Vec3d& position,
double gradient,
RimFaultReactivation::GridPart gridPart,
int timeStepIndex,
int frameIndex ) const;
RigFemScalarResultFrames* dataFrames( StressType stressType ) const;
@@ -86,12 +90,9 @@ private:
RigFemScalarResultFrames* m_s11Frames;
RigFemScalarResultFrames* m_s22Frames;
RigFemScalarResultFrames* m_s33Frames;
RigFemScalarResultFrames* m_porBarFrames;
const RigFemPart* m_femPart;
cvf::ref<RigWellPath> m_faceAWellPath;
cvf::ref<RigWellPath> m_faceBWellPath;
int m_partIndexA;
int m_partIndexB;
cvf::ref<RigGeoMechWellLogExtractor> m_extractorA;
cvf::ref<RigGeoMechWellLogExtractor> m_extractorB;
std::map<RimFaultReactivation::GridPart, cvf::ref<RigWellPath>> m_wellPaths;
std::map<RimFaultReactivation::GridPart, cvf::ref<RigGeoMechWellLogExtractor>> m_extractors;
};

View File

@@ -70,7 +70,17 @@ std::pair<double, cvf::Vec3d> RimFaultReactivationDataAccessorWellLogExtraction:
{
// Fill in missing values
fillInMissingValuesWithGradient( intersections, values, gradient );
return findValueAndPosition( intersections, values, position );
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 };
}
}
//--------------------------------------------------------------------------------------------------