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

Add optional max distance between curve points in WBS plots

Browse Source
This commit is contained in:
Magne Sjaastad 2023-03-31 16:27:14 +02:00 committed by GitHub
parent cde96a39f6
commit b2e8cc1663
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
11 changed files with 278 additions and 102 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
ApplicationLibCode/ProjectDataModel/RimWellBoreStabilityPlot.cpp
View File

@ -15,6 +15,7 @@
// for more details.
//
/////////////////////////////////////////////////////////////////////////////////
#include "RimWellBoreStabilityPlot.h"
#include "RiaDefines.h"

14
ApplicationLibCode/ProjectDataModel/WellLog/RimWellLogCurveCommonDataSource.cpp
View File

@ -106,6 +106,8 @@ RimWellLogCurveCommonDataSource::RimWellLogCurveCommonDataSource()
CAF_PDM_InitField( &m_wbsSmoothingThreshold, "WBSSmoothingThreshold", -1.0, "Smoothing Threshold" );
CAF_PDM_InitField( &m_maximumCurvePointInterval, "MaximumCurvePointInterval", std::make_pair( true, 10.0 ), "Maximum Curve Point Interval" );
CAF_PDM_InitFieldNoDefault( &m_rftTimeStep, "RftTimeStep", "RFT Time Step" );
CAF_PDM_InitFieldNoDefault( &m_rftWellName, "RftWellName", "RFT Well Name" );
CAF_PDM_InitFieldNoDefault( &m_rftSegmentBranchIndex, "SegmentBranchIndex", "RFT Branch" );
@ -590,15 +592,20 @@ void RimWellLogCurveCommonDataSource::applyDataSourceChanges( const std::vector<
if ( !wbsSmoothingToApply().isPartiallyTrue() )
{
wbsCurve->setSmoothCurve( wbsSmoothingToApply().isTrue() );
updatedSomething = true;
}
if ( wbsSmoothingThreshold() != 1.0 )
{
wbsCurve->setSmoothingThreshold( wbsSmoothingThreshold() );
updatedSomething = true;
}
wbsCurve->enableMaximumCurvePointInterval( m_maximumCurvePointInterval().first );
wbsCurve->setMaximumCurvePointInterval( m_maximumCurvePointInterval().second );
// Always do an update for wbs plots
updatedSomething = true;
}
if ( updatedSomething )
{
RimWellLogPlot* parentPlot = nullptr;
@ -1058,6 +1065,7 @@ void RimWellLogCurveCommonDataSource::defineUiOrdering( QString uiConfigName, ca
{
group->add( &m_wbsSmoothing );
group->add( &m_wbsSmoothingThreshold );
group->add( &m_maximumCurvePointInterval );
}
if ( !m_uniqueRftTimeSteps.empty() ) group->add( &m_rftTimeStep );
@ -1110,7 +1118,7 @@ void RimWellLogCurveCommonDataSource::defineEditorAttribute( const caf::PdmField
}
}
auto* uiDisplayStringAttr = dynamic_cast<caf::PdmUiLineEditorAttributeUiDisplayString*>( attribute );
if ( uiDisplayStringAttr && wbsSmoothingThreshold() == -1.0 )
if ( uiDisplayStringAttr && ( wbsSmoothingThreshold() == -1.0 ) && ( field == &m_wbsSmoothingThreshold ) )
{
QString displayString = "Mixed";

6
ApplicationLibCode/ProjectDataModel/WellLog/RimWellLogCurveCommonDataSource.h
View File

@ -125,8 +125,10 @@ private:
caf::PdmField<int> m_branchIndex;
caf::PdmField<caf::Tristate> m_branchDetection;
caf::PdmField<int> m_timeStep;
caf::PdmField<caf::Tristate> m_wbsSmoothing;
caf::PdmField<double> m_wbsSmoothingThreshold;
caf::PdmField<caf::Tristate> m_wbsSmoothing;
caf::PdmField<double> m_wbsSmoothingThreshold;
caf::PdmField<std::pair<bool, double>> m_maximumCurvePointInterval;
caf::PdmField<QDateTime> m_rftTimeStep;
caf::PdmField<QString> m_rftWellName;

56
ApplicationLibCode/ProjectDataModel/WellLog/RimWellLogExtractionCurve.cpp
View File

@ -421,13 +421,15 @@ void RimWellLogExtractionCurve::onLoadDataAndUpdate( bool updateParentPlot )
//--------------------------------------------------------------------------------------------------
void RimWellLogExtractionCurve::performDataExtraction( bool* isUsingPseudoLength )
{
extractData( isUsingPseudoLength );
extractData( isUsingPseudoLength, {}, {} );
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RimWellLogExtractionCurve::extractData( bool* isUsingPseudoLength, bool performDataSmoothing /*= false*/, double smoothingThreshold /*= -1.0 */ )
void RimWellLogExtractionCurve::extractData( bool* isUsingPseudoLength,
const std::optional<double>& smoothingThreshold,
const std::optional<double>& maxDistanceBetweenCurvePoints )
{
CAF_ASSERT( isUsingPseudoLength );
@ -447,13 +449,15 @@ void RimWellLogExtractionCurve::extractData( bool* isUsingPseudoLength, bool per
}
else if ( geomCase && geomCase->geoMechData() )
{
curveData = extractGeomData( geomCase, isUsingPseudoLength, performDataSmoothing, smoothingThreshold );
curveData = extractGeomData( geomCase, isUsingPseudoLength, smoothingThreshold, maxDistanceBetweenCurvePoints );
}
if ( !curveData.values.empty() && !curveData.measuredDepthValues.empty() )
{
bool useLogarithmicScale = false;
bool performDataSmoothing = smoothingThreshold.has_value();
RimWellLogTrack* track = nullptr;
firstAncestorOfType( track );
if ( track )
@ -592,17 +596,37 @@ RimWellLogExtractionCurve::WellLogExtractionCurveData RimWellLogExtractionCurve:
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
RimWellLogExtractionCurve::WellLogExtractionCurveData RimWellLogExtractionCurve::extractGeomData( RimGeoMechCase* geomCase,
bool* isUsingPseudoLength,
bool performDataSmoothing,
double smoothingThreshold )
RimWellLogExtractionCurve::WellLogExtractionCurveData
RimWellLogExtractionCurve::extractGeomData( RimGeoMechCase* geoMechCase,
bool* isUsingPseudoLength,
const std::optional<double>& smoothingThreshold,
const std::optional<double>& maxDistanceBetweenCurvePoints )
{
WellLogExtractionCurveData curveData;
RimWellLogPlotCollection* wellLogCollection = RimMainPlotCollection::current()->wellLogPlotCollection();
cvf::ref<RigGeoMechWellLogExtractor> wellExtractor = wellLogCollection->findOrCreateExtractor( m_wellPath, geomCase );
cvf::ref<RigGeoMechWellLogExtractor> refWellExtractor = wellLogCollection->findOrCreateExtractor( m_refWellPath, geomCase );
WellLogExtractionCurveData curveData;
RimWellLogPlotCollection* wellLogCollection = RimMainPlotCollection::current()->wellLogPlotCollection();
auto [timeStepIdx, frameIdx] = geomCase->geoMechData()->femPartResults()->stepListIndexToTimeStepAndDataFrameIndex( m_timeStep );
cvf::ref<RigGeoMechWellLogExtractor> wellExtractor;
if ( maxDistanceBetweenCurvePoints.has_value() && maxDistanceBetweenCurvePoints.value() > 0.0 )
{
RigGeoMechCaseData* caseData = geoMechCase->geoMechData();
auto wellPathGeometry = m_wellPath->wellPathGeometry();
if ( caseData && wellPathGeometry )
{
std::string errorIdName = ( m_wellPath->name() + " " + geoMechCase->caseUserDescription() ).toStdString();
wellExtractor = new RigGeoMechWellLogExtractor( caseData, wellPathGeometry, errorIdName );
// make sure the resampling of the well path is done before the extraction of the curve data
wellExtractor->resampleIntersections( maxDistanceBetweenCurvePoints.value() );
}
}
else
{
wellExtractor = wellLogCollection->findOrCreateExtractor( m_wellPath, geoMechCase );
}
cvf::ref<RigGeoMechWellLogExtractor> refWellExtractor = wellLogCollection->findOrCreateExtractor( m_refWellPath, geoMechCase );
auto [timeStepIdx, frameIdx] = geoMechCase->geoMechData()->femPartResults()->stepListIndexToTimeStepAndDataFrameIndex( m_timeStep );
if ( wellExtractor.notNull() )
{
@ -659,26 +683,26 @@ RimWellLogExtractionCurve::WellLogExtractionCurveData RimWellLogExtractionCurve:
refWellPropertyValues,
refWellIndexKValues );
}
if ( performDataSmoothing )
if ( smoothingThreshold.has_value() )
{
refWellExtractor->performCurveDataSmoothing( timeStepIdx,
frameIdx,
&curveData.measuredDepthValues,
&curveData.tvDepthValues,
&curveData.values,
smoothingThreshold );
smoothingThreshold.value() );
}
return curveData;
}
if ( wellExtractor.notNull() && performDataSmoothing )
if ( wellExtractor.notNull() && smoothingThreshold.has_value() )
{
wellExtractor->performCurveDataSmoothing( timeStepIdx,
frameIdx,
&curveData.measuredDepthValues,
&curveData.tvDepthValues,
&curveData.values,
smoothingThreshold );
smoothingThreshold.value() );
}
return curveData;
}

29
ApplicationLibCode/ProjectDataModel/WellLog/RimWellLogExtractionCurve.h
View File

@ -105,7 +105,9 @@ protected:
void connectCaseSignals( RimCase* rimCase );
virtual void performDataExtraction( bool* isUsingPseudoLength );
void extractData( bool* isUsingPseudoLength, bool performDataSmoothing = false, double smoothingThreshold = -1.0 );
void extractData( bool* isUsingPseudoLength,
const std::optional<double>& smoothingThreshold,
const std::optional<double>& maxDistanceBetweenCurvePoints );
void fieldChangedByUi( const caf::PdmFieldHandle* changedField, const QVariant& oldValue, const QVariant& newValue ) override;
void defineUiOrdering( QString uiConfigName, caf::PdmUiOrdering& uiOrdering ) override;
@ -154,16 +156,17 @@ private:
private:
WellLogExtractionCurveData extractEclipseData( RimEclipseCase* eclipseCase, bool* isUsingPseudoLength );
WellLogExtractionCurveData extractGeomData( RimGeoMechCase* geomCase,
bool* isUsingPseudoLength,
bool performDataSmoothing = false,
double smoothingThreshold = -1.0 );
void mapPropertyValuesFromReferenceWell( std::vector<double>& rMeasuredDepthValues,
std::vector<double>& rTvDepthValues,
std::vector<double>& rPropertyValues,
const std::vector<double>& indexKValues,
const std::vector<double>& refWellMeasuredDepthValues,
const std::vector<double>& refWellTvDepthValues,
const std::vector<double>& refWellPropertyValues,
const std::vector<double>& refWellIndexKValues );
WellLogExtractionCurveData extractGeomData( RimGeoMechCase* geoMechCase,
bool* isUsingPseudoLength,
const std::optional<double>& smoothingThreshold,
const std::optional<double>& maxDistanceBetweenCurvePoints );
void mapPropertyValuesFromReferenceWell( std::vector<double>& rMeasuredDepthValues,
std::vector<double>& rTvDepthValues,
std::vector<double>& rPropertyValues,
const std::vector<double>& indexKValues,
const std::vector<double>& refWellMeasuredDepthValues,
const std::vector<double>& refWellTvDepthValues,
const std::vector<double>& refWellPropertyValues,
const std::vector<double>& refWellIndexKValues );
};

26
ApplicationLibCode/ProjectDataModel/WellLog/RimWellLogWbsCurve.cpp
View File

@ -16,6 +16,8 @@ RimWellLogWbsCurve::RimWellLogWbsCurve()
CAF_PDM_InitField( &m_smoothCurve, "SmoothCurve", false, "Smooth Curve" );
CAF_PDM_InitField( &m_smoothingThreshold, "SmoothingThreshold", 0.002, "Smoothing Threshold" );
CAF_PDM_InitField( &m_maximumCurvePointInterval, "MaximumCurvePointInterval", std::make_pair( false, 10.0 ), "Maximum Curve Point Interval" );
}
//--------------------------------------------------------------------------------------------------
@ -50,12 +52,31 @@ void RimWellLogWbsCurve::setSmoothingThreshold( double threshold )
m_smoothingThreshold = threshold;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RimWellLogWbsCurve::enableMaximumCurvePointInterval( bool enable )
{
m_maximumCurvePointInterval.v().first = enable;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RimWellLogWbsCurve::setMaximumCurvePointInterval( double interval )
{
m_maximumCurvePointInterval.v().second = interval;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RimWellLogWbsCurve::performDataExtraction( bool* isUsingPseudoLength )
{
extractData( isUsingPseudoLength, m_smoothCurve(), m_smoothingThreshold() );
auto smoothingThreshold = createOptional( m_smoothCurve(), m_smoothingThreshold() );
auto maxDistanceBetweenPoints = createOptional( m_maximumCurvePointInterval() );
extractData( isUsingPseudoLength, smoothingThreshold, maxDistanceBetweenPoints );
}
//--------------------------------------------------------------------------------------------------
@ -68,6 +89,7 @@ void RimWellLogWbsCurve::defineUiOrdering( QString uiConfigName, caf::PdmUiOrder
caf::PdmUiGroup* dataGroup = uiOrdering.findGroup( dataSourceGroupKeyword() );
dataGroup->add( &m_smoothCurve );
dataGroup->add( &m_smoothingThreshold );
dataGroup->add( &m_maximumCurvePointInterval );
}
//--------------------------------------------------------------------------------------------------
@ -77,7 +99,7 @@ void RimWellLogWbsCurve::fieldChangedByUi( const caf::PdmFieldHandle* changedFie
{
RimWellLogExtractionCurve::fieldChangedByUi( changedField, oldValue, newValue );
if ( changedField == &m_smoothCurve || changedField == &m_smoothingThreshold )
if ( changedField == &m_smoothCurve || changedField == &m_smoothingThreshold || changedField == &m_maximumCurvePointInterval )
{
this->loadDataAndUpdate( true );
}

22
ApplicationLibCode/ProjectDataModel/WellLog/RimWellLogWbsCurve.h
View File

@ -19,6 +19,23 @@
#include "RimWellLogExtractionCurve.h"
//==================================================================================================
/// Helpers for creating optional values, move to stdOptionalTools when used in more places
//==================================================================================================
template <typename T>
std::optional<T> createOptional( bool enable, const T& value )
{
if ( !enable ) return {};
return std::optional<T>( value );
}
template <typename T>
std::optional<T> createOptional( const std::pair<bool, T>& value )
{
return createOptional( value.first, value.second );
}
//==================================================================================================
///
///
@ -36,12 +53,17 @@ public:
void setSmoothCurve( bool smooth );
void setSmoothingThreshold( double threshold );
void enableMaximumCurvePointInterval( bool enable );
void setMaximumCurvePointInterval( double interval );
protected:
void performDataExtraction( bool* isUsingPseudoLength ) override;
void defineUiOrdering( QString uiConfigName, caf::PdmUiOrdering& uiOrdering ) override;
void fieldChangedByUi( const caf::PdmFieldHandle* changedField, const QVariant& oldValue, const QVariant& newValue ) override;
protected:
caf::PdmField<std::pair<bool, double>> m_maximumCurvePointInterval;
caf::PdmField<bool> m_smoothCurve;
caf::PdmField<double> m_smoothingThreshold;
};

8
ApplicationLibCode/ReservoirDataModel/RigEclipseWellLogExtractor.cpp
View File

@ -177,12 +177,12 @@ void RigEclipseWellLogExtractor::calculateIntersection()
void RigEclipseWellLogExtractor::curveData( const RigResultAccessor* resultAccessor, std::vector<double>* values )
{
CVF_TIGHT_ASSERT( values );
values->resize( m_intersections.size() );
values->resize( intersections().size() );
for ( size_t cpIdx = 0; cpIdx < m_intersections.size(); ++cpIdx )
for ( size_t cpIdx = 0; cpIdx < intersections().size(); ++cpIdx )
{
size_t cellIdx = m_intersectedCellsGlobIdx[cpIdx];
cvf::StructGridInterface::FaceType cellFace = m_intersectedCellFaces[cpIdx];
size_t cellIdx = intersectedCellsGlobIdx()[cpIdx];
cvf::StructGridInterface::FaceType cellFace = intersectedCellFaces()[cpIdx];
( *values )[cpIdx] = resultAccessor->cellFaceScalarGlobIdx( cellIdx, cellFace );
}
}

105
ApplicationLibCode/ReservoirDataModel/RigGeoMechWellLogExtractor.cpp
View File

@ -99,7 +99,7 @@ void RigGeoMechWellLogExtractor::performCurveDataSmoothing( int
std::vector<double> interfaceShValuesDbl( interfaceShValues.size(), std::numeric_limits<double>::infinity() );
std::vector<double> interfacePorePressuresDbl( interfacePorePressures.size(), std::numeric_limits<double>::infinity() );
#pragma omp parallel for
for ( int64_t i = 0; i < int64_t( m_intersections.size() ); ++i )
for ( int64_t i = 0; i < static_cast<int64_t>( intersections().size() ); ++i )
{
double hydroStaticPorePressureBar = hydroStaticPorePressureForSegment( i );
interfaceShValuesDbl[i] = interfaceShValues[i] / hydroStaticPorePressureBar;
@ -219,7 +219,7 @@ QString RigGeoMechWellLogExtractor::curveData( const RigFemResultAddress& resAdd
values->resize( interfaceValues.size(), std::numeric_limits<double>::infinity() );
#pragma omp parallel for
for ( int64_t intersectionIdx = 0; intersectionIdx < (int64_t)m_intersections.size(); ++intersectionIdx )
for ( int64_t intersectionIdx = 0; intersectionIdx < static_cast<int64_t>( intersections().size() ); ++intersectionIdx )
{
( *values )[intersectionIdx] = static_cast<double>( interfaceValues[intersectionIdx] );
}
@ -241,7 +241,7 @@ std::vector<RigGeoMechWellLogExtractor::WbsParameterSource>
{
RigFemPartResultsCollection* resultCollection = m_caseData->femPartResults();
std::vector<WbsParameterSource> finalSourcesPerSegment( m_intersections.size(), RigWbsParameter::UNDEFINED );
std::vector<WbsParameterSource> finalSourcesPerSegment( intersections().size(), RigWbsParameter::UNDEFINED );
if ( primarySource == RigWbsParameter::UNDEFINED )
{
@ -264,10 +264,10 @@ std::vector<RigGeoMechWellLogExtractor::WbsParameterSource>
const std::vector<float>& unscaledResultValues = resultCollection->resultValues( nativeAddr, 0, timeStepIndex, frameIndex );
std::vector<float> interpolatedInterfaceValues =
interpolateInterfaceValues( nativeAddr, timeStepIndex, frameIndex, unscaledResultValues );
gridValues.resize( m_intersections.size(), std::numeric_limits<double>::infinity() );
gridValues.resize( intersections().size(), std::numeric_limits<double>::infinity() );
#pragma omp parallel for
for ( int64_t intersectionIdx = 0; intersectionIdx < (int64_t)m_intersections.size(); ++intersectionIdx )
for ( int64_t intersectionIdx = 0; intersectionIdx < static_cast<int64_t>( intersections().size() ); ++intersectionIdx )
{
float averageUnscaledValue = std::numeric_limits<float>::infinity();
averageIntersectionValuesToSegmentValue( intersectionIdx,
@ -287,7 +287,7 @@ std::vector<RigGeoMechWellLogExtractor::WbsParameterSource>
const std::vector<float>* elementPropertyValuesInput = nullptr;
std::vector<float> tvdRKBs;
for ( double tvdValue : m_intersectionTVDs )
for ( double tvdValue : cellIntersectionTVDs() )
{
tvdRKBs.push_back( tvdValue + m_wellPathGeometry->rkbDiff() );
}
@ -303,18 +303,19 @@ std::vector<RigGeoMechWellLogExtractor::WbsParameterSource>
}
}
std::vector<double> unscaledValues( m_intersections.size(), std::numeric_limits<double>::infinity() );
std::vector<double> unscaledValues( intersections().size(), std::numeric_limits<double>::infinity() );
double waterDensityGCM3 = m_userDefinedValues[RigWbsParameter::waterDensity()];
for ( int64_t intersectionIdx = 0; intersectionIdx < (int64_t)m_intersections.size(); ++intersectionIdx )
for ( int64_t intersectionIdx = 0; intersectionIdx < static_cast<int64_t>( intersections().size() ); ++intersectionIdx )
{
// Loop from primary source and out for each value
for ( auto it = primary_it; it != allSources.end(); ++it )
{
if ( *it == RigWbsParameter::GRID ) // Priority 0: Grid
{
if ( intersectionIdx < (int64_t)gridValues.size() && gridValues[intersectionIdx] != std::numeric_limits<double>::infinity() )
if ( intersectionIdx < static_cast<int64_t>( gridValues.size() ) &&
gridValues[intersectionIdx] != std::numeric_limits<double>::infinity() )
{
unscaledValues[intersectionIdx] = gridValues[intersectionIdx];
finalSourcesPerSegment[intersectionIdx] = RigWbsParameter::GRID;
@ -346,7 +347,7 @@ std::vector<RigGeoMechWellLogExtractor::WbsParameterSource>
{
if ( !elementPropertyValues.empty() )
{
size_t elmIdx = m_intersectedCellsGlobIdx[intersectionIdx];
size_t elmIdx = intersectedCellsGlobIdx()[intersectionIdx];
if ( elmIdx < elementPropertyValues.size() )
{
unscaledValues[intersectionIdx] = elementPropertyValues[elmIdx];
@ -376,7 +377,7 @@ std::vector<RigGeoMechWellLogExtractor::WbsParameterSource>
outputValues->resize( unscaledValues.size(), std::numeric_limits<double>::infinity() );
#pragma omp parallel for
for ( int64_t intersectionIdx = 0; intersectionIdx < (int64_t)m_intersections.size(); ++intersectionIdx )
for ( int64_t intersectionIdx = 0; intersectionIdx < static_cast<int64_t>( intersections().size() ); ++intersectionIdx )
{
RigWbsParameter::Source source = finalSourcesPerSegment[intersectionIdx];
@ -444,12 +445,12 @@ void RigGeoMechWellLogExtractor::wellPathAngles( const RigFemResultAddress& resA
{
CVF_ASSERT( values );
CVF_ASSERT( resAddr.fieldName == "Azimuth" || resAddr.fieldName == "Inclination" );
values->resize( m_intersections.size(), 0.0f );
values->resize( intersections().size(), 0.0f );
const double epsilon = 1.0e-6 * 360;
const cvf::Vec3d trueNorth( 0.0, 1.0, 0.0 );
const cvf::Vec3d up( 0.0, 0.0, 1.0 );
double previousAzimuth = 0.0;
for ( int64_t intersectionIdx = 0; intersectionIdx < (int64_t)m_intersections.size(); ++intersectionIdx )
for ( int64_t intersectionIdx = 0; intersectionIdx < static_cast<int64_t>( intersections().size() ); ++intersectionIdx )
{
cvf::Vec3d wellPathTangent = calculateWellPathTangent( intersectionIdx, TangentFollowWellPathSegments );
@ -512,14 +513,14 @@ std::vector<RigGeoMechWellLogExtractor::WbsParameterSource>
{
CVF_ASSERT( values );
values->resize( m_intersections.size(), std::numeric_limits<double>::infinity() );
std::vector<WbsParameterSource> sources( m_intersections.size(), RigWbsParameter::UNDEFINED );
values->resize( intersections().size(), std::numeric_limits<double>::infinity() );
std::vector<WbsParameterSource> sources( intersections().size(), RigWbsParameter::UNDEFINED );
if ( resAddr.fieldName == RiaResultNames::wbsPPResult().toStdString() )
{
// Las or element property table values
std::vector<double> ppSandValues( m_intersections.size(), std::numeric_limits<double>::infinity() );
std::vector<double> ppShaleValues( m_intersections.size(), std::numeric_limits<double>::infinity() );
std::vector<double> ppSandValues( intersections().size(), std::numeric_limits<double>::infinity() );
std::vector<double> ppShaleValues( intersections().size(), std::numeric_limits<double>::infinity() );
std::vector<WbsParameterSource> ppSandSources;
if ( forceGridSourceForPPReservoir )
@ -537,7 +538,7 @@ std::vector<RigGeoMechWellLogExtractor::WbsParameterSource>
calculateWbsParameterForAllSegments( RigWbsParameter::PP_NonReservoir(), 0, 0, &ppShaleValues, true );
#pragma omp parallel for
for ( int64_t intersectionIdx = 0; intersectionIdx < (int64_t)m_intersections.size(); ++intersectionIdx )
for ( int64_t intersectionIdx = 0; intersectionIdx < static_cast<int64_t>( intersections().size() ); ++intersectionIdx )
{
if ( ( *values )[intersectionIdx] == std::numeric_limits<double>::infinity() )
{
@ -603,20 +604,20 @@ void RigGeoMechWellLogExtractor::wellBoreWallCurveData( const RigFemResultAddres
std::vector<caf::Ten3d> interpolatedInterfaceStressBar =
interpolateInterfaceValues( stressResAddr, timeStepIndex, frameIndex, vertexStresses );
values->resize( m_intersections.size(), std::numeric_limits<float>::infinity() );
values->resize( intersections().size(), std::numeric_limits<float>::infinity() );
std::vector<double> ppSandAllSegments( m_intersections.size(), std::numeric_limits<double>::infinity() );
std::vector<double> ppSandAllSegments( intersections().size(), std::numeric_limits<double>::infinity() );
std::vector<WbsParameterSource> ppSources =
calculateWbsParameterForAllSegments( RigWbsParameter::PP_Reservoir(), RigWbsParameter::GRID, frameIndex, &ppSandAllSegments, false );
std::vector<double> poissonAllSegments( m_intersections.size(), std::numeric_limits<double>::infinity() );
std::vector<double> poissonAllSegments( intersections().size(), std::numeric_limits<double>::infinity() );
calculateWbsParameterForAllSegments( RigWbsParameter::poissonRatio(), timeStepIndex, frameIndex, &poissonAllSegments, false );
std::vector<double> ucsAllSegments( m_intersections.size(), std::numeric_limits<double>::infinity() );
std::vector<double> ucsAllSegments( intersections().size(), std::numeric_limits<double>::infinity() );
calculateWbsParameterForAllSegments( RigWbsParameter::UCS(), timeStepIndex, frameIndex, &ucsAllSegments, false );
#pragma omp parallel for
for ( int64_t intersectionIdx = 0; intersectionIdx < (int64_t)m_intersections.size(); ++intersectionIdx )
for ( int64_t intersectionIdx = 0; intersectionIdx < static_cast<int64_t>( intersections().size() ); ++intersectionIdx )
{
// FG is for sands, SFG for shale. Sands has valid PP, shale does not.
bool isFGregion = ppSources[intersectionIdx] == RigWbsParameter::GRID;
@ -673,7 +674,7 @@ void RigGeoMechWellLogExtractor::wellBoreWallCurveData( const RigFemResultAddres
//--------------------------------------------------------------------------------------------------
void RigGeoMechWellLogExtractor::wellBoreFGShale( int timeStepIndex, int frameIndex, std::vector<double>* values )
{
if ( values->empty() ) values->resize( m_intersections.size(), std::numeric_limits<double>::infinity() );
if ( values->empty() ) values->resize( intersections().size(), std::numeric_limits<double>::infinity() );
WbsParameterSource source = m_parameterSources.at( RigWbsParameter::FG_Shale() );
if ( source == RigWbsParameter::DERIVED_FROM_K0FG )
@ -688,7 +689,7 @@ void RigGeoMechWellLogExtractor::wellBoreFGShale( int timeStepIndex, int frameIn
calculateWbsParameterForAllSegments( RigWbsParameter::OBG0(), 0, 0, &OBG0, true );
#pragma omp parallel for
for ( int64_t intersectionIdx = 0; intersectionIdx < (int64_t)m_intersections.size(); ++intersectionIdx )
for ( int64_t intersectionIdx = 0; intersectionIdx < static_cast<int64_t>( intersections().size() ); ++intersectionIdx )
{
if ( !isValid( ( *values )[intersectionIdx] ) )
{
@ -704,11 +705,11 @@ void RigGeoMechWellLogExtractor::wellBoreFGShale( int timeStepIndex, int frameIn
{
std::vector<double> SH;
calculateWbsParameterForAllSegments( RigWbsParameter::SH(), timeStepIndex, frameIndex, &SH, true );
CVF_ASSERT( SH.size() == m_intersections.size() );
CVF_ASSERT( SH.size() == intersections().size() );
double multiplier = m_userDefinedValues.at( RigWbsParameter::FG_Shale() );
CVF_ASSERT( multiplier != std::numeric_limits<double>::infinity() );
#pragma omp parallel for
for ( int64_t intersectionIdx = 0; intersectionIdx < (int64_t)m_intersections.size(); ++intersectionIdx )
for ( int64_t intersectionIdx = 0; intersectionIdx < static_cast<int64_t>( intersections().size() ); ++intersectionIdx )
{
if ( !isValid( ( *values )[intersectionIdx] ) )
{
@ -738,10 +739,10 @@ void RigGeoMechWellLogExtractor::wellBoreSH_MatthewsKelly( int timeStepIndex, in
calculateWbsParameterForAllSegments( RigWbsParameter::OBG0(), 0, 0, &OBG0, true );
calculateWbsParameterForAllSegments( RigWbsParameter::DF(), timeStepIndex, frameIndex, &DF, true );
values->resize( m_intersections.size(), std::numeric_limits<double>::infinity() );
values->resize( intersections().size(), std::numeric_limits<double>::infinity() );
#pragma omp parallel for
for ( int64_t intersectionIdx = 0; intersectionIdx < (int64_t)m_intersections.size(); ++intersectionIdx )
for ( int64_t intersectionIdx = 0; intersectionIdx < static_cast<int64_t>( intersections().size() ); ++intersectionIdx )
{
if ( isValid( PP[intersectionIdx] ) && isValid( PP0[intersectionIdx] ) && isValid( OBG0[intersectionIdx] ) &&
isValid( K0_SH[intersectionIdx] ) && isValid( DF[intersectionIdx] ) )
@ -816,7 +817,7 @@ std::vector<double> RigGeoMechWellLogExtractor::porePressureSourceRegions( int t
std::vector<double> doubleSources( sources.size(), 0.0 );
#pragma omp parallel for
for ( int64_t intersectionIdx = 0; intersectionIdx < (int64_t)m_intersections.size(); ++intersectionIdx )
for ( int64_t intersectionIdx = 0; intersectionIdx < static_cast<int64_t>( intersections().size() ); ++intersectionIdx )
{
doubleSources[intersectionIdx] = static_cast<double>( sources[intersectionIdx] );
}
@ -833,7 +834,7 @@ std::vector<double> RigGeoMechWellLogExtractor::poissonSourceRegions( int timeSt
std::vector<double> doubleSources( sources.size(), 0.0 );
#pragma omp parallel for
for ( int64_t intersectionIdx = 0; intersectionIdx < (int64_t)m_intersections.size(); ++intersectionIdx )
for ( int64_t intersectionIdx = 0; intersectionIdx < static_cast<int64_t>( intersections().size() ); ++intersectionIdx )
{
doubleSources[intersectionIdx] = static_cast<double>( sources[intersectionIdx] );
}
@ -851,7 +852,7 @@ std::vector<double> RigGeoMechWellLogExtractor::ucsSourceRegions( int timeStepIn
std::vector<double> doubleSources( sources.size(), 0.0 );
#pragma omp parallel for
for ( int64_t intersectionIdx = 0; intersectionIdx < (int64_t)m_intersections.size(); ++intersectionIdx )
for ( int64_t intersectionIdx = 0; intersectionIdx < static_cast<int64_t>( intersections().size() ); ++intersectionIdx )
{
doubleSources[intersectionIdx] = static_cast<double>( sources[intersectionIdx] );
}
@ -869,7 +870,7 @@ T RigGeoMechWellLogExtractor::interpolateGridResultValue( RigFemResultPosEnum
const RigFemPart* femPart = m_caseData->femParts()->part( 0 );
const std::vector<cvf::Vec3f>& nodeCoords = femPart->nodes().coordinates;
size_t elmIdx = m_intersectedCellsGlobIdx[intersectionIdx];
size_t elmIdx = intersectedCellsGlobIdx()[intersectionIdx];
RigElementType elmType = femPart->elementType( elmIdx );
if ( !( elmType == HEX8 || elmType == HEX8P ) ) return T();
@ -884,7 +885,7 @@ T RigGeoMechWellLogExtractor::interpolateGridResultValue( RigFemResultPosEnum
return gridResultValues[elmIdx];
}
cvf::StructGridInterface::FaceType cellFace = m_intersectedCellFaces[intersectionIdx];
cvf::StructGridInterface::FaceType cellFace = intersectedCellFaces()[intersectionIdx];
if ( cellFace == cvf::StructGridInterface::NO_FACE )
{
@ -943,7 +944,7 @@ T RigGeoMechWellLogExtractor::interpolateGridResultValue( RigFemResultPosEnum
nodeResultValues[2],
v3,
nodeResultValues[3],
m_intersections[intersectionIdx] );
intersections()[intersectionIdx] );
return interpolatedValue;
}
@ -1065,7 +1066,7 @@ cvf::Vec3d RigGeoMechWellLogExtractor::calculateWellPathTangent( int64_t interse
if ( calculationType == TangentFollowWellPathSegments )
{
cvf::Vec3d segmentStart, segmentEnd;
m_wellPathGeometry->twoClosestPoints( m_intersections[intersectionIdx], &segmentStart, &segmentEnd );
m_wellPathGeometry->twoClosestPoints( intersections()[intersectionIdx], &segmentStart, &segmentEnd );
return ( segmentEnd - segmentStart ).getNormalized();
}
else
@ -1073,11 +1074,11 @@ cvf::Vec3d RigGeoMechWellLogExtractor::calculateWellPathTangent( int64_t interse
cvf::Vec3d wellPathTangent;
if ( intersectionIdx % 2 == 0 )
{
wellPathTangent = m_intersections[intersectionIdx + 1] - m_intersections[intersectionIdx];
wellPathTangent = intersections()[intersectionIdx + 1] - intersections()[intersectionIdx];
}
else
{
wellPathTangent = m_intersections[intersectionIdx] - m_intersections[intersectionIdx - 1];
wellPathTangent = intersections()[intersectionIdx] - intersections()[intersectionIdx - 1];
}
CVF_ASSERT( wellPathTangent.length() > 1.0e-7 );
return wellPathTangent.getNormalized();
@ -1107,7 +1108,7 @@ cvf::Vec3f RigGeoMechWellLogExtractor::cellCentroid( size_t intersectionIdx ) co
const RigFemPart* femPart = m_caseData->femParts()->part( 0 );
const std::vector<cvf::Vec3f>& nodeCoords = femPart->nodes().coordinates;
size_t elmIdx = m_intersectedCellsGlobIdx[intersectionIdx];
size_t elmIdx = intersectedCellsGlobIdx()[intersectionIdx];
RigElementType elmType = femPart->elementType( elmIdx );
int elementNodeCount = RigFemTypes::elementNodeCount( elmType );
@ -1129,7 +1130,7 @@ double RigGeoMechWellLogExtractor::getWellLogIntersectionValue( size_t
{
const double eps = 1.0e-4;
double intersection_md = m_intersectionMeasuredDepths[intersectionIdx];
double intersection_md = cellIntersectionMDs()[intersectionIdx];
for ( size_t i = 0; i < wellLogValues.size() - 1; ++i )
{
double las_md_i = wellLogValues[i].first;
@ -1203,16 +1204,16 @@ bool RigGeoMechWellLogExtractor::averageIntersectionValuesToSegmentValue( size_t
value1 = values[intersectionIdx];
value2 = values[intersectionIdx + 1];
dist1 = ( centroid - m_intersections[intersectionIdx] ).length();
dist2 = ( centroid - m_intersections[intersectionIdx + 1] ).length();
dist1 = ( centroid - intersections()[intersectionIdx] ).length();
dist2 = ( centroid - intersections()[intersectionIdx + 1] ).length();
}
else
{
value1 = values[intersectionIdx - 1];
value2 = values[intersectionIdx];
dist1 = ( centroid - m_intersections[intersectionIdx - 1] ).length();
dist2 = ( centroid - m_intersections[intersectionIdx] ).length();
dist1 = ( centroid - intersections()[intersectionIdx - 1] ).length();
dist2 = ( centroid - intersections()[intersectionIdx] ).length();
}
if ( invalidValue == value1 || invalidValue == value2 )
@ -1240,14 +1241,14 @@ std::vector<T> RigGeoMechWellLogExtractor::interpolateInterfaceValues( RigFemRes
const std::vector<T>& unscaledResultValues )
{
std::vector<T> interpolatedInterfaceValues;
initializeResultValues( interpolatedInterfaceValues, m_intersections.size() );
initializeResultValues( interpolatedInterfaceValues, intersections().size() );
const RigFemPart* femPart = m_caseData->femParts()->part( 0 );
#pragma omp parallel for
for ( int64_t intersectionIdx = 0; intersectionIdx < (int64_t)m_intersections.size(); ++intersectionIdx )
for ( int64_t intersectionIdx = 0; intersectionIdx < static_cast<int64_t>( intersections().size() ); ++intersectionIdx )
{
size_t elmIdx = m_intersectedCellsGlobIdx[intersectionIdx];
size_t elmIdx = intersectedCellsGlobIdx()[intersectionIdx];
RigElementType elmType = femPart->elementType( elmIdx );
if ( !( elmType == HEX8 || elmType == HEX8P ) ) continue;
@ -1287,7 +1288,7 @@ void RigGeoMechWellLogExtractor::smoothSegments( std::vector<double>*
double maxOriginalMd = ( *mds )[0];
double maxOriginalTvd = ( !tvds->empty() ) ? ( *tvds )[0] : 0.0;
for ( int64_t i = 1; i < int64_t( mds->size() - 1 ); ++i )
for ( int64_t i = 1; i < static_cast<int64_t>( mds->size() - 1 ); ++i )
{
double originalMD = ( *mds )[i];
double originalTVD = ( !tvds->empty() ) ? ( *tvds )[i] : 0.0;
@ -1351,7 +1352,7 @@ std::vector<unsigned char> RigGeoMechWellLogExtractor::determineFilteringOrSmoot
{
std::vector<unsigned char> smoothOrFilterSegments( porePressures.size(), false );
#pragma omp parallel for
for ( int64_t i = 1; i < int64_t( porePressures.size() - 1 ); ++i )
for ( int64_t i = 1; i < static_cast<int64_t>( porePressures.size() - 1 ); ++i )
{
bool validPP_im1 = porePressures[i - 1] >= 0.0 && porePressures[i - 1] != std::numeric_limits<double>::infinity();
bool validPP_i = porePressures[i] >= 0.0 && porePressures[i] != std::numeric_limits<double>::infinity();
@ -1367,7 +1368,7 @@ std::vector<unsigned char> RigGeoMechWellLogExtractor::determineFilteringOrSmoot
//--------------------------------------------------------------------------------------------------
double RigGeoMechWellLogExtractor::hydroStaticPorePressureForIntersection( size_t intersectionIdx, double waterDensityGCM3 ) const
{
double trueVerticalDepth = m_intersectionTVDs[intersectionIdx];
double trueVerticalDepth = cellIntersectionTVDs()[intersectionIdx];
double effectiveDepthMeters = trueVerticalDepth + m_wellPathGeometry->rkbDiff();
return hydroStaticPorePressureAtDepth( effectiveDepthMeters, waterDensityGCM3 );
}
@ -1436,7 +1437,7 @@ bool RigGeoMechWellLogExtractor::isValid( float value )
double RigGeoMechWellLogExtractor::calculateWaterDepth() const
{
// Need a well path with intersections to generate a precise water depth
if ( m_intersectionTVDs.empty() || m_wellPathGeometry->wellPathPoints().empty() )
if ( cellIntersectionTVDs().empty() || m_wellPathGeometry->wellPathPoints().empty() )
{
return std::numeric_limits<double>::infinity();
}
@ -1449,7 +1450,7 @@ double RigGeoMechWellLogExtractor::calculateWaterDepth() const
}
// Water depth is always the first intersection with model for geo mech models.
double waterDepth = m_intersectionTVDs.front();
double waterDepth = cellIntersectionTVDs().front();
return waterDepth;
}

89
ApplicationLibCode/ReservoirDataModel/RigWellLogExtractor.cpp
View File

@ -18,8 +18,10 @@
/////////////////////////////////////////////////////////////////////////////////
#include "RigWellLogExtractor.h"
#include "RiaLogging.h"
#include "RigWellPath.h"
#include "cvfTrace.h"
//--------------------------------------------------------------------------------------------------
@ -85,6 +87,74 @@ std::vector<WellPathCellIntersectionInfo> RigWellLogExtractor::cellIntersectionI
return infoVector;
}
//--------------------------------------------------------------------------------------------------
/// Insert additional intersections if the distance between two intersections is larger than maxDistanceBetweenIntersections
//--------------------------------------------------------------------------------------------------
void RigWellLogExtractor::resampleIntersections( double maxDistanceBetweenIntersections )
{
if ( maxDistanceBetweenIntersections <= 0.0 ) return;
std::vector<cvf::Vec3d> resampledIntersections;
std::vector<size_t> resampledIntersectedCellsGlobIdx;
std::vector<cvf::StructGridInterface::FaceType> resampledIntersectedCellFaces;
std::vector<double> resampledIntersectionMeasuredDepths;
std::vector<double> resampledIntersectionTVDs;
if ( m_intersections.size() > 1 )
{
resampledIntersections.push_back( m_intersections[0] );
resampledIntersectedCellsGlobIdx.push_back( m_intersectedCellsGlobIdx[0] );
resampledIntersectedCellFaces.push_back( m_intersectedCellFaces[0] );
resampledIntersectionMeasuredDepths.push_back( m_intersectionMeasuredDepths[0] );
resampledIntersectionTVDs.push_back( m_intersectionTVDs[0] );
for ( size_t i = 1; i < m_intersections.size(); i++ )
{
const cvf::Vec3d p1 = m_intersections[i - 1];
const cvf::Vec3d p2 = m_intersections[i];
const double distance = ( p2 - p1 ).length();
if ( distance > maxDistanceBetweenIntersections )
{
size_t numberOfPoints = static_cast<size_t>( distance / maxDistanceBetweenIntersections ) + 1;
for ( size_t j = 1; j < numberOfPoints; j++ )
{
// Insert data twice, as the curve visualization algorithm is based on pairs of intersections
resampledIntersectedCellsGlobIdx.push_back( m_intersectedCellsGlobIdx[i] );
resampledIntersectedCellsGlobIdx.push_back( m_intersectedCellsGlobIdx[i] );
resampledIntersectedCellFaces.push_back( m_intersectedCellFaces[i] );
resampledIntersectedCellFaces.push_back( m_intersectedCellFaces[i] );
const double t = static_cast<double>( j ) / static_cast<double>( numberOfPoints );
const cvf::Vec3d p = p1 + ( p2 - p1 ) * t;
resampledIntersections.push_back( p );
resampledIntersections.push_back( p );
const double md = m_intersectionMeasuredDepths[i - 1] +
( m_intersectionMeasuredDepths[i] - m_intersectionMeasuredDepths[i - 1] ) * t;
resampledIntersectionMeasuredDepths.push_back( md );
resampledIntersectionMeasuredDepths.push_back( md );
const double tvd = m_intersectionTVDs[i - 1] + ( m_intersectionTVDs[i] - m_intersectionTVDs[i - 1] ) * t;
resampledIntersectionTVDs.push_back( tvd );
resampledIntersectionTVDs.push_back( tvd );
}
}
resampledIntersections.push_back( p2 );
resampledIntersectedCellsGlobIdx.push_back( m_intersectedCellsGlobIdx[i] );
resampledIntersectedCellFaces.push_back( m_intersectedCellFaces[i] );
resampledIntersectionMeasuredDepths.push_back( m_intersectionMeasuredDepths[i] );
resampledIntersectionTVDs.push_back( m_intersectionTVDs[i] );
}
}
m_intersectionMeasuredDepths = resampledIntersectionMeasuredDepths;
m_intersectionTVDs = resampledIntersectionTVDs;
m_intersections = resampledIntersections;
m_intersectedCellsGlobIdx = resampledIntersectedCellsGlobIdx;
m_intersectedCellFaces = resampledIntersectedCellFaces;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
@ -93,6 +163,22 @@ const std::vector<size_t>& RigWellLogExtractor::intersectedCellsGlobIdx() const
return m_intersectedCellsGlobIdx;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
const std::vector<cvf::Vec3d>& RigWellLogExtractor::intersections() const
{
return m_intersections;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
const std::vector<cvf::StructGridInterface::FaceType>& RigWellLogExtractor::intersectedCellFaces() const
{
return m_intersectedCellFaces;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
@ -319,6 +405,9 @@ void RigWellLogExtractor::populateReturnArrays( std::map<RigMDCellIdxEnterLeaveK
}
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RigWellLogExtractor::appendIntersectionToArrays( double measuredDepth,
const HexIntersectionInfo& intersection,
gsl::not_null<QStringList*> errorMessages )

24
ApplicationLibCode/ReservoirDataModel/RigWellLogExtractor.h
View File

@ -61,14 +61,18 @@ public:
RigWellLogExtractor( gsl::not_null<const RigWellPath*> wellpath, const std::string& wellCaseErrorMsgName );
~RigWellLogExtractor() override;
const std::vector<double>& cellIntersectionMDs() const;
const std::vector<double>& cellIntersectionTVDs() const;
const std::vector<size_t>& intersectedCellsGlobIdx() const;
const std::vector<double>& cellIntersectionMDs() const;
const std::vector<double>& cellIntersectionTVDs() const;
const std::vector<size_t>& intersectedCellsGlobIdx() const;
const std::vector<cvf::Vec3d>& intersections() const;
const std::vector<cvf::StructGridInterface::FaceType>& intersectedCellFaces() const;
const RigWellPath* wellPathGeometry() const;
std::vector<WellPathCellIntersectionInfo> cellIntersectionInfosAlongWellPath() const;
void resampleIntersections( double maxDistanceBetweenIntersections );
protected:
static void insertIntersectionsInMap( const std::vector<HexIntersectionInfo>& intersections,
cvf::Vec3d p1,
@ -83,14 +87,14 @@ protected:
virtual cvf::Vec3d calculateLengthInCell( size_t cellIndex, const cvf::Vec3d& startPoint, const cvf::Vec3d& endPoint ) const = 0;
protected:
cvf::cref<RigWellPath> m_wellPathGeometry;
std::string m_wellCaseErrorMsgName;
private:
std::vector<cvf::Vec3d> m_intersections;
std::vector<size_t> m_intersectedCellsGlobIdx;
std::vector<cvf::StructGridInterface::FaceType> m_intersectedCellFaces;
cvf::cref<RigWellPath> m_wellPathGeometry;
std::vector<double> m_intersectionMeasuredDepths;
std::vector<double> m_intersectionTVDs;
std::string m_wellCaseErrorMsgName;
std::vector<double> m_intersectionMeasuredDepths;
std::vector<double> m_intersectionTVDs;
};