diff --git a/ApplicationCode/GeoMech/GeoMechDataModel/CMakeLists.txt b/ApplicationCode/GeoMech/GeoMechDataModel/CMakeLists.txt index 5dd34d35d7..4adf7f2cee 100644 --- a/ApplicationCode/GeoMech/GeoMechDataModel/CMakeLists.txt +++ b/ApplicationCode/GeoMech/GeoMechDataModel/CMakeLists.txt @@ -96,6 +96,8 @@ add_library( ${PROJECT_NAME} RigFemPartResultCalculatorPoreCompressibility.cpp RigFemPartResultCalculatorPorosityPermeability.h RigFemPartResultCalculatorPorosityPermeability.cpp + RigFemPartResultCalculatorInitialPorosity.h + RigFemPartResultCalculatorInitialPorosity.cpp RigFemPartResultCalculatorMudWeightWindow.h RigFemPartResultCalculatorMudWeightWindow.cpp RigFemPartResultCalculatorShearSlipIndicator.h diff --git a/ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultCalculatorInitialPorosity.cpp b/ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultCalculatorInitialPorosity.cpp new file mode 100644 index 0000000000..e24bffd11d --- /dev/null +++ b/ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultCalculatorInitialPorosity.cpp @@ -0,0 +1,137 @@ +///////////////////////////////////////////////////////////////////////////////// +// +// 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 "RigFemPartResultCalculatorInitialPorosity.h" + +#include "RiaEclipseUnitTools.h" +#include "RiaLogging.h" + +#include "RigFemPart.h" +#include "RigFemPartCollection.h" +#include "RigFemPartResultsCollection.h" +#include "RigFemResultAddress.h" +#include "RigFemScalarResultFrames.h" + +#include "cafProgressInfo.h" + +#include + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +RigFemPartResultCalculatorInitialPorosity::RigFemPartResultCalculatorInitialPorosity( RigFemPartResultsCollection& collection ) + : RigFemPartResultCalculator( collection ) +{ +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +RigFemPartResultCalculatorInitialPorosity::~RigFemPartResultCalculatorInitialPorosity() +{ +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +bool RigFemPartResultCalculatorInitialPorosity::isMatching( const RigFemResultAddress& resVarAddr ) const +{ + return ( resVarAddr.fieldName == "PORO-PERM" && resVarAddr.componentName == "PHI0" ); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +RigFemScalarResultFrames* RigFemPartResultCalculatorInitialPorosity::calculate( int partIndex, + const RigFemResultAddress& resVarAddr ) +{ + caf::ProgressInfo frameCountProgress( m_resultCollection->frameCount() * 2, "" ); + frameCountProgress.setProgressDescription( "Calculating Initial Porosity" ); + + frameCountProgress.setNextProgressIncrement( m_resultCollection->frameCount() ); + + RigFemScalarResultFrames* voidRatioFrames = + m_resultCollection->findOrLoadScalarResult( partIndex, + RigFemResultAddress( resVarAddr.resultPosType, "VOIDR", "" ) ); + + RigFemScalarResultFrames* porosityFrames = + m_resultCollection->createScalarResult( partIndex, + RigFemResultAddress( resVarAddr.resultPosType, resVarAddr.fieldName, "PHI0" ) ); + frameCountProgress.incrementProgress(); + + const RigFemPart* femPart = m_resultCollection->parts()->part( partIndex ); + float inf = std::numeric_limits::infinity(); + + frameCountProgress.setNextProgressIncrement( 1u ); + + int referenceFrameIdx = m_resultCollection->referenceTimeStep(); + + int frameCount = voidRatioFrames->frameCount(); + for ( int fIdx = 0; fIdx < frameCount; ++fIdx ) + { + const std::vector& voidRatioData = voidRatioFrames->frameData( 0 ); + std::vector& porosityFrameData = porosityFrames->frameData( fIdx ); + + size_t valCount = voidRatioData.size(); + porosityFrameData.resize( valCount ); + + int elementCount = femPart->elementCount(); + +#pragma omp parallel for + for ( int elmIdx = 0; elmIdx < elementCount; ++elmIdx ) + { + RigElementType elmType = femPart->elementType( elmIdx ); + + int elmNodeCount = RigFemTypes::elementNodeCount( femPart->elementType( elmIdx ) ); + + if ( elmType == HEX8P ) + { + for ( int elmNodIdx = 0; elmNodIdx < elmNodeCount; ++elmNodIdx ) + { + size_t elmNodResIdx = femPart->elementNodeResultIdx( elmIdx, elmNodIdx ); + if ( elmNodResIdx < voidRatioData.size() ) + { + int nodeIdx = femPart->nodeIdxFromElementNodeResultIdx( elmNodResIdx ); + + // Calculate initial porosity + double voidr = voidRatioData[elmNodResIdx]; + double initialPorosity = voidr / ( 1.0 + voidr ); + + porosityFrameData[elmNodResIdx] = initialPorosity; + } + } + } + else + { + for ( int elmNodIdx = 0; elmNodIdx < elmNodeCount; ++elmNodIdx ) + { + size_t elmNodResIdx = femPart->elementNodeResultIdx( elmIdx, elmNodIdx ); + if ( elmNodResIdx < voidRatioData.size() ) + { + porosityFrameData[elmNodResIdx] = inf; + } + } + } + } + + frameCountProgress.incrementProgress(); + } + + RigFemScalarResultFrames* requestedResultFrames = m_resultCollection->findOrLoadScalarResult( partIndex, resVarAddr ); + return requestedResultFrames; +} diff --git a/ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultCalculatorInitialPorosity.h b/ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultCalculatorInitialPorosity.h new file mode 100644 index 0000000000..ac875a0348 --- /dev/null +++ b/ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultCalculatorInitialPorosity.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 RigFemPartResultCalculatorInitialPorosity : public RigFemPartResultCalculator +{ +public: + explicit RigFemPartResultCalculatorInitialPorosity( RigFemPartResultsCollection& collection ); + virtual ~RigFemPartResultCalculatorInitialPorosity(); + 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 fc9e9c407c..bd42cfea4f 100644 --- a/ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.cpp +++ b/ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.cpp @@ -40,6 +40,7 @@ #include "RigFemPartResultCalculatorFOS.h" #include "RigFemPartResultCalculatorFormationIndices.h" #include "RigFemPartResultCalculatorGamma.h" +#include "RigFemPartResultCalculatorInitialPorosity.h" #include "RigFemPartResultCalculatorMudWeightWindow.h" #include "RigFemPartResultCalculatorNE.h" #include "RigFemPartResultCalculatorNodalGradients.h" @@ -184,6 +185,8 @@ RigFemPartResultsCollection::RigFemPartResultsCollection( RifGeoMechReaderInterf std::unique_ptr( new RigFemPartResultCalculatorPoreCompressibility( *this ) ) ); m_resultCalculators.push_back( std::unique_ptr( new RigFemPartResultCalculatorPorosityPermeability( *this ) ) ); + m_resultCalculators.push_back( + std::unique_ptr( new RigFemPartResultCalculatorInitialPorosity( *this ) ) ); m_resultCalculators.push_back( std::unique_ptr( new RigFemPartResultCalculatorMudWeightWindow( *this ) ) ); m_resultCalculators.push_back( @@ -661,6 +664,7 @@ std::map> fieldCompNames["COMPRESSIBILITY"].push_back( "VERTICAL" ); fieldCompNames["COMPRESSIBILITY"].push_back( "VERTICAL-RATIO" ); + fieldCompNames["PORO-PERM"].push_back( "PHI0" ); fieldCompNames["PORO-PERM"].push_back( "PHI" ); fieldCompNames["PORO-PERM"].push_back( "DPHI" ); fieldCompNames["PORO-PERM"].push_back( "PERM" ); @@ -750,6 +754,7 @@ std::map> fieldCompNames["COMPRESSIBILITY"].push_back( "VERTICAL" ); fieldCompNames["COMPRESSIBILITY"].push_back( "VERTICAL-RATIO" ); + fieldCompNames["PORO-PERM"].push_back( "PHI0" ); fieldCompNames["PORO-PERM"].push_back( "PHI" ); fieldCompNames["PORO-PERM"].push_back( "DPHI" ); fieldCompNames["PORO-PERM"].push_back( "PERM" );