StimPlan Model Calculator: Add more error checking for stress calculation.

This commit is contained in:
Kristian Bendiksen
2021-05-19 14:22:09 +02:00
committed by Magne Sjaastad
parent 1465b179c5
commit 3b383f1a57
2 changed files with 84 additions and 36 deletions

View File

@@ -364,37 +364,48 @@ bool RimStimPlanModelCalculator::calculateStressWithGradients( std::vector<doubl
int timeStep = m_stimPlanModel->timeStep(); int timeStep = m_stimPlanModel->timeStep();
// Biot coefficient std::vector<RiaDefines::CurveProperty> inputProperties = { RiaDefines::CurveProperty::BIOT_COEFFICIENT,
std::vector<double> biotData = extractValues( RiaDefines::CurveProperty::BIOT_COEFFICIENT, timeStep ); RiaDefines::CurveProperty::K0,
RiaDefines::CurveProperty::PRESSURE,
// K0 RiaDefines::CurveProperty::INITIAL_PRESSURE,
std::vector<double> k0Data = extractValues( RiaDefines::CurveProperty::K0, timeStep ); RiaDefines::CurveProperty::POISSONS_RATIO,
RiaDefines::CurveProperty::PRESSURE_GRADIENT };
// Pressure at the give time step std::map<RiaDefines::CurveProperty, std::vector<double>> inputData;
std::vector<double> timeStepPressureData = extractValues( RiaDefines::CurveProperty::PRESSURE, timeStep ); for ( auto inputProperty : inputProperties )
// Initial pressure
std::vector<double> initialPressureData = extractValues( RiaDefines::CurveProperty::INITIAL_PRESSURE, timeStep );
// Poissons ratio
std::vector<double> poissonsRatioData = extractValues( RiaDefines::CurveProperty::POISSONS_RATIO, timeStep );
// Pressure difference
std::vector<double> pressureDiffData = extractValues( RiaDefines::CurveProperty::PRESSURE_GRADIENT, timeStep );
// Check that we have data from all curves
if ( biotData.empty() || k0Data.empty() || timeStepPressureData.empty() || initialPressureData.empty() ||
poissonsRatioData.empty() || pressureDiffData.empty() )
{ {
return false; QString propertyName = caf::AppEnum<RiaDefines::CurveProperty>::uiText( inputProperty );
} std::vector<double> values = extractValues( inputProperty, timeStep );
// Check that we have data for the curve
if ( values.empty() )
{
RiaLogging::error(
QString( "Unexpected empty input data for %1 in stress calculation for StimPlan Model: %2" )
.arg( propertyName )
.arg( m_stimPlanModel->name() ) );
return false;
}
if ( biotData.size() < layerBoundaryIndexes.size() || k0Data.size() < layerBoundaryIndexes.size() || // Check that there is enough data
timeStepPressureData.size() < layerBoundaryIndexes.size() || if ( values.size() < layerBoundaryIndexes.size() )
initialPressureData.size() < layerBoundaryIndexes.size() || {
poissonsRatioData.size() < layerBoundaryIndexes.size() || pressureDiffData.size() < layerBoundaryIndexes.size() ) RiaLogging::error(
{ QString( "Unexpected input data size for %1 in stress calculation for StimPlan Model: %2" )
return false; .arg( propertyName )
.arg( m_stimPlanModel->name() ) );
return false;
}
// Pressure gradient is allowed to have infs.
if ( inputProperty != RiaDefines::CurveProperty::PRESSURE_GRADIENT &&
!isValidInputData( values, propertyName, layerBoundaryIndexes, layerBoundaryDepths ) )
{
RiaLogging::error( QString( "Unexpected bad data for %1 in stress calculation for StimPlan Model: %2." )
.arg( propertyName )
.arg( m_stimPlanModel->name() ) );
return false;
}
inputData[inputProperty] = values;
} }
// Calculate the stress // Calculate the stress
@@ -403,11 +414,15 @@ bool RimStimPlanModelCalculator::calculateStressWithGradients( std::vector<doubl
double depthTopOfZone = layerBoundaryDepths[i].first; double depthTopOfZone = layerBoundaryDepths[i].first;
// 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( inputData[RiaDefines::CurveProperty::K0], layerBoundaryIndexes, i );
double biot = findValueAtTopOfLayer( biotData, layerBoundaryIndexes, i ); double biot =
double poissonsRatio = findValueAtTopOfLayer( poissonsRatioData, layerBoundaryIndexes, i ); findValueAtTopOfLayer( inputData[RiaDefines::CurveProperty::BIOT_COEFFICIENT], layerBoundaryIndexes, i );
double initialPressure = findValueAtTopOfLayer( initialPressureData, layerBoundaryIndexes, i ); double poissonsRatio =
double timeStepPressure = findValueAtTopOfLayer( timeStepPressureData, layerBoundaryIndexes, i ); findValueAtTopOfLayer( inputData[RiaDefines::CurveProperty::POISSONS_RATIO], layerBoundaryIndexes, i );
double initialPressure =
findValueAtTopOfLayer( inputData[RiaDefines::CurveProperty::INITIAL_PRESSURE], layerBoundaryIndexes, i );
double timeStepPressure =
findValueAtTopOfLayer( inputData[RiaDefines::CurveProperty::PRESSURE], layerBoundaryIndexes, i );
// Vertical stress // Vertical stress
// Use difference between reference depth and depth of top of zone // Use difference between reference depth and depth of top of zone
@@ -434,8 +449,8 @@ bool RimStimPlanModelCalculator::calculateStressWithGradients( std::vector<doubl
depthTopOfZone, depthTopOfZone,
depthBottomOfZone, depthBottomOfZone,
Sv, Sv,
initialPressureData, inputData[RiaDefines::CurveProperty::INITIAL_PRESSURE],
pressureDiffData, inputData[RiaDefines::CurveProperty::PRESSURE_GRADIENT],
stressDepthRef, stressDepthRef,
verticalStressRef, verticalStressRef,
verticalStressGradientRef, verticalStressGradientRef,
@@ -686,3 +701,31 @@ double RimStimPlanModelCalculator::calculateStressAtDepth( double depth,
double stress = verticalStressRef + verticalStressGradientRef * depthDiff; double stress = verticalStressRef + verticalStressGradientRef * depthDiff;
return stress; return stress;
} }
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RimStimPlanModelCalculator::isValidInputData( const std::vector<double>& values,
const QString& propertyName,
const std::vector<std::pair<size_t, size_t>>& layerBoundaryIndexes,
const std::vector<std::pair<double, double>>& layerBoundaryDepths
)
{
for ( size_t i = 0; i < layerBoundaryIndexes.size(); i++ )
{
double value = findValueAtTopOfLayer( values, layerBoundaryIndexes, i );
if ( std::isinf( value ) || std::isnan( value ) )
{
double depthTopOfZone = layerBoundaryDepths[i].first;
RiaLogging::error( QString( "Found invalid value for property %1 top of zone %2, depth %3: %4." )
.arg( propertyName )
.arg( i )
.arg( depthTopOfZone )
.arg( value ) );
return false;
}
}
return true;
}

View File

@@ -116,6 +116,11 @@ protected:
double verticalStressRef, double verticalStressRef,
double verticalStressGradientRef ); double verticalStressGradientRef );
static bool isValidInputData( const std::vector<double>& values,
const QString& propertyName,
const std::vector<std::pair<size_t, size_t>>& layerBoundaryIndexes,
const std::vector<std::pair<double, double>>& layerBoundaryDepths );
private: private:
RimStimPlanModel* m_stimPlanModel; RimStimPlanModel* m_stimPlanModel;
std::vector<std::unique_ptr<RimStimPlanModelPropertyCalculator>> m_resultCalculators; std::vector<std::unique_ptr<RimStimPlanModelPropertyCalculator>> m_resultCalculators;