#5542 Fix scaling issues causing wrong FG results in some cases.

This commit is contained in:
Gaute Lindkvist 2020-02-13 10:22:22 +01:00
parent a35c395f54
commit d5351585f0
3 changed files with 41 additions and 33 deletions

View File

@ -129,7 +129,7 @@ double RigGeoMechBoreHoleStressCalculator::solveSecant( MemberFunc fn, double* t
double f_x0 = ( this->*fn )( x_0, &theta ); double f_x0 = ( this->*fn )( x_0, &theta );
double x_1 = m_porePressure; double x_1 = m_porePressure;
double f_x1 = ( this->*fn )( x_1, &theta ); double f_x1 = ( this->*fn )( x_1, &theta );
double x = 0.0; double x = m_porePressure;
double f_x = 0.0; double f_x = 0.0;
int i = 0; int i = 0;
for ( ; i <= N && std::abs( f_x1 - f_x0 ) > epsilon; ++i ) for ( ; i <= N && std::abs( f_x1 - f_x0 ) > epsilon; ++i )
@ -145,7 +145,7 @@ double RigGeoMechBoreHoleStressCalculator::solveSecant( MemberFunc fn, double* t
f_x1 = f_x; f_x1 = f_x;
} }
if ( i == N || std::abs( f_x ) > epsilon * m_porePressure ) if ( i == 0 || i == N || std::abs( f_x ) > epsilon * m_porePressure )
{ {
// Fallback to bisection if secant doesn't converge or converged to a wrong solution. // Fallback to bisection if secant doesn't converge or converged to a wrong solution.
return solveBisection( 0.0, m_porePressure * 2.0, fn, thetaOut ); return solveBisection( 0.0, m_porePressure * 2.0, fn, thetaOut );

View File

@ -171,7 +171,7 @@ QString RigGeoMechWellLogExtractor::curveData( const RigFemResultAddress& resAdd
{ {
frameIndex = 0; frameIndex = 0;
} }
calculateWbsParameterForAllSegments( param, frameIndex, values ); calculateWbsParameterForAllSegments( param, frameIndex, values, true );
if ( param == RigWbsParameter::UCS() ) // UCS is reported as UCS/100 if ( param == RigWbsParameter::UCS() ) // UCS is reported as UCS/100
{ {
for ( double& value : *values ) for ( double& value : *values )
@ -220,7 +220,8 @@ std::vector<RigGeoMechWellLogExtractor::WbsParameterSource>
RigGeoMechWellLogExtractor::calculateWbsParameterForAllSegments( const RigWbsParameter& parameter, RigGeoMechWellLogExtractor::calculateWbsParameterForAllSegments( const RigWbsParameter& parameter,
WbsParameterSource primarySource, WbsParameterSource primarySource,
int frameIndex, int frameIndex,
std::vector<double>* outputValues ) std::vector<double>* outputValues,
bool allowNormalization )
{ {
RigFemPartResultsCollection* resultCollection = m_caseData->femPartResults(); RigFemPartResultsCollection* resultCollection = m_caseData->femPartResults();
@ -332,7 +333,7 @@ std::vector<RigGeoMechWellLogExtractor::WbsParameterSource>
} }
} }
if ( parameter.normalizeByHydrostaticPP() ) if ( allowNormalization && parameter.normalizeByHydrostaticPP() )
{ {
outputValues->resize( unscaledValues.size(), std::numeric_limits<double>::infinity() ); outputValues->resize( unscaledValues.size(), std::numeric_limits<double>::infinity() );
@ -366,9 +367,14 @@ std::vector<RigGeoMechWellLogExtractor::WbsParameterSource>
std::vector<RigGeoMechWellLogExtractor::WbsParameterSource> std::vector<RigGeoMechWellLogExtractor::WbsParameterSource>
RigGeoMechWellLogExtractor::calculateWbsParameterForAllSegments( const RigWbsParameter& parameter, RigGeoMechWellLogExtractor::calculateWbsParameterForAllSegments( const RigWbsParameter& parameter,
int frameIndex, int frameIndex,
std::vector<double>* outputValues ) std::vector<double>* outputValues,
bool allowNormalization )
{ {
return calculateWbsParameterForAllSegments( parameter, m_parameterSources.at( parameter ), frameIndex, outputValues ); return calculateWbsParameterForAllSegments( parameter,
m_parameterSources.at( parameter ),
frameIndex,
outputValues,
allowNormalization );
} }
//-------------------------------------------------------------------------------------------------- //--------------------------------------------------------------------------------------------------
@ -377,7 +383,8 @@ std::vector<RigGeoMechWellLogExtractor::WbsParameterSource>
std::vector<RigGeoMechWellLogExtractor::WbsParameterSource> std::vector<RigGeoMechWellLogExtractor::WbsParameterSource>
RigGeoMechWellLogExtractor::calculateWbsParametersForAllSegments( const RigFemResultAddress& resAddr, RigGeoMechWellLogExtractor::calculateWbsParametersForAllSegments( const RigFemResultAddress& resAddr,
int frameIndex, int frameIndex,
std::vector<double>* values ) std::vector<double>* values,
bool allowNormalization )
{ {
CVF_ASSERT( values ); CVF_ASSERT( values );
@ -387,7 +394,7 @@ std::vector<RigGeoMechWellLogExtractor::WbsParameterSource>
CVF_ASSERT( false && "wbsParameters() called on something that isn't a wbs parameter" ); CVF_ASSERT( false && "wbsParameters() called on something that isn't a wbs parameter" );
} }
return calculateWbsParameterForAllSegments( param, m_userDefinedValues.at( param ), values ); return calculateWbsParameterForAllSegments( param, m_userDefinedValues.at( param ), values, allowNormalization );
} }
//-------------------------------------------------------------------------------------------------- //--------------------------------------------------------------------------------------------------
@ -474,11 +481,9 @@ std::vector<RigGeoMechWellLogExtractor::WbsParameterSource>
std::vector<double> ppShaleValues( m_intersections.size(), std::numeric_limits<double>::infinity() ); std::vector<double> ppShaleValues( m_intersections.size(), std::numeric_limits<double>::infinity() );
std::vector<WbsParameterSource> ppSandSources = std::vector<WbsParameterSource> ppSandSources =
calculateWbsParameterForAllSegments( RigWbsParameter::PP_Reservoir(), frameIndex, &ppSandValues ); calculateWbsParameterForAllSegments( RigWbsParameter::PP_Reservoir(), frameIndex, &ppSandValues, true );
std::vector<WbsParameterSource> ppShaleSources = std::vector<WbsParameterSource> ppShaleSources =
calculateWbsParameterForAllSegments( RigWbsParameter::PP_NonReservoir(), 0, &ppShaleValues ); calculateWbsParameterForAllSegments( RigWbsParameter::PP_NonReservoir(), 0, &ppShaleValues, true );
double waterDensityGCM3 = m_userDefinedValues[RigWbsParameter::waterDensity()];
#pragma omp parallel for #pragma omp parallel for
for ( int64_t intersectionIdx = 0; intersectionIdx < (int64_t)m_intersections.size(); ++intersectionIdx ) for ( int64_t intersectionIdx = 0; intersectionIdx < (int64_t)m_intersections.size(); ++intersectionIdx )
@ -497,8 +502,7 @@ std::vector<RigGeoMechWellLogExtractor::WbsParameterSource>
} }
else else
{ {
( *values )[intersectionIdx] = ( *values )[intersectionIdx] = 1.0;
hydroStaticPorePressureForIntersection( intersectionIdx, waterDensityGCM3 );
sources[intersectionIdx] = RigWbsParameter::HYDROSTATIC; sources[intersectionIdx] = RigWbsParameter::HYDROSTATIC;
} }
} }
@ -506,11 +510,11 @@ std::vector<RigGeoMechWellLogExtractor::WbsParameterSource>
} }
else if ( resAddr.fieldName == RiaDefines::wbsOBGResult().toStdString() ) else if ( resAddr.fieldName == RiaDefines::wbsOBGResult().toStdString() )
{ {
sources = calculateWbsParameterForAllSegments( RigWbsParameter::OBG(), frameIndex, values ); sources = calculateWbsParameterForAllSegments( RigWbsParameter::OBG(), frameIndex, values, true );
} }
else else
{ {
sources = calculateWbsParameterForAllSegments( RigWbsParameter::SH(), frameIndex, values ); sources = calculateWbsParameterForAllSegments( RigWbsParameter::SH(), frameIndex, values, true );
} }
return sources; return sources;
@ -555,19 +559,20 @@ void RigGeoMechWellLogExtractor::wellBoreWallCurveData( const RigFemResultAddres
std::vector<caf::Ten3d> interpolatedInterfaceStressBar = std::vector<caf::Ten3d> interpolatedInterfaceStressBar =
interpolateInterfaceValues( stressResAddr, frameIndex, vertexStresses ); interpolateInterfaceValues( stressResAddr, frameIndex, vertexStresses );
values->resize( m_intersections.size(), 0.0f ); values->resize( m_intersections.size(), std::numeric_limits<float>::infinity() );
std::vector<double> ppSandAllSegments( m_intersections.size(), std::numeric_limits<double>::infinity() ); std::vector<double> ppSandAllSegments( m_intersections.size(), std::numeric_limits<double>::infinity() );
std::vector<WbsParameterSource> ppSources = calculateWbsParameterForAllSegments( RigWbsParameter::PP_Reservoir(), std::vector<WbsParameterSource> ppSources = calculateWbsParameterForAllSegments( RigWbsParameter::PP_Reservoir(),
RigWbsParameter::GRID, RigWbsParameter::GRID,
frameIndex, frameIndex,
&ppSandAllSegments ); &ppSandAllSegments,
false );
std::vector<double> poissonAllSegments( m_intersections.size(), std::numeric_limits<double>::infinity() ); std::vector<double> poissonAllSegments( m_intersections.size(), std::numeric_limits<double>::infinity() );
calculateWbsParameterForAllSegments( RigWbsParameter::poissonRatio(), frameIndex, &poissonAllSegments ); calculateWbsParameterForAllSegments( RigWbsParameter::poissonRatio(), frameIndex, &poissonAllSegments, false );
std::vector<double> ucsAllSegments( m_intersections.size(), std::numeric_limits<double>::infinity() ); std::vector<double> ucsAllSegments( m_intersections.size(), std::numeric_limits<double>::infinity() );
calculateWbsParameterForAllSegments( RigWbsParameter::UCS(), frameIndex, &ucsAllSegments ); calculateWbsParameterForAllSegments( RigWbsParameter::UCS(), frameIndex, &ucsAllSegments, false );
#pragma omp parallel for #pragma omp parallel for
for ( int64_t intersectionIdx = 0; intersectionIdx < (int64_t)m_intersections.size(); ++intersectionIdx ) for ( int64_t intersectionIdx = 0; intersectionIdx < (int64_t)m_intersections.size(); ++intersectionIdx )
@ -638,8 +643,8 @@ void RigGeoMechWellLogExtractor::wellBoreFGShale( int frameIndex, std::vector<do
curveData( ppAddr, 0, &PP0 ); curveData( ppAddr, 0, &PP0 );
calculateWbsParameterForAllSegments( RigWbsParameter::K0_FG(), frameIndex, &K0_FG ); calculateWbsParameterForAllSegments( RigWbsParameter::K0_FG(), frameIndex, &K0_FG, true );
calculateWbsParameterForAllSegments( RigWbsParameter::OBG0(), 0, &OBG0 ); calculateWbsParameterForAllSegments( RigWbsParameter::OBG0(), 0, &OBG0, true );
#pragma omp parallel for #pragma omp parallel for
for ( int64_t intersectionIdx = 0; intersectionIdx < (int64_t)m_intersections.size(); ++intersectionIdx ) for ( int64_t intersectionIdx = 0; intersectionIdx < (int64_t)m_intersections.size(); ++intersectionIdx )
@ -658,7 +663,7 @@ void RigGeoMechWellLogExtractor::wellBoreFGShale( int frameIndex, std::vector<do
else else
{ {
std::vector<double> SH; std::vector<double> SH;
calculateWbsParameterForAllSegments( RigWbsParameter::SH(), frameIndex, &SH ); calculateWbsParameterForAllSegments( RigWbsParameter::SH(), frameIndex, &SH, true );
CVF_ASSERT( SH.size() == m_intersections.size() ); CVF_ASSERT( SH.size() == m_intersections.size() );
double multiplier = m_userDefinedValues.at( RigWbsParameter::FG_Shale() ); double multiplier = m_userDefinedValues.at( RigWbsParameter::FG_Shale() );
CVF_ASSERT( multiplier != std::numeric_limits<double>::infinity() ); CVF_ASSERT( multiplier != std::numeric_limits<double>::infinity() );
@ -689,9 +694,9 @@ void RigGeoMechWellLogExtractor::wellBoreSH_MatthewsKelly( int frameIndex, std::
curveData( ppAddr, frameIndex, &PP ); curveData( ppAddr, frameIndex, &PP );
curveData( ppAddr, 0, &PP0 ); curveData( ppAddr, 0, &PP0 );
calculateWbsParameterForAllSegments( RigWbsParameter::K0_SH(), frameIndex, &K0_SH ); calculateWbsParameterForAllSegments( RigWbsParameter::K0_SH(), frameIndex, &K0_SH, true );
calculateWbsParameterForAllSegments( RigWbsParameter::OBG0(), 0, &OBG0 ); calculateWbsParameterForAllSegments( RigWbsParameter::OBG0(), 0, &OBG0, true );
calculateWbsParameterForAllSegments( RigWbsParameter::DF(), frameIndex, &DF ); calculateWbsParameterForAllSegments( RigWbsParameter::DF(), frameIndex, &DF, true );
values->resize( m_intersections.size(), std::numeric_limits<double>::infinity() ); values->resize( m_intersections.size(), std::numeric_limits<double>::infinity() );
@ -779,7 +784,7 @@ std::vector<double> RigGeoMechWellLogExtractor::poissonSourceRegions( int frameI
{ {
std::vector<double> outputValues; std::vector<double> outputValues;
std::vector<WbsParameterSource> sources = std::vector<WbsParameterSource> sources =
calculateWbsParameterForAllSegments( RigWbsParameter::poissonRatio(), frameIndex, &outputValues ); calculateWbsParameterForAllSegments( RigWbsParameter::poissonRatio(), frameIndex, &outputValues, false );
std::vector<double> doubleSources( sources.size(), 0.0 ); std::vector<double> doubleSources( sources.size(), 0.0 );
#pragma omp parallel for #pragma omp parallel for
@ -797,7 +802,7 @@ std::vector<double> RigGeoMechWellLogExtractor::ucsSourceRegions( int frameIndex
{ {
std::vector<double> outputValues; std::vector<double> outputValues;
std::vector<WbsParameterSource> sources = std::vector<WbsParameterSource> sources =
calculateWbsParameterForAllSegments( RigWbsParameter::UCS(), frameIndex, &outputValues ); calculateWbsParameterForAllSegments( RigWbsParameter::UCS(), frameIndex, &outputValues, true );
std::vector<double> doubleSources( sources.size(), 0.0 ); std::vector<double> doubleSources( sources.size(), 0.0 );
#pragma omp parallel for #pragma omp parallel for

View File

@ -91,13 +91,16 @@ private:
std::vector<WbsParameterSource> calculateWbsParameterForAllSegments( const RigWbsParameter& parameter, std::vector<WbsParameterSource> calculateWbsParameterForAllSegments( const RigWbsParameter& parameter,
WbsParameterSource primarySource, WbsParameterSource primarySource,
int frameIndex, int frameIndex,
std::vector<double>* outputValues ); std::vector<double>* outputValues,
bool allowNormalization );
std::vector<WbsParameterSource> calculateWbsParameterForAllSegments( const RigWbsParameter& parameter, std::vector<WbsParameterSource> calculateWbsParameterForAllSegments( const RigWbsParameter& parameter,
int frameIndex, int frameIndex,
std::vector<double>* outputValues ); std::vector<double>* outputValues,
bool allowNormalization );
std::vector<WbsParameterSource> calculateWbsParametersForAllSegments( const RigFemResultAddress& resAddr, std::vector<WbsParameterSource> calculateWbsParametersForAllSegments( const RigFemResultAddress& resAddr,
int frameIndex, int frameIndex,
std::vector<double>* values ); std::vector<double>* values,
bool allowNormalization );
void wellPathAngles( const RigFemResultAddress& resAddr, std::vector<double>* values ); void wellPathAngles( const RigFemResultAddress& resAddr, std::vector<double>* values );
std::vector<WbsParameterSource> std::vector<WbsParameterSource>