Well target mapping: option to avoid growing across faults

Add AvoidFaults option that blocks cluster growth across fault faces and
fault NNCs, and generate a WELL_TARGET_FAULT_FACES result with the number
of fault faces per cell.

ApplicationLibCode/ProjectDataModel/RimWellTargetMapping.cpp

@@ -172,6 +172,14 @@ RimWellTargetMapping::RimWellTargetMapping()
                                 "",
                                 "Limits where clusters can grow. If empty, clusters can grow anywhere in the model." );
 
+    CAF_PDM_InitField( &m_avoidFaults,
+                       "AvoidFaults",
+                       false,
+                       "Avoid Growing Across Faults",
+                       "",
+                       "Do not grow clusters across fault faces or fault NNCs. Also generates a result with the number of fault faces per "
+                       "cell." );
+
     CAF_PDM_InitFieldNoDefault( &m_ensembleStatisticsCase, "EnsembleStatisticsCase", "Ensemble Statistics Case" );
 
     m_minimumSaturationOil = cvf::UNDEFINED_DOUBLE;
@@ -510,6 +518,7 @@ void RimWellTargetMapping::defineUiOrdering( QString uiConfigName, caf::PdmUiOrd
     auto hasEnsembleParent = firstAncestorOrThisOfType<RimEclipseCaseEnsemble>() != nullptr || gridEnsemble != nullptr;
     uiOrdering.add( &m_startDataFilter );
     uiOrdering.add( &m_dataFilter );
+    uiOrdering.add( &m_avoidFaults );
 
     caf::PdmUiGroup* minimumCellValuesGroup = uiOrdering.addNewGroup( "Minimum Cell Values" );
     if ( showOilOptions ) minimumCellValuesGroup->add( &m_saturationOil );
@@ -570,6 +579,7 @@ RigWellTargetMapping::ClusteringLimits RimWellTargetMapping::getClusteringLimits
              .transmissibility = m_transmissibility,
              .maxNumTargets    = m_maxNumTargets,
              .maxIterations    = m_maxIterations,
+             .avoidFaults      = m_avoidFaults,
              .filterAddress    = m_resultDefinition->eclipseResultAddress(),
              .startDataFilter  = m_startDataFilter(),
              .growDataFilter   = m_dataFilter() };

ApplicationLibCode/ProjectDataModel/RimWellTargetMapping.h

@@ -110,6 +110,8 @@ private:
     caf::PdmPtrField<RimCellFilter*> m_startDataFilter;
     caf::PdmPtrField<RimCellFilter*> m_dataFilter; // grow filter
 
+    caf::PdmField<bool> m_avoidFaults;
+
     caf::PdmChildField<RimEclipseCase*> m_ensembleStatisticsCase;
 
     double m_minimumSaturationOil;

ApplicationLibCode/ReservoirDataModel/Well/RigWellTargetMapping.cpp

@@ -146,6 +146,8 @@ void RigWellTargetMapping::generateCandidates( RimEclipseCase*            eclips
         data.transmissibilityNNC = nullptr;
     }
 
+    if ( limits.avoidFaults && mainGrid != nullptr ) data.isFaultNnc = RigWellTargetMappingTools::findFaultNncs( mainGrid );
+
     // Evaluate for each realization after opening its grid, using the mapping's time step.
     auto computeFilterVector = [&]( RimCellFilter* dataFilter ) -> std::vector<double>
     {
@@ -245,6 +247,15 @@ void RigWellTargetMapping::generateCandidates( RimEclipseCase*            eclips
     QString resultName = RigWellTargetMapping::wellTargetResultName();
     RigWellTargetMappingTools::createResultVector( *eclipseCase, resultName, clusters, timeStepIdx );
 
+    if ( limits.avoidFaults && mainGrid != nullptr )
+    {
+        RigWellTargetMappingTools::createResultVector( *eclipseCase,
+                                                       wellTargetFaultFacesResultName(),
+                                                       RigWellTargetMappingTools::countFaultFacesPerActiveCell( mainGrid,
+                                                                                                                resultsData->activeCellInfo() ),
+                                                       timeStepIdx );
+    }
+
     std::vector<RigWellTargetMappingTools::ClusterStatistics> statistics =
         RigWellTargetMappingTools::generateStatistics( eclipseCase, data.pressure, data.permeabilityX, numClustersFound, timeStepIdx, resultName );
     std::vector<double> totalPorvSoil( clusters.size(), std::numeric_limits<double>::infinity() );
@@ -520,3 +531,11 @@ QString RigWellTargetMapping::wellTargetResultName()
 {
     return "WELL_TARGET_NUM";
 }
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+QString RigWellTargetMapping::wellTargetFaultFacesResultName()
+{
+    return "WELL_TARGET_FAULT_FACES";
+}

ApplicationLibCode/ReservoirDataModel/Well/RigWellTargetMapping.h

@@ -67,6 +67,7 @@ public:
         double                  transmissibility;
         int                     maxNumTargets;
         int                     maxIterations;
+        bool                    avoidFaults = false; // do not grow across fault faces and fault NNCs
         RigEclipseResultAddress filterAddress;
         RimCellFilter*          startDataFilter = nullptr; // falls back to growDataFilter if null
         RimCellFilter*          growDataFilter  = nullptr;
@@ -102,4 +103,5 @@ public:
                                                            double                              expandBoundingBoxZPercent );
 
     static QString wellTargetResultName();
+    static QString wellTargetFaultFacesResultName();
 };

ApplicationLibCode/ReservoirDataModel/Well/RigWellTargetMappingTools.cpp

@@ -24,10 +24,12 @@
 #include "RiaResultNames.h"
 #include "RiaWeightedMeanCalculator.h"
 
+#include "RigActiveCellInfo.h"
 #include "RigCaseCellResultsData.h"
 #include "RigCell.h"
 #include "RigEclipseCaseData.h"
 #include "RigEclipseResultInfo.h"
+#include "RigFault.h"
 #include "RigFloodingSettings.h"
 #include "RigGridBase.h"
 #include "RigHydrocarbonFlowTools.h"
@@ -159,6 +161,46 @@ void RigWellTargetMappingTools::assignClusterIdToCells( const RigActiveCellInfo&
     }
 }
 
+//--------------------------------------------------------------------------------------------------
+/// Flags the NNCs that connect cells across a fault
+//--------------------------------------------------------------------------------------------------
+std::vector<bool> RigWellTargetMappingTools::findFaultNncs( RigMainGrid* mainGrid )
+{
+    if ( mainGrid->nncData() == nullptr ) return {};
+
+    std::vector<bool> isFaultNnc( mainGrid->nncData()->allConnections().size(), false );
+    for ( const auto& fault : mainGrid->faults() )
+    {
+        for ( size_t nncIdx : fault->connectionIndices() )
+        {
+            if ( nncIdx < isFaultNnc.size() ) isFaultNnc[nncIdx] = true;
+        }
+    }
+    return isFaultNnc;
+}
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+std::vector<int> RigWellTargetMappingTools::countFaultFacesPerActiveCell( const RigMainGrid* mainGrid, const RigActiveCellInfo* activeCellInfo )
+{
+    const std::vector<CellFaceType> faces =
+        { CellFaceType::POS_I, CellFaceType::NEG_I, CellFaceType::POS_J, CellFaceType::NEG_J, CellFaceType::POS_K, CellFaceType::NEG_K };
+
+    std::vector<int> result( activeCellInfo->reservoirActiveCellCount(), 0 );
+    for ( size_t reservoirIdx = 0; reservoirIdx < activeCellInfo->reservoirCellCount(); reservoirIdx++ )
+    {
+        ActiveCellIndex resultIdx = activeCellInfo->cellResultIndex( ReservoirCellIndex( reservoirIdx ) );
+        if ( resultIdx.value() == cvf::UNDEFINED_SIZE_T ) continue;
+
+        for ( CellFaceType face : faces )
+        {
+            if ( mainGrid->findFaultFromCellIndexAndCellFace( reservoirIdx, face ) != nullptr ) result[resultIdx.value()]++;
+        }
+    }
+    return result;
+}
+
 //--------------------------------------------------------------------------------------------------
 ///
 //--------------------------------------------------------------------------------------------------
@@ -442,6 +484,8 @@ std::vector<ReservoirCellIndex> RigWellTargetMappingTools::findCandidates( RimEc
                 ActiveCellIndex    neighborResultIndex = resultsData->activeCellInfo()->cellResultIndex( neighborResvCellIdx );
                 if ( neighborResultIndex.value() != cvf::UNDEFINED_SIZE_T && clusters[neighborResultIndex.value()] == 0 )
                 {
+                    if ( limits.avoidFaults && mainGrid->findFaultFromCellIndexAndCellFace( cellIdx.value(), face ) != nullptr ) continue;
+
                     double permeability     = data.permeabilityX[neighborResultIndex.value()];
                     double transmissibility = getTransmissibilityValueForFace( data.transmissibilityX,
                                                                                data.transmissibilityY,
@@ -471,7 +515,8 @@ std::vector<ReservoirCellIndex> RigWellTargetMappingTools::findCandidates( RimEc
             for ( auto& [cellInfo, nncResultIdx] : nncCells )
             {
                 auto& [otherCellIdx, face] = cellInfo;
-                double transmissibility    = data.transmissibilityNNC->at( nncResultIdx );
+                if ( limits.avoidFaults && nncResultIdx < data.isFaultNnc.size() && data.isFaultNnc[nncResultIdx] ) continue;
+                double transmissibility = data.transmissibilityNNC->at( nncResultIdx );
 
                 ActiveCellIndex otherResultIndex = resultsData->activeCellInfo()->cellResultIndex( otherCellIdx );
                 if ( otherResultIndex.value() == cvf::UNDEFINED_SIZE_T ) continue;

ApplicationLibCode/ReservoirDataModel/Well/RigWellTargetMappingTools.h

@@ -63,6 +63,7 @@ public:
         std::vector<double>        transmissibilityY;
         std::vector<double>        transmissibilityZ;
         const std::vector<double>* transmissibilityNNC;
+        std::vector<bool>          isFaultNnc; // empty if fault barriers are not used
     };
 
     class ClusterStatistics
@@ -100,6 +101,9 @@ public:
         double pressure;
     };
 
+    static std::vector<bool> findFaultNncs( RigMainGrid* mainGrid );
+    static std::vector<int>  countFaultFacesPerActiveCell( const RigMainGrid* mainGrid, const RigActiveCellInfo* activeCellInfo );
+
     static double getValueForFace( const std::vector<double>& x,
                                    const std::vector<double>& y,
                                    const std::vector<double>& z,

ApplicationLibCode/UnitTests/RimWellTargetMapping-Test.cpp

@@ -126,6 +126,33 @@ TEST( RimWellTargetMappingTest, StartFilterOptionsUseOwningCase )
     EXPECT_EQ( caf::PdmOptionItemInfo( "Case filter", filter ).value(), options[1].value() );
 }
 
+TEST( RimWellTargetMappingTest, AvoidFaultsGeneratesFaultFaceResult )
+{
+    RigEclipseResultAddress address( RiaDefines::ResultCatType::GENERATED, RigWellTargetMapping::wellTargetFaultFacesResultName() );
+
+    for ( bool avoidFaults : { false, true } )
+    {
+        auto eclipseCase = openCase( "Real0", "BRUGGE_0000.EGRID" );
+        ASSERT_NE( nullptr, eclipseCase );
+
+        auto mapping = new TestableWellTargetMapping;
+        eclipseCase->addWellTargetMapping( mapping );
+        setField( *mapping, "TimeStep", 1 );
+        setField( *mapping, "Iterations", 10 );
+        setField( *mapping, "MaxNumTargets", 1 );
+        setField( *mapping, "AvoidFaults", avoidFaults );
+        mapping->initAfterRead();
+
+        auto results = eclipseCase->results( RiaDefines::PorosityModelType::MATRIX_MODEL );
+        ASSERT_EQ( avoidFaults, results->hasResultEntry( address ) );
+        if ( !avoidFaults ) continue;
+
+        const auto& faultFaces = results->cellScalarResults( address, 1 );
+        ASSERT_EQ( results->activeCellInfo()->reservoirActiveCellCount(), faultFaces.size() );
+        EXPECT_GT( std::count_if( faultFaces.begin(), faultFaces.end(), []( double v ) { return v > 0.0; } ), 0 );
+    }
+}
+
 TEST( RimWellTargetMappingTest, FilterOptionsUseOwningEnsemble )
 {
     RimReservoirGridEnsemble ensemble;