mirror of
https://github.com/OPM/ResInsight.git
synced 2025-01-08 23:23:01 -06:00
446 lines
20 KiB
C++
446 lines
20 KiB
C++
/////////////////////////////////////////////////////////////////////////////////
|
|
//
|
|
// Copyright (C) 2020- Equinor ASA
|
|
//
|
|
// ResInsight is free software: you can redistribute it and/or modify
|
|
// it under the terms of the GNU General Public License as published by
|
|
// the Free Software Foundation, either version 3 of the License, or
|
|
// (at your option) any later version.
|
|
//
|
|
// ResInsight is distributed in the hope that it will be useful, but WITHOUT ANY
|
|
// WARRANTY; without even the implied warranty of MERCHANTABILITY or
|
|
// FITNESS FOR A PARTICULAR PURPOSE.
|
|
//
|
|
// See the GNU General Public License at <http://www.gnu.org/licenses/gpl.html>
|
|
// for more details.
|
|
//
|
|
/////////////////////////////////////////////////////////////////////////////////
|
|
#include "RimStimPlanModelWellLogCalculator.h"
|
|
|
|
#include "RiaDefines.h"
|
|
#include "RiaInterpolationTools.h"
|
|
#include "RiaLogging.h"
|
|
#include "RiaStimPlanModelDefines.h"
|
|
|
|
#include "RigEclipseCaseData.h"
|
|
#include "RigEclipseWellLogExtractor.h"
|
|
#include "RigResultAccessor.h"
|
|
#include "RigResultAccessorFactory.h"
|
|
#include "RigWellLogCurveData.h"
|
|
#include "RigWellPath.h"
|
|
|
|
#include "RimCase.h"
|
|
#include "RimEclipseCase.h"
|
|
#include "RimEclipseInputProperty.h"
|
|
#include "RimEclipseInputPropertyCollection.h"
|
|
#include "RimEclipseResultDefinition.h"
|
|
#include "RimModeledWellPath.h"
|
|
#include "RimNonNetLayers.h"
|
|
#include "RimStimPlanModel.h"
|
|
#include "RimStimPlanModelCalculator.h"
|
|
#include "RimStimPlanModelTemplate.h"
|
|
#include "RimWellPath.h"
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
RimStimPlanModelWellLogCalculator::RimStimPlanModelWellLogCalculator( RimStimPlanModelCalculator* stimPlanModelCalculator )
|
|
: m_stimPlanModelCalculator( stimPlanModelCalculator )
|
|
{
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
bool RimStimPlanModelWellLogCalculator::isMatching( RiaDefines::CurveProperty curveProperty ) const
|
|
{
|
|
std::vector<RiaDefines::CurveProperty> matching = {
|
|
RiaDefines::CurveProperty::FACIES,
|
|
RiaDefines::CurveProperty::POROSITY,
|
|
RiaDefines::CurveProperty::POROSITY_UNSCALED,
|
|
RiaDefines::CurveProperty::PERMEABILITY_X,
|
|
RiaDefines::CurveProperty::PERMEABILITY_Z,
|
|
RiaDefines::CurveProperty::INITIAL_PRESSURE,
|
|
RiaDefines::CurveProperty::PRESSURE,
|
|
RiaDefines::CurveProperty::NET_TO_GROSS,
|
|
};
|
|
|
|
return std::find( matching.begin(), matching.end(), curveProperty ) != matching.end();
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
bool RimStimPlanModelWellLogCalculator::calculate( RiaDefines::CurveProperty curveProperty,
|
|
const RimStimPlanModel* stimPlanModel,
|
|
int timeStep,
|
|
std::vector<double>& values,
|
|
std::vector<double>& measuredDepthValues,
|
|
std::vector<double>& tvDepthValues,
|
|
double& rkbDiff ) const
|
|
{
|
|
RimEclipseCase* eclipseCase = stimPlanModel->eclipseCaseForProperty( curveProperty );
|
|
if ( !eclipseCase )
|
|
{
|
|
return false;
|
|
}
|
|
|
|
if ( !stimPlanModel->thicknessDirectionWellPath() )
|
|
{
|
|
return false;
|
|
}
|
|
|
|
RigWellPath* wellPathGeometry = stimPlanModel->thicknessDirectionWellPath()->wellPathGeometry();
|
|
if ( !wellPathGeometry )
|
|
{
|
|
RiaLogging::error( "No well path geometry found for well log exctration" );
|
|
return false;
|
|
}
|
|
|
|
RigEclipseWellLogExtractor eclExtractor( eclipseCase->eclipseCaseData(), wellPathGeometry, "fracture model" );
|
|
|
|
measuredDepthValues = eclExtractor.cellIntersectionMDs();
|
|
tvDepthValues = eclExtractor.cellIntersectionTVDs();
|
|
rkbDiff = eclExtractor.wellPathGeometry()->rkbDiff();
|
|
|
|
RimEclipseResultDefinition eclipseResultDefinition;
|
|
eclipseResultDefinition.setEclipseCase( eclipseCase );
|
|
eclipseResultDefinition.setResultType( stimPlanModel->eclipseResultCategory( curveProperty ) );
|
|
eclipseResultDefinition.setPorosityModel( RiaDefines::PorosityModelType::MATRIX_MODEL );
|
|
eclipseResultDefinition.setResultVariable( stimPlanModel->eclipseResultVariable( curveProperty ) );
|
|
|
|
eclipseResultDefinition.loadResult();
|
|
|
|
if ( stimPlanModel->eclipseResultCategory( curveProperty ) != RiaDefines::ResultCatType::DYNAMIC_NATIVE ||
|
|
curveProperty == RiaDefines::CurveProperty::INITIAL_PRESSURE )
|
|
{
|
|
timeStep = 0;
|
|
}
|
|
|
|
cvf::ref<RigResultAccessor> resAcc =
|
|
RigResultAccessorFactory::createFromResultDefinition( eclipseCase->eclipseCaseData(),
|
|
0,
|
|
timeStep,
|
|
&eclipseResultDefinition );
|
|
|
|
if ( resAcc.notNull() )
|
|
{
|
|
eclExtractor.curveData( resAcc.p(), &values );
|
|
}
|
|
else
|
|
{
|
|
RiaLogging::error( QString( "No result found for %1" ).arg( eclipseResultDefinition.resultVariable() ) );
|
|
return false;
|
|
}
|
|
|
|
double overburdenHeight = stimPlanModel->overburdenHeight();
|
|
if ( overburdenHeight > 0.0 )
|
|
{
|
|
addOverburden( curveProperty, stimPlanModel, tvDepthValues, measuredDepthValues, values );
|
|
}
|
|
|
|
double underburdenHeight = stimPlanModel->underburdenHeight();
|
|
if ( underburdenHeight > 0.0 )
|
|
{
|
|
addUnderburden( curveProperty, stimPlanModel, tvDepthValues, measuredDepthValues, values );
|
|
}
|
|
|
|
if ( hasMissingValues( values ) )
|
|
{
|
|
if ( stimPlanModel->missingValueStrategy( curveProperty ) == RimStimPlanModel::MissingValueStrategy::DEFAULT_VALUE )
|
|
{
|
|
// Try to locate a backup accessor (e.g. PORO_1 for PORO)
|
|
cvf::ref<RigResultAccessor> backupResAcc = findMissingValuesAccessor( eclipseCase->eclipseCaseData(),
|
|
eclipseCase->inputPropertyCollection(),
|
|
0,
|
|
timeStep,
|
|
&eclipseResultDefinition );
|
|
|
|
if ( backupResAcc.notNull() )
|
|
{
|
|
RiaLogging::info( QString( "Reading missing values from input properties for %1." )
|
|
.arg( eclipseResultDefinition.resultVariable() ) );
|
|
std::vector<double> replacementValues;
|
|
eclExtractor.curveData( backupResAcc.p(), &replacementValues );
|
|
|
|
double overburdenHeight = stimPlanModel->overburdenHeight();
|
|
if ( overburdenHeight > 0.0 )
|
|
{
|
|
double defaultOverburdenValue = std::numeric_limits<double>::infinity();
|
|
if ( stimPlanModel->burdenStrategy( curveProperty ) == RimStimPlanModel::BurdenStrategy::DEFAULT_VALUE )
|
|
{
|
|
defaultOverburdenValue = stimPlanModel->getDefaultForMissingOverburdenValue( curveProperty );
|
|
}
|
|
|
|
replacementValues.insert( replacementValues.begin(), defaultOverburdenValue );
|
|
replacementValues.insert( replacementValues.begin(), defaultOverburdenValue );
|
|
}
|
|
|
|
double underburdenHeight = stimPlanModel->underburdenHeight();
|
|
if ( underburdenHeight > 0.0 )
|
|
{
|
|
double defaultUnderburdenValue = std::numeric_limits<double>::infinity();
|
|
if ( stimPlanModel->burdenStrategy( curveProperty ) == RimStimPlanModel::BurdenStrategy::DEFAULT_VALUE )
|
|
{
|
|
defaultUnderburdenValue = stimPlanModel->getDefaultForMissingUnderburdenValue( curveProperty );
|
|
}
|
|
|
|
replacementValues.push_back( defaultUnderburdenValue );
|
|
replacementValues.push_back( defaultUnderburdenValue );
|
|
}
|
|
|
|
replaceMissingValues( values, replacementValues );
|
|
}
|
|
|
|
// If the backup accessor is not found, or does not provide all the missing values:
|
|
// use default value from the fracture model
|
|
if ( !backupResAcc.notNull() || hasMissingValues( values ) )
|
|
{
|
|
double defaultValue = stimPlanModel->getDefaultForMissingValue( curveProperty );
|
|
replaceMissingValues( values, defaultValue );
|
|
}
|
|
}
|
|
else if ( stimPlanModel->missingValueStrategy( curveProperty ) ==
|
|
RimStimPlanModel::MissingValueStrategy::LINEAR_INTERPOLATION )
|
|
{
|
|
RiaLogging::info(
|
|
QString( "Interpolating missing values for %1" ).arg( eclipseResultDefinition.resultVariable() ) );
|
|
RiaInterpolationTools::interpolateMissingValues( measuredDepthValues, values );
|
|
}
|
|
else
|
|
{
|
|
// Get the missing data from other curve
|
|
RiaDefines::CurveProperty replacementProperty =
|
|
stimPlanModel->getDefaultPropertyForMissingValues( curveProperty );
|
|
|
|
std::vector<double> initialValues;
|
|
std::vector<double> initialMds;
|
|
std::vector<double> initialTvds;
|
|
double initialRkbDiff = -1.0;
|
|
calculate( replacementProperty, stimPlanModel, timeStep, initialValues, initialMds, initialTvds, initialRkbDiff );
|
|
|
|
if ( initialValues.empty() )
|
|
{
|
|
RiaLogging::error( QString( "Empty replacement data found for fracture model curve." ) );
|
|
return false;
|
|
}
|
|
|
|
CVF_ASSERT( values.size() == initialValues.size() );
|
|
replaceMissingValues( values, initialValues );
|
|
}
|
|
}
|
|
|
|
if ( stimPlanModel->isScaledByNetToGross( curveProperty ) )
|
|
{
|
|
std::vector<double> netToGross =
|
|
m_stimPlanModelCalculator->extractValues( RiaDefines::CurveProperty::NET_TO_GROSS, timeStep );
|
|
|
|
scaleByNetToGross( stimPlanModel, netToGross, values );
|
|
}
|
|
|
|
return true;
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
bool RimStimPlanModelWellLogCalculator::hasMissingValues( const std::vector<double>& values )
|
|
{
|
|
for ( double v : values )
|
|
{
|
|
if ( v == std::numeric_limits<double>::infinity() )
|
|
{
|
|
return true;
|
|
}
|
|
}
|
|
|
|
return false;
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
void RimStimPlanModelWellLogCalculator::replaceMissingValues( std::vector<double>& values, double defaultValue )
|
|
{
|
|
for ( double& v : values )
|
|
{
|
|
if ( v == std::numeric_limits<double>::infinity() )
|
|
{
|
|
v = defaultValue;
|
|
}
|
|
}
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
void RimStimPlanModelWellLogCalculator::replaceMissingValues( std::vector<double>& values,
|
|
const std::vector<double>& replacementValues )
|
|
{
|
|
CVF_ASSERT( values.size() == replacementValues.size() );
|
|
for ( size_t i = 0; i < values.size(); i++ )
|
|
{
|
|
if ( values[i] == std::numeric_limits<double>::infinity() )
|
|
{
|
|
values[i] = replacementValues[i];
|
|
}
|
|
}
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
cvf::ref<RigResultAccessor>
|
|
RimStimPlanModelWellLogCalculator::findMissingValuesAccessor( RigEclipseCaseData* caseData,
|
|
RimEclipseInputPropertyCollection* inputPropertyCollection,
|
|
int gridIndex,
|
|
int timeStepIndex,
|
|
RimEclipseResultDefinition* eclipseResultDefinition ) const
|
|
{
|
|
QString resultName = eclipseResultDefinition->resultVariable();
|
|
|
|
for ( RimEclipseInputProperty* inputProperty : inputPropertyCollection->inputProperties() )
|
|
{
|
|
// Look for input properties starting with the same name as result definition
|
|
if ( inputProperty && inputProperty->resultName().startsWith( resultName ) )
|
|
{
|
|
RiaLogging::info(
|
|
QString( "Found missing values result for %1: %2" ).arg( resultName ).arg( inputProperty->resultName() ) );
|
|
|
|
RigEclipseResultAddress resultAddress( RiaDefines::ResultCatType::INPUT_PROPERTY, inputProperty->resultName() );
|
|
caseData->results( eclipseResultDefinition->porosityModel() )->ensureKnownResultLoaded( resultAddress );
|
|
cvf::ref<RigResultAccessor> resAcc =
|
|
RigResultAccessorFactory::createFromResultAddress( caseData,
|
|
gridIndex,
|
|
eclipseResultDefinition->porosityModel(),
|
|
timeStepIndex,
|
|
resultAddress );
|
|
|
|
return resAcc;
|
|
}
|
|
}
|
|
|
|
return nullptr;
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
void RimStimPlanModelWellLogCalculator::addOverburden( RiaDefines::CurveProperty curveProperty,
|
|
const RimStimPlanModel* stimPlanModel,
|
|
std::vector<double>& tvDepthValues,
|
|
std::vector<double>& measuredDepthValues,
|
|
std::vector<double>& values ) const
|
|
{
|
|
if ( !values.empty() )
|
|
{
|
|
double overburdenHeight = stimPlanModel->overburdenHeight();
|
|
double tvdOverburdenBottom = tvDepthValues[0];
|
|
double tvdOverburdenTop = tvdOverburdenBottom - overburdenHeight;
|
|
|
|
double overburdenTopValue = std::numeric_limits<double>::infinity();
|
|
double overburdenBottomValue = std::numeric_limits<double>::infinity();
|
|
if ( stimPlanModel->burdenStrategy( curveProperty ) == RimStimPlanModel::BurdenStrategy::DEFAULT_VALUE )
|
|
{
|
|
overburdenTopValue = stimPlanModel->getDefaultForMissingOverburdenValue( curveProperty );
|
|
overburdenBottomValue = overburdenTopValue;
|
|
}
|
|
else
|
|
{
|
|
double gradient = stimPlanModel->getOverburdenGradient( curveProperty );
|
|
overburdenBottomValue = values[0];
|
|
overburdenTopValue = overburdenBottomValue + gradient * -overburdenHeight;
|
|
}
|
|
|
|
// Prepend the new "fake" depth for start of overburden
|
|
tvDepthValues.insert( tvDepthValues.begin(), tvdOverburdenBottom );
|
|
tvDepthValues.insert( tvDepthValues.begin(), tvdOverburdenTop );
|
|
|
|
// TODO: this is not always correct
|
|
double mdTop = measuredDepthValues[0];
|
|
measuredDepthValues.insert( measuredDepthValues.begin(), mdTop );
|
|
measuredDepthValues.insert( measuredDepthValues.begin(), mdTop - overburdenHeight );
|
|
|
|
values.insert( values.begin(), overburdenBottomValue );
|
|
values.insert( values.begin(), overburdenTopValue );
|
|
}
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
void RimStimPlanModelWellLogCalculator::addUnderburden( RiaDefines::CurveProperty curveProperty,
|
|
const RimStimPlanModel* stimPlanModel,
|
|
std::vector<double>& tvDepthValues,
|
|
std::vector<double>& measuredDepthValues,
|
|
std::vector<double>& values ) const
|
|
{
|
|
if ( !values.empty() )
|
|
{
|
|
size_t lastIndex = tvDepthValues.size() - 1;
|
|
|
|
double underburdenHeight = stimPlanModel->underburdenHeight();
|
|
double tvdUnderburdenTop = tvDepthValues[lastIndex];
|
|
double tvdUnderburdenBottom = tvdUnderburdenTop + underburdenHeight;
|
|
|
|
double underburdenTopValue = std::numeric_limits<double>::infinity();
|
|
double underburdenBottomValue = std::numeric_limits<double>::infinity();
|
|
if ( stimPlanModel->burdenStrategy( curveProperty ) == RimStimPlanModel::BurdenStrategy::DEFAULT_VALUE )
|
|
{
|
|
underburdenTopValue = stimPlanModel->getDefaultForMissingUnderburdenValue( curveProperty );
|
|
underburdenBottomValue = underburdenTopValue;
|
|
}
|
|
else
|
|
{
|
|
double gradient = stimPlanModel->getUnderburdenGradient( curveProperty );
|
|
underburdenTopValue = values[lastIndex];
|
|
underburdenBottomValue = underburdenTopValue + gradient * underburdenHeight;
|
|
}
|
|
|
|
// Append the new "fake" depth for start of underburden
|
|
tvDepthValues.push_back( tvdUnderburdenTop );
|
|
tvDepthValues.push_back( tvdUnderburdenBottom );
|
|
|
|
// Append the new "fake" md
|
|
// TODO: check if this is correct???
|
|
double mdBottom = measuredDepthValues[lastIndex];
|
|
measuredDepthValues.push_back( mdBottom );
|
|
measuredDepthValues.push_back( mdBottom + underburdenHeight );
|
|
|
|
values.push_back( underburdenTopValue );
|
|
values.push_back( underburdenBottomValue );
|
|
}
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
void RimStimPlanModelWellLogCalculator::scaleByNetToGross( const RimStimPlanModel* stimPlanModel,
|
|
const std::vector<double>& netToGross,
|
|
std::vector<double>& values )
|
|
{
|
|
if ( netToGross.size() != values.size() )
|
|
{
|
|
RiaLogging::error( QString( "Different sizes for net to gross calculation." ) );
|
|
return;
|
|
}
|
|
|
|
double netToGrossCutoff = 1.0;
|
|
bool useNetToGross = false;
|
|
|
|
if ( stimPlanModel->stimPlanModelTemplate() && stimPlanModel->stimPlanModelTemplate()->nonNetLayers() )
|
|
{
|
|
netToGrossCutoff = stimPlanModel->stimPlanModelTemplate()->nonNetLayers()->cutOff();
|
|
useNetToGross = !netToGross.empty() && stimPlanModel->stimPlanModelTemplate()->nonNetLayers()->isChecked();
|
|
}
|
|
|
|
for ( size_t i = 0; i < values.size(); i++ )
|
|
{
|
|
double ntg = netToGross[i];
|
|
if ( useNetToGross && ntg <= netToGrossCutoff )
|
|
{
|
|
values[i] = ntg * values[i];
|
|
}
|
|
}
|
|
}
|