#7400 StimPlanModel: Search vertically for EQLNUM for inactive regions.

This commit is contained in:
Kristian Bendiksen 2021-02-24 19:43:56 +01:00 committed by Magne Sjaastad
parent 8ce620323e
commit 5c5b9e7f34
5 changed files with 158 additions and 4 deletions

View File

@ -1335,6 +1335,10 @@ std::deque<RimStimPlanModel::MissingValueStrategy>
{
if ( curveProperty == RiaDefines::CurveProperty::PRESSURE )
return { RimStimPlanModel::MissingValueStrategy::OTHER_CURVE_PROPERTY };
else if ( curveProperty == RiaDefines::CurveProperty::EQLNUM )
return { RimStimPlanModel::MissingValueStrategy::DEFAULT_VALUE,
RimStimPlanModel::MissingValueStrategy::CELLS_ABOVE,
RimStimPlanModel::MissingValueStrategy::CELLS_BELOW };
else
return { RimStimPlanModel::MissingValueStrategy::DEFAULT_VALUE };
}

View File

@ -71,7 +71,9 @@ public:
{
DEFAULT_VALUE,
LINEAR_INTERPOLATION,
OTHER_CURVE_PROPERTY
OTHER_CURVE_PROPERTY,
CELLS_ABOVE,
CELLS_BELOW
};
enum class BurdenStrategy

View File

@ -516,9 +516,12 @@ bool RimStimPlanModelPressureCalculator::interpolateInitialPressureByEquilibrati
{
int eqlNum = static_cast<int>( eqlNumValues[i] );
DepthValuePairVector depthValuePairs = valuesPerEqlNum[eqlNum];
double depth = tvDepthValues[i];
double value = interpolatePressure( depthValuePairs, depth, eqlNum );
values[i] = value;
if ( !depthValuePairs.empty() )
{
double depth = tvDepthValues[i];
double value = interpolatePressure( depthValuePairs, depth, eqlNum );
values[i] = value;
}
}
}

View File

@ -22,8 +22,10 @@
#include "RiaLogging.h"
#include "RiaStimPlanModelDefines.h"
#include "RigActiveCellInfo.h"
#include "RigEclipseCaseData.h"
#include "RigEclipseWellLogExtractor.h"
#include "RigMainGrid.h"
#include "RigResultAccessor.h"
#include "RigResultAccessorFactory.h"
#include "RigWellLogCurveData.h"
@ -121,6 +123,24 @@ bool RimStimPlanModelWellLogCalculator::calculate( RiaDefines::CurveProperty cur
return false;
}
}
else if ( strategy == RimStimPlanModel::MissingValueStrategy::CELLS_ABOVE )
{
// K-1 is up
int kDirection = -1;
if ( !replaceMissingValuesWithOtherKLayer( curveProperty, stimPlanModel, timeStep, measuredDepthValues, values, kDirection ) )
{
return false;
}
}
else if ( strategy == RimStimPlanModel::MissingValueStrategy::CELLS_BELOW )
{
// K+1 is down
int kDirection = 1;
if ( !replaceMissingValuesWithOtherKLayer( curveProperty, stimPlanModel, timeStep, measuredDepthValues, values, kDirection ) )
{
return false;
}
}
}
if ( stimPlanModel->isScaledByNetToGross( curveProperty ) )
@ -486,6 +506,9 @@ bool RimStimPlanModelWellLogCalculator::replaceMissingValuesWithDefault( RiaDefi
return true;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RimStimPlanModelWellLogCalculator::replaceMissingValuesWithOtherProperty( RiaDefines::CurveProperty curveProperty,
const RimStimPlanModel* stimPlanModel,
int timeStep,
@ -510,3 +533,111 @@ bool RimStimPlanModelWellLogCalculator::replaceMissingValuesWithOtherProperty( R
replaceMissingValues( values, initialValues );
return true;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
const std::vector<double>& RimStimPlanModelWellLogCalculator::loadResults( RigEclipseCaseData* caseData,
RiaDefines::PorosityModelType porosityModel,
RiaDefines::ResultCatType resultType,
const QString& propertyName )
{
// TODO: is this always enough?
auto resultData = caseData->results( porosityModel );
int timeStepIndex = 0;
RigEclipseResultAddress resultAddress( resultType, propertyName );
resultData->ensureKnownResultLoaded( resultAddress );
return caseData->results( porosityModel )->cellScalarResults( resultAddress, timeStepIndex );
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RimStimPlanModelWellLogCalculator::replaceMissingValuesWithOtherKLayer( RiaDefines::CurveProperty curveProperty,
const RimStimPlanModel* stimPlanModel,
int timeStep,
const std::vector<double>& measuredDepths,
std::vector<double>& values,
int moveDirection ) const
{
RimEclipseCase* eclipseCase = stimPlanModel->eclipseCaseForProperty( curveProperty );
if ( !eclipseCase ) return false;
if ( !stimPlanModel->thicknessDirectionWellPath() ) return false;
RigWellPath* wellPathGeometry = stimPlanModel->thicknessDirectionWellPath()->wellPathGeometry();
if ( !wellPathGeometry ) return false;
const RigMainGrid* mainGrid = eclipseCase->mainGrid();
RigEclipseCaseData* caseData = eclipseCase->eclipseCaseData();
RiaDefines::PorosityModelType porosityModel = RiaDefines::PorosityModelType::MATRIX_MODEL;
RiaDefines::ResultCatType resultType = stimPlanModel->eclipseResultCategory( curveProperty );
QString resultName = stimPlanModel->eclipseResultVariable( curveProperty );
const std::vector<double>& cellValues = loadResults( caseData, porosityModel, resultType, resultName );
auto activeCellInfo = caseData->activeCellInfo( porosityModel );
for ( size_t idx = 0; idx < values.size(); idx++ )
{
if ( std::isinf( values[idx] ) )
{
double measuredDepth = measuredDepths[idx];
cvf::Vec3d position = wellPathGeometry->interpolatedPointAlongWellPath( measuredDepth );
size_t cellIdx = mainGrid->findReservoirCellIndexFromPoint( position );
if ( cellIdx != cvf::UNDEFINED_SIZE_T )
{
size_t i;
size_t j;
size_t k;
bool isValid = mainGrid->ijkFromCellIndex( cellIdx, &i, &j, &k );
if ( isValid )
{
RiaLogging::info( QString( "Missing value at MD: %1. Cell [%2, %3, %4]" )
.arg( measuredDepth )
.arg( i + 1 )
.arg( j + 1 )
.arg( k + 1 ) );
int neighborK = k + moveDirection;
const int minK = 1;
const int maxK = mainGrid->cellCountK();
bool isFound = false;
while ( !isFound && neighborK >= minK && neighborK <= maxK )
{
size_t neighborCellIdx = mainGrid->cellIndexFromIJK( i, j, neighborK );
size_t resultIdx = activeCellInfo->cellResultIndex( neighborCellIdx );
if ( neighborCellIdx != cvf::UNDEFINED_SIZE_T && resultIdx < cellValues.size() )
{
double neighborValue = cellValues[resultIdx];
if ( !std::isinf( neighborValue ) )
{
RiaLogging::info( QString( "Found value. [%1, %2, %3]. Idx=%4 Value=%5." )
.arg( i + 1 )
.arg( j + 1 )
.arg( neighborK + 1 )
.arg( neighborCellIdx )
.arg( neighborValue ) );
isFound = true;
values[idx] = neighborValue;
}
}
neighborK += moveDirection;
}
}
}
}
}
return true;
}

View File

@ -20,6 +20,8 @@
#include "RimStimPlanModelCalculator.h"
#include "RimStimPlanModelPropertyCalculator.h"
#include "RiaDefines.h"
#include "RiaPorosityModel.h"
#include "RiaStimPlanModelDefines.h"
#include "cvfObject.h"
@ -46,6 +48,11 @@ public:
bool isMatching( RiaDefines::CurveProperty curveProperty ) const override;
static const std::vector<double>& loadResults( RigEclipseCaseData* caseData,
RiaDefines::PorosityModelType porosityModel,
RiaDefines::ResultCatType resultType,
const QString& propertyName );
protected:
static bool hasMissingValues( const std::vector<double>& values );
static void replaceMissingValues( std::vector<double>& values, double defaultValue );
@ -89,5 +96,12 @@ protected:
int timeStep,
std::vector<double>& values ) const;
bool replaceMissingValuesWithOtherKLayer( RiaDefines::CurveProperty curveProperty,
const RimStimPlanModel* stimPlanModel,
int timeStep,
const std::vector<double>& measuredDepths,
std::vector<double>& values,
int kDirection ) const;
RimStimPlanModelCalculator* m_stimPlanModelCalculator;
};