diff --git a/ApplicationCode/GeoMech/GeoMechDataModel/CMakeLists.txt b/ApplicationCode/GeoMech/GeoMechDataModel/CMakeLists.txt index 0b1ad06d9e..5dd34d35d7 100644 --- a/ApplicationCode/GeoMech/GeoMechDataModel/CMakeLists.txt +++ b/ApplicationCode/GeoMech/GeoMechDataModel/CMakeLists.txt @@ -98,6 +98,8 @@ add_library( ${PROJECT_NAME} RigFemPartResultCalculatorPorosityPermeability.cpp RigFemPartResultCalculatorMudWeightWindow.h RigFemPartResultCalculatorMudWeightWindow.cpp + RigFemPartResultCalculatorShearSlipIndicator.h + RigFemPartResultCalculatorShearSlipIndicator.cpp RimGeoMechGeometrySelectionItem.h RimGeoMechGeometrySelectionItem.cpp ) diff --git a/ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultCalculatorShearSlipIndicator.cpp b/ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultCalculatorShearSlipIndicator.cpp new file mode 100644 index 0000000000..0a1917e231 --- /dev/null +++ b/ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultCalculatorShearSlipIndicator.cpp @@ -0,0 +1,163 @@ +///////////////////////////////////////////////////////////////////////////////// +// +// 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 "RigFemPartResultCalculatorShearSlipIndicator.h" + +#include "RigFemPart.h" +#include "RigFemPartCollection.h" +#include "RigFemPartGrid.h" +#include "RigFemPartResultsCollection.h" +#include "RigFemResultAddress.h" +#include "RigFemScalarResultFrames.h" +#include "RigGeoMechWellLogExtractor.h" + +#include "cafProgressInfo.h" + +#include + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +RigFemPartResultCalculatorShearSlipIndicator::RigFemPartResultCalculatorShearSlipIndicator( RigFemPartResultsCollection& collection ) + : RigFemPartResultCalculator( collection ) +{ +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +RigFemPartResultCalculatorShearSlipIndicator::~RigFemPartResultCalculatorShearSlipIndicator() +{ +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +bool RigFemPartResultCalculatorShearSlipIndicator::isMatching( const RigFemResultAddress& resVarAddr ) const +{ + return ( resVarAddr.fieldName == "ST" || resVarAddr.componentName == "DPN" ); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +RigFemScalarResultFrames* + RigFemPartResultCalculatorShearSlipIndicator::calculate( int partIndex, const RigFemResultAddress& resVarAddr ) +{ + CVF_ASSERT( isMatching( resVarAddr ) ); + + caf::ProgressInfo frameCountProgress( m_resultCollection->frameCount() * 3, "" ); + frameCountProgress.setProgressDescription( "Calculating Shear Slip Indicator." ); + + // Pore pressure + frameCountProgress.setNextProgressIncrement( m_resultCollection->frameCount() ); + RigFemScalarResultFrames* porePressureDataFrames = + m_resultCollection->findOrLoadScalarResult( partIndex, + RigFemResultAddress( resVarAddr.resultPosType, "POR-Bar", "" ) ); + frameCountProgress.incrementProgress(); + + // Total vertical stress (ST.S33) + frameCountProgress.setNextProgressIncrement( m_resultCollection->frameCount() ); + RigFemScalarResultFrames* stressDataFrames = + m_resultCollection->findOrLoadScalarResult( partIndex, + RigFemResultAddress( resVarAddr.resultPosType, "ST", "S33" ) ); + frameCountProgress.incrementProgress(); + + frameCountProgress.setNextProgressIncrement( m_resultCollection->frameCount() ); + + RigFemScalarResultFrames* shearSlipIndicatorFrames = + m_resultCollection->createScalarResult( partIndex, RigFemResultAddress( resVarAddr.resultPosType, "ST", "DPN" ) ); + + const RigFemPart* femPart = m_resultCollection->parts()->part( partIndex ); + const RigFemPartGrid* femPartGrid = femPart->getOrCreateStructGrid(); + + float inf = std::numeric_limits::infinity(); + + frameCountProgress.setNextProgressIncrement( 1u ); + + int frameCount = stressDataFrames->frameCount(); + for ( int fIdx = 0; fIdx < frameCount; ++fIdx ) + { + const std::vector& porFrameData = porePressureDataFrames->frameData( fIdx ); + + const std::vector& stressFrameData = stressDataFrames->frameData( fIdx ); + + std::vector& shearSlipIndicatorFrameData = shearSlipIndicatorFrames->frameData( fIdx ); + + size_t valCount = stressFrameData.size(); + shearSlipIndicatorFrameData.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 ( femPart->isHexahedron( elmIdx ) ) + { + for ( int elmNodIdx = 0; elmNodIdx < elmNodeCount; ++elmNodIdx ) + { + // Use hydrostatic pressure from cell centroid. + // Use centroid to avoid intra-element differences + cvf::Vec3d cellCentroid = femPartGrid->cellCentroid( elmIdx ); + double cellCentroidTvdMSL = -cellCentroid.z(); + + double waterDensity = m_resultCollection->waterDensityShearSlipIndicator(); + double cellCenterHydroStaticPressure = + RigGeoMechWellLogExtractor::hydroStaticPorePressureAtDepth( cellCentroidTvdMSL, waterDensity ); + + size_t elmNodResIdx = femPart->elementNodeResultIdx( elmIdx, elmNodIdx ); + if ( elmNodResIdx < stressFrameData.size() ) + { + // Pore pressure (unit: Bar) + float porePressureBar = porFrameData[elmNodResIdx]; + float totalVerticalStress = stressFrameData[elmNodResIdx]; + + float shearSlipIndicator = inf; + if ( porePressureBar != inf && totalVerticalStress - cellCenterHydroStaticPressure != 0.0 ) + { + shearSlipIndicator = ( porePressureBar - cellCenterHydroStaticPressure ) / + ( totalVerticalStress - cellCenterHydroStaticPressure ); + } + + shearSlipIndicatorFrameData[elmNodResIdx] = shearSlipIndicator; + } + } + } + else + { + for ( int elmNodIdx = 0; elmNodIdx < elmNodeCount; ++elmNodIdx ) + { + size_t elmNodResIdx = femPart->elementNodeResultIdx( elmIdx, elmNodIdx ); + if ( elmNodResIdx < stressFrameData.size() ) + { + shearSlipIndicatorFrameData[elmNodResIdx] = inf; + } + } + } + } + + frameCountProgress.incrementProgress(); + } + + RigFemScalarResultFrames* requestedResultFrames = m_resultCollection->findOrLoadScalarResult( partIndex, resVarAddr ); + return requestedResultFrames; +} diff --git a/ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultCalculatorShearSlipIndicator.h b/ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultCalculatorShearSlipIndicator.h new file mode 100644 index 0000000000..8da3fb16f2 --- /dev/null +++ b/ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultCalculatorShearSlipIndicator.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 RigFemPartResultCalculatorShearSlipIndicator : public RigFemPartResultCalculator +{ +public: + explicit RigFemPartResultCalculatorShearSlipIndicator( RigFemPartResultsCollection& collection ); + virtual ~RigFemPartResultCalculatorShearSlipIndicator(); + 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 9f582fba10..fc9e9c407c 100644 --- a/ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.cpp +++ b/ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.cpp @@ -55,6 +55,7 @@ #include "RigFemPartResultCalculatorSM.h" #include "RigFemPartResultCalculatorShearSE.h" #include "RigFemPartResultCalculatorShearST.h" +#include "RigFemPartResultCalculatorShearSlipIndicator.h" #include "RigFemPartResultCalculatorStressAnisotropy.h" #include "RigFemPartResultCalculatorStressGradients.h" #include "RigFemPartResultCalculatorSurfaceAlignedStress.h" @@ -126,6 +127,8 @@ RigFemPartResultsCollection::RigFemPartResultsCollection( RifGeoMechReaderInterf m_nonReservoirPorePressureAddressMudWeightWindow = ""; m_hydrostaticMultiplierPPNonResMudWeightWindow = 1.0; + m_waterDensityShearSlipIndicator = 1.03; + m_resultCalculators.push_back( std::unique_ptr( new RigFemPartResultCalculatorTimeLapse( *this ) ) ); m_resultCalculators.push_back( @@ -183,6 +186,8 @@ RigFemPartResultsCollection::RigFemPartResultsCollection( RifGeoMechReaderInterf std::unique_ptr( new RigFemPartResultCalculatorPorosityPermeability( *this ) ) ); m_resultCalculators.push_back( std::unique_ptr( new RigFemPartResultCalculatorMudWeightWindow( *this ) ) ); + m_resultCalculators.push_back( + std::unique_ptr( new RigFemPartResultCalculatorShearSlipIndicator( *this ) ) ); m_resultCalculators.push_back( std::unique_ptr( new RigFemPartResultCalculatorFormationIndices( *this ) ) ); } @@ -614,6 +619,7 @@ std::map> fieldCompNames["ST"].push_back( "SM" ); fieldCompNames["ST"].push_back( "Q" ); + fieldCompNames["ST"].push_back( "DPN" ); for ( auto& s : stressComponentNames ) { @@ -697,6 +703,7 @@ std::map> fieldCompNames["ST"].push_back( "SM" ); fieldCompNames["ST"].push_back( "Q" ); + fieldCompNames["ST"].push_back( "DPN" ); fieldCompNames["ST"].push_back( "S11" ); fieldCompNames["ST"].push_back( "S22" ); @@ -1736,3 +1743,24 @@ RimMudWeightWindowParameters::FractureGradientCalculationType { return m_fractureGradientCalculationTypeMudWeightWindow; } + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +double RigFemPartResultsCollection::waterDensityShearSlipIndicator() const +{ + return m_waterDensityShearSlipIndicator; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +void RigFemPartResultsCollection::setWaterDensityShearSlipIndicator( double waterDensity ) +{ + m_waterDensityShearSlipIndicator = waterDensity; + + for ( auto elementType : {RIG_ELEMENT_NODAL, RIG_INTEGRATION_POINT} ) + { + deleteResult( RigFemResultAddress( elementType, "ST", "DPN", RigFemResultAddress::allTimeLapsesValue() ) ); + } +} diff --git a/ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.h b/ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.h index 79206bb121..bc17c5bb69 100644 --- a/ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.h +++ b/ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.h @@ -111,6 +111,9 @@ public: RimMudWeightWindowParameters::FractureGradientCalculationType fractureGradientCalculationTypeMudWeightWindow() const; + double waterDensityShearSlipIndicator() const; + void setWaterDensityShearSlipIndicator( double waterDensity ); + std::map> scalarFieldAndComponentNames( RigFemResultPosEnum resPos ); std::vector filteredStepNames() const; bool assertResultsLoaded( const RigFemResultAddress& resVarAddr ); @@ -224,6 +227,8 @@ private: std::map parameterAddresses; std::map parameterValues; + double m_waterDensityShearSlipIndicator; + std::vector> m_resultCalculators; RigStatisticsDataCache* statistics( const RigFemResultAddress& resVarAddr ); diff --git a/ApplicationCode/ProjectDataModel/RimGeoMechCase.cpp b/ApplicationCode/ProjectDataModel/RimGeoMechCase.cpp index f1a423a20f..8fe9b309d8 100644 --- a/ApplicationCode/ProjectDataModel/RimGeoMechCase.cpp +++ b/ApplicationCode/ProjectDataModel/RimGeoMechCase.cpp @@ -177,6 +177,9 @@ RimGeoMechCase::RimGeoMechCase( void ) CAF_PDM_InitField( &m_permeabilityExponent, "PermeabilityExponent", 1.0, "Permeability Exponent", "", "", "" ); m_permeabilityExponent.uiCapability()->setUiEditorTypeName( caf::PdmUiDoubleValueEditor::uiEditorTypeName() ); + CAF_PDM_InitField( &m_waterDensityShearSlipIndicator, "WaterDensityShearSlipIndicator", 1.03, "Water Density", "", "", "" ); + m_waterDensityShearSlipIndicator.uiCapability()->setUiEditorTypeName( caf::PdmUiDoubleValueEditor::uiEditorTypeName() ); + CAF_PDM_InitFieldNoDefault( &m_contourMapCollection, "ContourMaps", "2d Contour Maps", "", "", "" ); m_contourMapCollection = new RimGeoMechContourMapViewCollection; m_contourMapCollection.uiCapability()->setUiTreeHidden( true ); @@ -378,6 +381,7 @@ RimGeoMechCase::CaseOpenStatus RimGeoMechCase::openGeoMechCase( std::string* err } geoMechCaseData->femPartResults()->addElementPropertyFiles( fileNames ); geoMechCaseData->femPartResults()->setCalculationParameters( m_cohesion, cvf::Math::toRadians( m_frictionAngleDeg() ) ); + geoMechCaseData->femPartResults()->setWaterDensityShearSlipIndicator( m_waterDensityShearSlipIndicator ); m_geoMechCaseData = geoMechCaseData; @@ -870,6 +874,11 @@ void RimGeoMechCase::fieldChangedByUi( const caf::PdmFieldHandle* changedField, updateConnectedViews(); } + else if ( changedField == &m_waterDensityShearSlipIndicator ) + { + rigCaseData->femPartResults()->setWaterDensityShearSlipIndicator( m_waterDensityShearSlipIndicator ); + updateConnectedViews(); + } else if ( changedField == &m_reloadElementPropertyFileCommand ) { m_reloadElementPropertyFileCommand = false; @@ -1134,6 +1143,9 @@ void RimGeoMechCase::defineUiOrdering( QString uiConfigName, caf::PdmUiOrdering& m_initialPermeabilityType != RimGeoMechCase::InitialPermeabilityType::INITIAL_PERMEABILITY_PER_ELEMENT ); permeabilityGroup->add( &m_permeabilityExponent ); + caf::PdmUiGroup* shearSlipIndicatorGroup = uiOrdering.addNewGroup( "Shear Slip Indicator" ); + shearSlipIndicatorGroup->add( &m_waterDensityShearSlipIndicator ); + caf::PdmUiGroup* timeStepFilterGroup = uiOrdering.addNewGroup( "Time Step Filter" ); timeStepFilterGroup->setCollapsedByDefault( true ); m_timeStepFilter->uiOrdering( uiConfigName, *timeStepFilterGroup ); diff --git a/ApplicationCode/ProjectDataModel/RimGeoMechCase.h b/ApplicationCode/ProjectDataModel/RimGeoMechCase.h index 1c528aff17..1cc99737b0 100644 --- a/ApplicationCode/ProjectDataModel/RimGeoMechCase.h +++ b/ApplicationCode/ProjectDataModel/RimGeoMechCase.h @@ -163,6 +163,8 @@ private: caf::PdmField m_initialPermeabilityResultAddress; caf::PdmField m_permeabilityExponent; + caf::PdmField m_waterDensityShearSlipIndicator; + caf::PdmChildField m_contourMapCollection; caf::PdmChildField m_mudWeightWindowParameters;