OilfieldToolsNet/ResInsight
4
0
Fork 0
mirror of https://github.com/OPM/ResInsight.git synced 2025-02-25 18:55:39 -06:00
Code Issues Packages Projects Releases Wiki Activity

#7452 Refactor: Extract method for stress gradient calculation.

Browse Source
This commit is contained in:
Kristian Bendiksen
2021-03-12 15:03:46 +01:00
parent c83af62aae
commit efd3aaeb4e
2 changed files with 84 additions and 42 deletions
Download Patch File Download Diff File Expand all files Collapse all files

114
ApplicationLibCode/ProjectDataModel/StimPlanModel/RimStimPlanModelCalculator.cpp
View File

@@ -400,8 +400,7 @@ bool RimStimPlanModelCalculator::calculateStressWithGradients( std::vector<doubl
// Calculate the stress // Calculate the stress
for ( size_t i = 0; i < layerBoundaryDepths.size(); i++ ) for ( size_t i = 0; i < layerBoundaryDepths.size(); i++ )
{ {
double depthTopOfZone = layerBoundaryDepths[i].first; double depthTopOfZone = layerBoundaryDepths[i].first;
double depthBottomOfZone = layerBoundaryDepths[i].second;
// Data from curves at the top zone depth // Data from curves at the top zone depth
double k0 = findValueAtTopOfLayer( k0Data, layerBoundaryIndexes, i ); double k0 = findValueAtTopOfLayer( k0Data, layerBoundaryIndexes, i );
@@ -428,50 +427,81 @@ bool RimStimPlanModelCalculator::calculateStressWithGradients( std::vector<doubl
initialStress.push_back( RiaEclipseUnitTools::barToPsi( Sh_init ) ); initialStress.push_back( RiaEclipseUnitTools::barToPsi( Sh_init ) );
// Use the bottom of the last layer to compute gradient for last layer // Calculate the stress gradient inside each layer
double bottomInitialPressure = findValueAtBottomOfLayer( initialPressureData, layerBoundaryIndexes, i ); double depthBottomOfZone = layerBoundaryDepths[i].second;
double stressGradient = calculateStressGradientForLayer( i,
double bottomDepthDiff = depthBottomOfZone - stressDepthRef; layerBoundaryIndexes,
double bottomSv = verticalStressRef + verticalStressGradientRef * bottomDepthDiff; depthTopOfZone,
depthBottomOfZone,
double lengthOfLayer = depthBottomOfZone - depthTopOfZone; Sv,
double diffStressLayer = bottomSv - Sv; initialPressureData,
double diffPressureLayer = bottomInitialPressure - initialPressure; pressureDiffData,
stressDepthRef,
// The pressure difference result is only defined where there was no pressure data. verticalStressRef,
// Diff pressure is interpolated in the equilibration region. verticalStressGradientRef,
double diffPressureEQLTop = findValueAtTopOfLayer( pressureDiffData, layerBoundaryIndexes, i ); k0 );
double diffPressureEQLBottom = findValueAtBottomOfLayer( pressureDiffData, layerBoundaryIndexes, i );
if ( !std::isinf( diffPressureEQLTop ) )
{
double offset = RimStimPlanModelPressureCalculator::pressureDifferenceInterpolationOffset();
lengthOfLayer = offset * 2.0;
diffPressureLayer = diffPressureEQLTop;
diffStressLayer = calculateStressDifferenceAtDepth( depthTopOfZone,
offset,
stressDepthRef,
verticalStressRef,
verticalStressGradientRef );
}
else if ( !std::isinf( diffPressureEQLBottom ) )
{
double offset = RimStimPlanModelPressureCalculator::pressureDifferenceInterpolationOffset();
lengthOfLayer = offset * 2.0;
diffPressureLayer = diffPressureEQLBottom;
diffStressLayer = calculateStressDifferenceAtDepth( depthBottomOfZone,
offset,
stressDepthRef,
verticalStressRef,
verticalStressGradientRef );
}
double stressGradient = ( diffStressLayer * k0 + diffPressureLayer * ( 1.0 - k0 ) ) / lengthOfLayer;
stressGradients.push_back( RiaEclipseUnitTools::barPerMeterToPsiPerFeet( stressGradient ) ); stressGradients.push_back( RiaEclipseUnitTools::barPerMeterToPsiPerFeet( stressGradient ) );
} }
return true; return true;
} }
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double RimStimPlanModelCalculator::calculateStressGradientForLayer( size_t i,
std::vector<std::pair<size_t, size_t>> layerBoundaryIndexes,
double depthTopOfZone,
double depthBottomOfZone,
double topSv,
const std::vector<double>& initialPressureData,
const std::vector<double>& pressureDiffData,
double stressDepthRef,
double verticalStressRef,
double verticalStressGradientRef,
double k0 ) const
{
double bottomInitialPressure = findValueAtBottomOfLayer( initialPressureData, layerBoundaryIndexes, i );
double topInitialPressure = findValueAtBottomOfLayer( initialPressureData, layerBoundaryIndexes, i );
double bottomDepthDiff = depthBottomOfZone - stressDepthRef;
double bottomSv = verticalStressRef + verticalStressGradientRef * bottomDepthDiff;
double lengthOfLayer = depthBottomOfZone - depthTopOfZone;
double diffStressLayer = bottomSv - topSv;
double diffPressureLayer = bottomInitialPressure - topInitialPressure;
// The pressure difference result is only defined where there was no pressure data.
// Diff pressure is interpolated in the equilibration region.
double diffPressureEQLTop = findValueAtTopOfLayer( pressureDiffData, layerBoundaryIndexes, i );
double diffPressureEQLBottom = findValueAtBottomOfLayer( pressureDiffData, layerBoundaryIndexes, i );
if ( !std::isinf( diffPressureEQLTop ) )
{
double offset = RimStimPlanModelPressureCalculator::pressureDifferenceInterpolationOffset();
lengthOfLayer = offset * 2.0;
diffPressureLayer = diffPressureEQLTop;
diffStressLayer = calculateStressDifferenceAtDepth( depthTopOfZone,
offset,
stressDepthRef,
verticalStressRef,
verticalStressGradientRef );
}
else if ( !std::isinf( diffPressureEQLBottom ) )
{
double offset = RimStimPlanModelPressureCalculator::pressureDifferenceInterpolationOffset();
lengthOfLayer = offset * 2.0;
diffPressureLayer = diffPressureEQLBottom;
diffStressLayer = calculateStressDifferenceAtDepth( depthBottomOfZone,
offset,
stressDepthRef,
verticalStressRef,
verticalStressGradientRef );
}
double stressGradient = ( diffStressLayer * k0 + diffPressureLayer * ( 1.0 - k0 ) ) / lengthOfLayer;
return stressGradient;
}
//-------------------------------------------------------------------------------------------------- //--------------------------------------------------------------------------------------------------
/// ///
//-------------------------------------------------------------------------------------------------- //--------------------------------------------------------------------------------------------------
@@ -481,8 +511,8 @@ double RimStimPlanModelCalculator::calculateStressDifferenceAtDepth( double dept
double verticalStressRef, double verticalStressRef,
double verticalStressGradientRef ) double verticalStressGradientRef )
{ {
return calculateStressAtDepth( depth - offset, stressDepthRef, verticalStressRef, verticalStressGradientRef ) - return calculateStressAtDepth( depth + offset, stressDepthRef, verticalStressRef, verticalStressGradientRef ) -
calculateStressAtDepth( depth + offset, stressDepthRef, verticalStressRef, verticalStressGradientRef ); calculateStressAtDepth( depth - offset, stressDepthRef, verticalStressRef, verticalStressGradientRef );
} }
//-------------------------------------------------------------------------------------------------- //--------------------------------------------------------------------------------------------------

12
ApplicationLibCode/ProjectDataModel/StimPlanModel/RimStimPlanModelCalculator.h
View File

@@ -93,6 +93,18 @@ protected:
const std::vector<double>& inputVector, const std::vector<double>& inputVector,
std::vector<double>& result ); std::vector<double>& result );
double calculateStressGradientForLayer( size_t i,
std::vector<std::pair<size_t, size_t>> layerBoundaryIndexes,
double depthTopOfZone,
double depthBottomOfZone,
double topSv,
const std::vector<double>& initialPressureData,
const std::vector<double>& pressureDiffData,
double stressDepthRef,
double verticalStressRef,
double verticalStressGradientRef,
double k0 ) const;
static double calculateStressDifferenceAtDepth( double depth, static double calculateStressDifferenceAtDepth( double depth,
double offset, double offset,
double stressDepthRef, double stressDepthRef,