diff --git a/ApplicationCode/GeoMech/GeoMechDataModel/CMakeLists.txt b/ApplicationCode/GeoMech/GeoMechDataModel/CMakeLists.txt index a0313dc6d4..943026882a 100644 --- a/ApplicationCode/GeoMech/GeoMechDataModel/CMakeLists.txt +++ b/ApplicationCode/GeoMech/GeoMechDataModel/CMakeLists.txt @@ -94,6 +94,8 @@ add_library( ${PROJECT_NAME} RigFemPartResultCalculatorNodalGradients.cpp RigFemPartResultCalculatorStressAnisotropy.h RigFemPartResultCalculatorStressAnisotropy.cpp + RigFemPartResultCalculatorPoreCompressibility.h + RigFemPartResultCalculatorPoreCompressibility.cpp RimGeoMechGeometrySelectionItem.h RimGeoMechGeometrySelectionItem.cpp ) diff --git a/ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultCalculatorPoreCompressibility.cpp b/ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultCalculatorPoreCompressibility.cpp new file mode 100644 index 0000000000..65a055940c --- /dev/null +++ b/ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultCalculatorPoreCompressibility.cpp @@ -0,0 +1,282 @@ +///////////////////////////////////////////////////////////////////////////////// +// +// 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 +// for more details. +// +///////////////////////////////////////////////////////////////////////////////// + +#include "RigFemPartResultCalculatorPoreCompressibility.h" + +#include "RiaLogging.h" + +#include "RigFemPart.h" +#include "RigFemPartCollection.h" +#include "RigFemPartResultsCollection.h" +#include "RigFemResultAddress.h" +#include "RigFemScalarResultFrames.h" + +#include "cafProgressInfo.h" + +#include + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +RigFemPartResultCalculatorPoreCompressibility::RigFemPartResultCalculatorPoreCompressibility( RigFemPartResultsCollection& collection ) + : RigFemPartResultCalculator( collection ) +{ +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +RigFemPartResultCalculatorPoreCompressibility::~RigFemPartResultCalculatorPoreCompressibility() +{ +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +bool RigFemPartResultCalculatorPoreCompressibility::isMatching( const RigFemResultAddress& resVarAddr ) const +{ + return ( resVarAddr.fieldName == "PORE-COMPRESSIBILITY" && + ( resVarAddr.componentName == "PORE" || resVarAddr.componentName == "VERTICAL" || + resVarAddr.componentName == "VERTICAL-RATIO" ) ); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +RigFemScalarResultFrames* + RigFemPartResultCalculatorPoreCompressibility::calculate( int partIndex, const RigFemResultAddress& resVarAddr ) +{ + caf::ProgressInfo frameCountProgress( m_resultCollection->frameCount() * 6, "" ); + frameCountProgress.setProgressDescription( "Calculating Pore Compressibility" ); + + frameCountProgress.setNextProgressIncrement( m_resultCollection->frameCount() ); + RigFemScalarResultFrames* srcPORDataFrames = + m_resultCollection->findOrLoadScalarResult( partIndex, RigFemResultAddress( RIG_NODAL, "POR-Bar", "" ) ); + frameCountProgress.incrementProgress(); + + // Volumetric Strain + frameCountProgress.setNextProgressIncrement( m_resultCollection->frameCount() ); + RigFemScalarResultFrames* srcEVDataFrames = + m_resultCollection->findOrLoadScalarResult( partIndex, RigFemResultAddress( resVarAddr.resultPosType, "NE", "EV" ) ); + + frameCountProgress.incrementProgress(); + + // Vertical Strain + frameCountProgress.setNextProgressIncrement( m_resultCollection->frameCount() ); + RigFemScalarResultFrames* verticalStrainDataFrames = + m_resultCollection->findOrLoadScalarResult( partIndex, + RigFemResultAddress( resVarAddr.resultPosType, "NE", "E33" ) ); + + frameCountProgress.incrementProgress(); + + // Biot porelastic coeffisient (alpha) + frameCountProgress.setNextProgressIncrement( m_resultCollection->frameCount() ); + RigFemScalarResultFrames* biotCoefficient = nullptr; + if ( !m_resultCollection->biotResultAddress().isEmpty() ) + { + biotCoefficient = + m_resultCollection + ->findOrLoadScalarResult( partIndex, + RigFemResultAddress( RIG_ELEMENT, + m_resultCollection->biotResultAddress().toStdString(), + "" ) ); + } + frameCountProgress.incrementProgress(); + + frameCountProgress.setNextProgressIncrement( m_resultCollection->frameCount() ); + RigFemScalarResultFrames* youngsModuliFrames = + m_resultCollection->findOrLoadScalarResult( partIndex, RigFemResultAddress( RIG_ELEMENT, "MODULUS", "" ) ); + if ( youngsModuliFrames->frameData( 0 ).empty() ) + { + RiaLogging::error( "Missing Youngs Moduli element data (MODULUS)." ); + return nullptr; + } + + RigFemScalarResultFrames* poissonRatioFrames = + m_resultCollection->findOrLoadScalarResult( partIndex, RigFemResultAddress( RIG_ELEMENT, "RATIO", "" ) ); + if ( youngsModuliFrames->frameData( 0 ).empty() ) + { + RiaLogging::error( "Missing Poisson Ratio element data (RATIO)." ); + return nullptr; + } + + RigFemScalarResultFrames* voidRatioFrames = + m_resultCollection->findOrLoadScalarResult( partIndex, + RigFemResultAddress( resVarAddr.resultPosType, "VOIDR", "" ) ); + + RigFemScalarResultFrames* poreCompressibilityFrames = + m_resultCollection->createScalarResult( partIndex, + RigFemResultAddress( resVarAddr.resultPosType, resVarAddr.fieldName, "PORE" ) ); + + RigFemScalarResultFrames* verticalCompressibilityFrames = + m_resultCollection->createScalarResult( partIndex, + RigFemResultAddress( resVarAddr.resultPosType, + resVarAddr.fieldName, + "VERTICAL" ) ); + RigFemScalarResultFrames* verticalCompressibilityRatioFrames = + m_resultCollection->createScalarResult( partIndex, + RigFemResultAddress( resVarAddr.resultPosType, + resVarAddr.fieldName, + "VERTICAL-RATIO" ) ); + frameCountProgress.incrementProgress(); + + const RigFemPart* femPart = m_resultCollection->parts()->part( partIndex ); + float inf = std::numeric_limits::infinity(); + + frameCountProgress.setNextProgressIncrement( 1u ); + + int frameCount = srcEVDataFrames->frameCount(); + for ( int fIdx = 0; fIdx < frameCount; ++fIdx ) + { + const std::vector& evData = srcEVDataFrames->frameData( fIdx ); + const std::vector& verticalStrainData = verticalStrainDataFrames->frameData( fIdx ); + const std::vector& youngsModuliData = youngsModuliFrames->frameData( fIdx ); + const std::vector& poissonRatioData = poissonRatioFrames->frameData( fIdx ); + const std::vector& voidRatioData = voidRatioFrames->frameData( fIdx ); + const std::vector& initialPorFrameData = srcPORDataFrames->frameData( 0 ); + const std::vector& porFrameData = srcPORDataFrames->frameData( fIdx ); + + std::vector& poreCompressibilityFrameData = poreCompressibilityFrames->frameData( fIdx ); + std::vector& verticalCompressibilityFrameData = verticalCompressibilityFrames->frameData( fIdx ); + std::vector& verticalCompressibilityRatioFrameData = verticalCompressibilityRatioFrames->frameData( fIdx ); + + size_t valCount = evData.size(); + poreCompressibilityFrameData.resize( valCount ); + verticalCompressibilityFrameData.resize( valCount ); + verticalCompressibilityRatioFrameData.resize( valCount ); + + int elementCount = femPart->elementCount(); + + std::vector biotData; + if ( biotCoefficient ) + { + biotData = biotCoefficient->frameData( fIdx ); + if ( !m_resultCollection->isValidBiotData( biotData, elementCount ) ) + { + m_resultCollection->deleteResult( resVarAddr ); + return nullptr; + } + } + +#pragma omp parallel for + for ( int elmIdx = 0; elmIdx < elementCount; ++elmIdx ) + { + RigElementType elmType = femPart->elementType( elmIdx ); + + int elmNodeCount = RigFemTypes::elmentNodeCount( femPart->elementType( elmIdx ) ); + + if ( elmType == HEX8P ) + { + for ( int elmNodIdx = 0; elmNodIdx < elmNodeCount; ++elmNodIdx ) + { + size_t elmNodResIdx = femPart->elementNodeResultIdx( elmIdx, elmNodIdx ); + if ( elmNodResIdx < evData.size() ) + { + if ( fIdx == 0 ) + { + // Geostatic step: result not defined + poreCompressibilityFrameData[elmNodResIdx] = inf; + verticalCompressibilityFrameData[elmNodResIdx] = inf; + verticalCompressibilityRatioFrameData[elmNodResIdx] = inf; + } + else + { + // Use biot coefficient for all other (not Geostatic) timesteps + double biotCoefficient = 1.0; + if ( biotData.empty() ) + { + biotCoefficient = m_resultCollection->biotFixedFactor(); + } + else + { + // Use coefficient from element property table + biotCoefficient = biotData[elmIdx]; + } + + int nodeIdx = femPart->nodeIdxFromElementNodeResultIdx( elmNodResIdx ); + + // Calculate bulk modulus for solids (grains) + double poissonRatio = poissonRatioData[elmIdx]; + double youngsModuli = youngsModuliData[elmIdx]; + double bulkModulusFrame = youngsModuli / ( 3.0 * ( 1.0 - 2.0 * poissonRatio ) ); + double bulkModulus = bulkModulusFrame / ( 1.0 - biotCoefficient ); + + // Calculate porosity + double voidr = voidRatioData[elmNodResIdx]; + 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; + + // Calculate pore compressibility + double poreCompressibility = inf; + if ( deltaPorePressure != 0.0 && porosity != 0.0 ) + { + poreCompressibility = + ( biotCoefficient * evData[elmNodResIdx] ) / ( 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 ) + { + poreCompressibility += ( 1.0 / bulkModulus ) * ( biotCoefficient / porosity - 1.0 ); + } + } + poreCompressibilityFrameData[elmNodResIdx] = poreCompressibility; + + double verticalCompressibility = inf; + double verticalCompressibilityRatio = inf; + + if ( biotCoefficient != 0.0 && deltaPorePressure != 0.0 ) + { + // Calculate vertical compressibility + verticalCompressibility = + -verticalStrainData[elmNodResIdx] / ( biotCoefficient * deltaPorePressure ); + + // Calculate vertical compressibility ratio + verticalCompressibilityRatio = + ( verticalCompressibility * youngsModuli * ( 1.0 - poissonRatio ) ) / + ( ( 1.0 + poissonRatio ) * ( 1.0 - 2.0 * poissonRatio ) ); + } + + verticalCompressibilityFrameData[elmNodResIdx] = verticalCompressibility; + verticalCompressibilityRatioFrameData[elmNodResIdx] = verticalCompressibilityRatio; + } + } + } + } + else + { + for ( int elmNodIdx = 0; elmNodIdx < elmNodeCount; ++elmNodIdx ) + { + size_t elmNodResIdx = femPart->elementNodeResultIdx( elmIdx, elmNodIdx ); + if ( elmNodResIdx < poreCompressibilityFrameData.size() ) + { + poreCompressibilityFrameData[elmNodResIdx] = inf; + } + } + } + } + + frameCountProgress.incrementProgress(); + } + + RigFemScalarResultFrames* requestedResultFrames = m_resultCollection->findOrLoadScalarResult( partIndex, resVarAddr ); + return requestedResultFrames; +} diff --git a/ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultCalculatorPoreCompressibility.h b/ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultCalculatorPoreCompressibility.h new file mode 100644 index 0000000000..e4d433367e --- /dev/null +++ b/ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultCalculatorPoreCompressibility.h @@ -0,0 +1,37 @@ +///////////////////////////////////////////////////////////////////////////////// +// +// 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 +// for more details. +// +///////////////////////////////////////////////////////////////////////////////// + +#pragma once + +#include "RigFemPartResultCalculator.h" + +class RigFemPartResultsCollection; +class RigFemScalarResultFrames; +class RigFemResultAddress; + +//================================================================================================== +/// +//================================================================================================== +class RigFemPartResultCalculatorPoreCompressibility : public RigFemPartResultCalculator +{ +public: + explicit RigFemPartResultCalculatorPoreCompressibility( RigFemPartResultsCollection& collection ); + virtual ~RigFemPartResultCalculatorPoreCompressibility(); + bool isMatching( const RigFemResultAddress& resVarAddr ) const override; + RigFemScalarResultFrames* calculate( int partIndex, const RigFemResultAddress& resVarAddr ) override; +}; diff --git a/ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.cpp b/ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.cpp index 1b9051d2ba..37b74c1322 100644 --- a/ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.cpp +++ b/ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.cpp @@ -23,6 +23,7 @@ #include "RifElementPropertyReader.h" #include "RifGeoMechReaderInterface.h" +#include "RigFemPartResultCalculatorPoreCompressibility.h" #ifdef USE_ODB_API #include "RifOdbReader.h" @@ -76,7 +77,7 @@ #include #include -#include +#include //-------------------------------------------------------------------------------------------------- /// @@ -164,6 +165,8 @@ RigFemPartResultsCollection::RigFemPartResultsCollection( RifGeoMechReaderInterf std::unique_ptr( new RigFemPartResultCalculatorPrincipalStress( *this ) ) ); m_resultCalculators.push_back( std::unique_ptr( new RigFemPartResultCalculatorStressAnisotropy( *this ) ) ); + m_resultCalculators.push_back( + std::unique_ptr( new RigFemPartResultCalculatorPoreCompressibility( *this ) ) ); m_resultCalculators.push_back( std::unique_ptr( new RigFemPartResultCalculatorFormationIndices( *this ) ) ); } @@ -351,6 +354,17 @@ void RigFemPartResultsCollection::setBiotCoefficientParameters( double biotFixed } } + deleteResult( + RigFemResultAddress( elementType, "PORE-COMPRESSIBILITY", "PORE", RigFemResultAddress::allTimeLapsesValue() ) ); + deleteResult( RigFemResultAddress( elementType, + "PORE-COMPRESSIBILITY", + "VERTICAL", + RigFemResultAddress::allTimeLapsesValue() ) ); + deleteResult( RigFemResultAddress( elementType, + "PORE-COMPRESSIBILITY", + "VERTICAL-RATIO", + RigFemResultAddress::allTimeLapsesValue() ) ); + // SE only: depends on SE.S1 and SE.S3 deleteResult( RigFemResultAddress( elementType, "SE", "SFI", RigFemResultAddress::allTimeLapsesValue() ) ); deleteResult( RigFemResultAddress( elementType, "SE", "DSM", RigFemResultAddress::allTimeLapsesValue() ) ); @@ -600,6 +614,10 @@ std::map> fieldCompNames["NE"].push_back( "E1" ); fieldCompNames["NE"].push_back( "E2" ); fieldCompNames["NE"].push_back( "E3" ); + + fieldCompNames["PORE-COMPRESSIBILITY"].push_back( "PORE" ); + fieldCompNames["PORE-COMPRESSIBILITY"].push_back( "VERTICAL" ); + fieldCompNames["PORE-COMPRESSIBILITY"].push_back( "VERTICAL-RATIO" ); } else if ( resPos == RIG_INTEGRATION_POINT ) { @@ -675,6 +693,10 @@ std::map> fieldCompNames["NE"].push_back( "E1" ); fieldCompNames["NE"].push_back( "E2" ); fieldCompNames["NE"].push_back( "E3" ); + + fieldCompNames["PORE-COMPRESSIBILITY"].push_back( "PORE" ); + fieldCompNames["PORE-COMPRESSIBILITY"].push_back( "VERTICAL" ); + fieldCompNames["PORE-COMPRESSIBILITY"].push_back( "VERTICAL-RATIO" ); } else if ( resPos == RIG_ELEMENT_NODAL_FACE ) {