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

#6031 Use reference case for Pore Compressibility calculation.

Browse Source
This commit is contained in:
Kristian Bendiksen 2020-06-18 09:28:04 +02:00
parent 034f75f64e
commit d17faca18a
5 changed files with 147 additions and 28 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

40
ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultCalculatorPoreCompressibility.cpp
View File

@ -140,16 +140,20 @@ RigFemScalarResultFrames*
frameCountProgress.setNextProgressIncrement( 1u );
int referenceFrameIdx = m_resultCollection->referenceTimeStep();
int frameCount = srcEVDataFrames->frameCount();
for ( int fIdx = 0; fIdx < frameCount; ++fIdx )
{
const std::vector<float>& evData = srcEVDataFrames->frameData( fIdx );
const std::vector<float>& verticalStrainData = verticalStrainDataFrames->frameData( fIdx );
const std::vector<float>& youngsModuliData = youngsModuliFrames->frameData( fIdx );
const std::vector<float>& poissonRatioData = poissonRatioFrames->frameData( fIdx );
const std::vector<float>& voidRatioData = voidRatioFrames->frameData( fIdx );
const std::vector<float>& initialPorFrameData = srcPORDataFrames->frameData( 0 );
const std::vector<float>& porFrameData = srcPORDataFrames->frameData( fIdx );
const std::vector<float>& evData = srcEVDataFrames->frameData( fIdx );
const std::vector<float>& referenceEvData = srcEVDataFrames->frameData( referenceFrameIdx );
const std::vector<float>& verticalStrainData = verticalStrainDataFrames->frameData( fIdx );
const std::vector<float>& referenceVerticalStrainData = verticalStrainDataFrames->frameData( referenceFrameIdx );
const std::vector<float>& youngsModuliData = youngsModuliFrames->frameData( fIdx );
const std::vector<float>& poissonRatioData = poissonRatioFrames->frameData( fIdx );
const std::vector<float>& voidRatioData = voidRatioFrames->frameData( referenceFrameIdx );
const std::vector<float>& referencePorFrameData = srcPORDataFrames->frameData( referenceFrameIdx );
const std::vector<float>& porFrameData = srcPORDataFrames->frameData( fIdx );
std::vector<float>& poreCompressibilityFrameData = poreCompressibilityFrames->frameData( fIdx );
std::vector<float>& verticalCompressibilityFrameData = verticalCompressibilityFrames->frameData( fIdx );
@ -187,16 +191,16 @@ RigFemScalarResultFrames*
size_t elmNodResIdx = femPart->elementNodeResultIdx( elmIdx, elmNodIdx );
if ( elmNodResIdx < evData.size() )
{
if ( fIdx == 0 )
if ( fIdx == referenceFrameIdx )
{
// Geostatic step: result not defined
// The time step and the reference time step are the same: results undefined
poreCompressibilityFrameData[elmNodResIdx] = inf;
verticalCompressibilityFrameData[elmNodResIdx] = inf;
verticalCompressibilityRatioFrameData[elmNodResIdx] = inf;
}
else
{
// Use biot coefficient for all other (not Geostatic) timesteps
// Use biot coefficient for all timesteps
double biotCoefficient = 1.0;
if ( biotData.empty() )
{
@ -221,16 +225,16 @@ RigFemScalarResultFrames*
double porosity = voidr / ( 1.0 + voidr );
// Calculate difference in pore pressure between reference state and this state
double initialPorePressure = initialPorFrameData[nodeIdx];
double framePorePressure = porFrameData[nodeIdx];
double deltaPorePressure = framePorePressure - initialPorePressure;
double referencePorePressure = referencePorFrameData[nodeIdx];
double framePorePressure = porFrameData[nodeIdx];
double deltaPorePressure = framePorePressure - referencePorePressure;
// Calculate pore compressibility
double poreCompressibility = inf;
if ( deltaPorePressure != 0.0 && porosity != 0.0 )
{
poreCompressibility =
-( biotCoefficient * evData[elmNodResIdx] ) / ( deltaPorePressure * porosity );
double deltaEv = evData[elmNodResIdx] - referenceEvData[elmNodResIdx];
poreCompressibility = -( biotCoefficient * deltaEv ) / ( deltaPorePressure * porosity );
// Guard against divide by zero: second term can be ignored when bulk modulus is zero,
// which can happens when biot coefficient is 1.0
if ( biotCoefficient != 1.0 && porosity != 1.0 )
@ -245,9 +249,11 @@ RigFemScalarResultFrames*
if ( biotCoefficient != 0.0 && deltaPorePressure != 0.0 )
{
double deltaStrain = verticalStrainData[elmNodResIdx] -
referenceVerticalStrainData[elmNodResIdx];
// Calculate vertical compressibility
verticalCompressibility =
-verticalStrainData[elmNodResIdx] / ( biotCoefficient * deltaPorePressure );
verticalCompressibility = -deltaStrain / ( biotCoefficient * deltaPorePressure );
// Calculate vertical compressibility ratio
verticalCompressibilityRatio =

76
ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.cpp
View File

@ -111,6 +111,8 @@ RigFemPartResultsCollection::RigFemPartResultsCollection( RifGeoMechReaderInterf
m_biotFixedFactor = 1.0;
m_biotResultAddress = "";
m_referenceTimeStep = 0;
m_resultCalculators.push_back(
std::unique_ptr<RigFemPartResultCalculator>( new RigFemPartResultCalculatorTimeLapse( *this ) ) );
m_resultCalculators.push_back(
@ -395,6 +397,28 @@ void RigFemPartResultsCollection::setBiotCoefficientParameters( double biotFixed
}
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RigFemPartResultsCollection::setReferenceTimeStep( int referenceTimeStep )
{
m_referenceTimeStep = referenceTimeStep;
std::set<RigFemResultAddress> results = referenceCaseDependentResults();
for ( auto result : results )
{
deleteResult( result );
}
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
int RigFemPartResultsCollection::referenceTimeStep() const
{
return m_referenceTimeStep;
}
//--------------------------------------------------------------------------------------------------
/// Will always return a valid object, but it can be empty
//--------------------------------------------------------------------------------------------------
@ -1150,6 +1174,23 @@ std::vector<RigFemResultAddress>
return addresses;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RigFemPartResultsCollection::isResultInSet( const RigFemResultAddress& result,
const std::set<RigFemResultAddress>& results )
{
for ( auto res : results )
{
if ( res.resultPosType == result.resultPosType && res.fieldName == result.fieldName &&
res.componentName == result.componentName )
{
return true;
}
}
return false;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
@ -1182,15 +1223,36 @@ std::set<RigFemResultAddress> RigFemPartResultsCollection::normalizedResults()
//--------------------------------------------------------------------------------------------------
bool RigFemPartResultsCollection::isNormalizableResult( const RigFemResultAddress& result )
{
for ( auto normRes : normalizedResults() )
return isResultInSet( result, normalizedResults() );
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::set<RigFemResultAddress> RigFemPartResultsCollection::referenceCaseDependentResults()
{
std::set<RigFemResultAddress> results;
for ( auto elementType : {RIG_ELEMENT_NODAL, RIG_INTEGRATION_POINT} )
{
if ( normRes.resultPosType == result.resultPosType && normRes.fieldName == result.fieldName &&
normRes.componentName == result.componentName )
{
return true;
}
results.insert(
RigFemResultAddress( elementType, "COMPRESSIBILITY", "PORE", RigFemResultAddress::allTimeLapsesValue() ) );
results.insert(
RigFemResultAddress( elementType, "COMPRESSIBILITY", "VERTICAL", RigFemResultAddress::allTimeLapsesValue() ) );
results.insert( RigFemResultAddress( elementType,
"COMPRESSIBILITY",
"VERTICAL-RATIO",
RigFemResultAddress::allTimeLapsesValue() ) );
}
return false;
return results;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RigFemPartResultsCollection::isReferenceCaseDependentResult( const RigFemResultAddress& result )
{
return isResultInSet( result, referenceCaseDependentResults() );
}
//--------------------------------------------------------------------------------------------------

10
ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.h
View File

@ -120,6 +120,8 @@ public:
double* globalPosClosestToZero,
double* globalNegClosestToZero );
static bool isResultInSet( const RigFemResultAddress& result, const std::set<RigFemResultAddress>& results );
static std::vector<RigFemResultAddress> tensorComponentAddresses( const RigFemResultAddress& resVarAddr );
static std::vector<RigFemResultAddress> tensorPrincipalComponentAdresses( const RigFemResultAddress& resVarAddr );
static std::set<RigFemResultAddress> normalizedResults();
@ -128,6 +130,12 @@ public:
void setNormalizationAirGap( double normalizationAirGap );
double normalizationAirGap() const;
void setReferenceTimeStep( int referenceTimeStep );
int referenceTimeStep() const;
static std::set<RigFemResultAddress> referenceCaseDependentResults();
static bool isReferenceCaseDependentResult( const RigFemResultAddress& result );
RigFemScalarResultFrames* findOrLoadScalarResult( int partIndex, const RigFemResultAddress& resVarAddr );
RigFemScalarResultFrames* createScalarResult( int partIndex, const RigFemResultAddress& resVarAddr );
@ -154,6 +162,8 @@ private:
double m_biotFixedFactor;
QString m_biotResultAddress;
int m_referenceTimeStep;
std::vector<std::unique_ptr<RigFemPartResultCalculator>> m_resultCalculators;
RigStatisticsDataCache* statistics( const RigFemResultAddress& resVarAddr );

47
ApplicationCode/ProjectDataModel/RimGeoMechResultDefinition.cpp
View File

@ -93,6 +93,7 @@ RimGeoMechResultDefinition::RimGeoMechResultDefinition( void )
"",
"",
"" );
CAF_PDM_InitField( &m_referenceTimeStep, "ReferenceTimeStep", 0, "Reference Time Step", "", "", "" );
CAF_PDM_InitField( &m_compactionRefLayer, "CompactionRefLayer", 0, "Compaction Ref Layer", "", "", "" );
m_compactionRefLayer.uiCapability()->setUiHidden( true );
@ -170,8 +171,17 @@ void RimGeoMechResultDefinition::defineUiOrdering( QString uiConfigName, caf::Pd
if ( m_resultPositionTypeUiField() != RIG_FORMATION_NAMES )
{
caf::PdmUiGroup* timeLapseGr = uiOrdering.addNewGroup( "Difference Options" );
timeLapseGr->add( &m_timeLapseBaseTimestep );
bool isReferenceCaseDependent = referenceCaseDependentResultSelected();
if ( isReferenceCaseDependent )
{
caf::PdmUiGroup* referenceCaseGr = uiOrdering.addNewGroup( "Reference Case" );
referenceCaseGr->add( &m_referenceTimeStep );
}
else
{
caf::PdmUiGroup* timeLapseGr = uiOrdering.addNewGroup( "Difference Options" );
timeLapseGr->add( &m_timeLapseBaseTimestep );
}
}
if ( m_resultPositionTypeUiField() == RIG_NODAL )
@ -264,6 +274,21 @@ QList<caf::PdmOptionItemInfo>
static_cast<int>( stepIdx ) ) );
}
}
else if ( &m_referenceTimeStep == fieldNeedingOptions )
{
std::vector<std::string> stepNames;
if ( m_geomCase->geoMechData() )
{
stepNames = m_geomCase->geoMechData()->femPartResults()->filteredStepNames();
}
for ( size_t stepIdx = 0; stepIdx < stepNames.size(); ++stepIdx )
{
options.push_back(
caf::PdmOptionItemInfo( QString( "%1 (#%2)" ).arg( QString::fromStdString( stepNames[stepIdx] ) ).arg( stepIdx ),
static_cast<int>( stepIdx ) ) );
}
}
else if ( &m_compactionRefLayerUiField == fieldNeedingOptions )
{
if ( m_geomCase->geoMechData() )
@ -319,7 +344,8 @@ void RimGeoMechResultDefinition::fieldChangedByUi( const caf::PdmFieldHandle* ch
}
}
if ( &m_resultPositionTypeUiField == changedField || &m_timeLapseBaseTimestep == changedField )
if ( &m_resultPositionTypeUiField == changedField || &m_timeLapseBaseTimestep == changedField ||
&m_referenceTimeStep == changedField )
{
std::map<std::string, std::vector<std::string>> fieldCompNames = getResultMetaDataForUIFieldSetting();
QStringList uiVarNames;
@ -336,6 +362,11 @@ void RimGeoMechResultDefinition::fieldChangedByUi( const caf::PdmFieldHandle* ch
{
m_resultVariableUiField = "";
}
if ( &m_referenceTimeStep == changedField )
{
m_geomCase->geoMechData()->femPartResults()->setReferenceTimeStep( m_referenceTimeStep() );
}
}
if ( &m_normalizeByHydrostaticPressure == changedField && m_normalizationAirGap == 0.0 )
{
@ -357,7 +388,7 @@ void RimGeoMechResultDefinition::fieldChangedByUi( const caf::PdmFieldHandle* ch
if ( &m_resultVariableUiField == changedField || &m_compactionRefLayerUiField == changedField ||
&m_timeLapseBaseTimestep == changedField || &m_normalizeByHydrostaticPressure == changedField ||
&m_normalizationAirGap == changedField )
&m_normalizationAirGap == changedField || &m_referenceTimeStep == changedField )
{
QStringList fieldComponentNames = m_resultVariableUiField().split( QRegExp( "\\s+" ) );
if ( fieldComponentNames.size() > 0 )
@ -854,6 +885,14 @@ bool RimGeoMechResultDefinition::normalizableResultSelected() const
return RigFemPartResultsCollection::isNormalizableResult( this->resultAddress() );
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RimGeoMechResultDefinition::referenceCaseDependentResultSelected() const
{
return RigFemPartResultsCollection::isReferenceCaseDependentResult( this->resultAddress() );
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------

2
ApplicationCode/ProjectDataModel/RimGeoMechResultDefinition.h
View File

@ -112,6 +112,7 @@ private:
static QString convertToUiResultFieldName( QString resultFieldName );
bool normalizableResultSelected() const;
bool referenceCaseDependentResultSelected() const;
// Data Fields
@ -122,6 +123,7 @@ private:
caf::PdmField<int> m_compactionRefLayer;
caf::PdmField<bool> m_normalizeByHydrostaticPressure;
caf::PdmField<double> m_normalizationAirGap;
caf::PdmField<int> m_referenceTimeStep;
// UI Fields only
friend class RimGeoMechPropertyFilter; // Property filter needs the ui fields