Merge pull request #6344 from OPM/shear-slip-indicator-6308

#6308 GeoMech: Add shear slip indicator calculation.

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
 )

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 <http://www.gnu.org/licenses/gpl.html>
+//  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 <QString>
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+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<float>::infinity();
+
+    frameCountProgress.setNextProgressIncrement( 1u );
+
+    int frameCount = stressDataFrames->frameCount();
+    for ( int fIdx = 0; fIdx < frameCount; ++fIdx )
+    {
+        const std::vector<float>& porFrameData = porePressureDataFrames->frameData( fIdx );
+
+        const std::vector<float>& stressFrameData = stressDataFrames->frameData( fIdx );
+
+        std::vector<float>& 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;
+}

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 <http://www.gnu.org/licenses/gpl.html>
+//  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;
+};

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<RigFemPartResultCalculator>( new RigFemPartResultCalculatorTimeLapse( *this ) ) );
     m_resultCalculators.push_back(
@@ -183,6 +186,8 @@ RigFemPartResultsCollection::RigFemPartResultsCollection( RifGeoMechReaderInterf
         std::unique_ptr<RigFemPartResultCalculator>( new RigFemPartResultCalculatorPorosityPermeability( *this ) ) );
     m_resultCalculators.push_back(
         std::unique_ptr<RigFemPartResultCalculator>( new RigFemPartResultCalculatorMudWeightWindow( *this ) ) );
+    m_resultCalculators.push_back(
+        std::unique_ptr<RigFemPartResultCalculator>( new RigFemPartResultCalculatorShearSlipIndicator( *this ) ) );
     m_resultCalculators.push_back(
         std::unique_ptr<RigFemPartResultCalculator>( new RigFemPartResultCalculatorFormationIndices( *this ) ) );
 }
@@ -614,6 +619,7 @@ std::map<std::string, std::vector<std::string>>
 
             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<std::string, std::vector<std::string>>
 
             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() ) );
+    }
+}

ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.h

@@ -111,6 +111,9 @@ public:
 
     RimMudWeightWindowParameters::FractureGradientCalculationType fractureGradientCalculationTypeMudWeightWindow() const;
 
+    double waterDensityShearSlipIndicator() const;
+    void   setWaterDensityShearSlipIndicator( double waterDensity );
+
     std::map<std::string, std::vector<std::string>> scalarFieldAndComponentNames( RigFemResultPosEnum resPos );
     std::vector<std::string>                        filteredStepNames() const;
     bool                                            assertResultsLoaded( const RigFemResultAddress& resVarAddr );
@@ -224,6 +227,8 @@ private:
     std::map<RimMudWeightWindowParameters::ParameterType, QString> parameterAddresses;
     std::map<RimMudWeightWindowParameters::ParameterType, double>  parameterValues;
 
+    double m_waterDensityShearSlipIndicator;
+
     std::vector<std::unique_ptr<RigFemPartResultCalculator>> m_resultCalculators;
 
     RigStatisticsDataCache*          statistics( const RigFemResultAddress& resVarAddr );

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 );

ApplicationCode/ProjectDataModel/RimGeoMechCase.h

@@ -163,6 +163,8 @@ private:
     caf::PdmField<QString>                               m_initialPermeabilityResultAddress;
     caf::PdmField<double>                                m_permeabilityExponent;
 
+    caf::PdmField<double> m_waterDensityShearSlipIndicator;
+
     caf::PdmChildField<RimGeoMechContourMapViewCollection*> m_contourMapCollection;
     caf::PdmChildField<RimMudWeightWindowParameters*>       m_mudWeightWindowParameters;