Use NNC info for well target generation and streamlines (#12279)

* Add support for NNC when growing clusters
* Trace streamlines through NNCs, too

ApplicationLibCode/ProjectDataModel/Streamlines/RimStreamlineDataAccess.cpp

@@ -23,6 +23,7 @@
 #include "RigEclipseCaseData.h"
 #include "RigEclipseResultAddress.h"
 #include "RigMainGrid.h"
+#include "RigNNCData.h"
 #include "RigResultAccessor.h"
 #include "RigResultAccessorFactory.h"
 
@@ -75,6 +76,24 @@ bool RimStreamlineDataAccess::setupDataAccess( RigMainGrid* grid, RigEclipseCase
         }
     }
 
+    auto nncData = grid->nncData();
+    m_nncData.clear();
+    if ( nncData != nullptr )
+    {
+        m_nncData[RiaDefines::PhaseType::WATER_PHASE] =
+            nncData->dynamicConnectionScalarResultByName( RiaDefines::propertyNameFluxWat(), timeIdx );
+        m_nncData[RiaDefines::PhaseType::OIL_PHASE] =
+            nncData->dynamicConnectionScalarResultByName( RiaDefines::propertyNameFluxOil(), timeIdx );
+        m_nncData[RiaDefines::PhaseType::GAS_PHASE] =
+            nncData->dynamicConnectionScalarResultByName( RiaDefines::propertyNameFluxGas(), timeIdx );
+    }
+    else
+    {
+        m_nncData[RiaDefines::PhaseType::WATER_PHASE] = nullptr;
+        m_nncData[RiaDefines::PhaseType::OIL_PHASE]   = nullptr;
+        m_nncData[RiaDefines::PhaseType::GAS_PHASE]   = nullptr;
+    }
+
     return true;
 }
 
@@ -229,3 +248,50 @@ double RimStreamlineDataAccess::combinedFaceRate( RigCell
 
     return retValue;
 }
+
+//--------------------------------------------------------------------------------------------------
+/// Return the nnc scalar value for the given connection, by combining flow for all specified phases
+/// Positive values is flow from cell1 to cell2, negative values is the other way
+//--------------------------------------------------------------------------------------------------
+double RimStreamlineDataAccess::combinedNNCRate( size_t                           resultIdx,
+                                                 std::list<RiaDefines::PhaseType> phases,
+                                                 double                           direction,
+                                                 RiaDefines::PhaseType&           outDominantPhase ) const
+{
+    double retValue  = 0.0;
+    outDominantPhase = phases.front();
+
+    double max = 0.0;
+
+    for ( auto phase : phases )
+    {
+        double tmp = 0.0;
+
+        if ( m_nncData.at( phase ) != nullptr )
+        {
+            if ( resultIdx < m_nncData.at( phase )->size() )
+            {
+                tmp = m_nncData.at( phase )->at( resultIdx );
+            }
+        }
+
+        if ( tmp * direction > max )
+        {
+            outDominantPhase = phase;
+            max              = std::abs( tmp );
+        }
+
+        retValue += tmp;
+    }
+
+    return retValue;
+}
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+const RigConnection& RimStreamlineDataAccess::nncConnection( size_t idx ) const
+{
+    auto nncData = m_grid->nncData();
+    return nncData->availableConnections()[idx];
+}

ApplicationLibCode/ProjectDataModel/Streamlines/RimStreamlineDataAccess.h

@@ -33,6 +33,7 @@ class RigMainGrid;
 class RigGridBase;
 class RigResultAccessor;
 class RigEclipseCaseData;
+class RigConnection;
 
 //--------------------------------------------------------------------------------------------------
 /// Specialized data access for streamline generation. Operates using flow rate in meters/day
@@ -54,7 +55,14 @@ public:
                              double                             direction,
                              RiaDefines::PhaseType&             dominantPhaseOut ) const;
 
-    const RigMainGrid* grid() const { return m_grid; }
+    double combinedNNCRate( size_t                           resultIdx,
+                            std::list<RiaDefines::PhaseType> phases,
+                            double                           direction,
+                            RiaDefines::PhaseType&           outDominantPhase ) const;
+
+    RigMainGrid* grid() const { return m_grid; }
+
+    const RigConnection& nncConnection( size_t idx ) const;
 
 protected:
     cvf::ref<RigResultAccessor> getDataAccessor( cvf::StructGridInterface::FaceType faceIdx, RiaDefines::PhaseType phase, int timeIdx );
@@ -68,4 +76,5 @@ private:
     RigEclipseCaseData* m_data;
 
     std::map<RiaDefines::PhaseType, std::vector<cvf::ref<RigResultAccessor>>> m_dataAccess;
+    std::map<RiaDefines::PhaseType, const std::vector<double>*>               m_nncData;
 };

ApplicationLibCode/ProjectDataModel/Streamlines/RimStreamlineGenerator.cpp

@@ -20,6 +20,8 @@
 
 #include "RigCell.h"
 #include "RigMainGrid.h"
+#include "RigNNCData.h"
+#include "RigNncConnection.h"
 
 #include "RimStreamline.h"
 #include "RimStreamlineDataAccess.h"
@@ -32,18 +34,22 @@
 class StreamlineSeedPoint
 {
 public:
-    StreamlineSeedPoint( double rate, size_t cellIdx, cvf::StructGridInterface::FaceType faceIdx )
+    StreamlineSeedPoint( double rate, size_t index, RimStreamlineGenerator::CellFaceType faceIdx )
         : m_rate( rate )
-        , m_cellIdx( cellIdx )
+        , m_index( index )
         , m_faceIdx( faceIdx ){};
     ~StreamlineSeedPoint(){};
 
-    bool operator<( const StreamlineSeedPoint& other ) const { return m_rate < other.m_rate; };
+    bool operator<( const StreamlineSeedPoint& other ) const { return rate() < other.rate(); };
 
-public:
+    double                               rate() const { return m_rate; };
-    double                             m_rate;
+    size_t                               index() const { return m_index; };
-    size_t                             m_cellIdx;
+    RimStreamlineGenerator::CellFaceType face() const { return m_faceIdx; };
-    cvf::StructGridInterface::FaceType m_faceIdx;
+
+private:
+    double                               m_rate;
+    size_t                               m_index;
+    RimStreamlineGenerator::CellFaceType m_faceIdx;
 };
 
 //--------------------------------------------------------------------------------------------------
@@ -77,19 +83,27 @@ void RimStreamlineGenerator::generateTracer( RigCell cell, double direction, QSt
         double rate = m_dataAccess->combinedFaceRate( cell, faceIdx, m_phases, direction, dominantPhase ) * direction;
         if ( rate > m_flowThreshold )
         {
-            m_seeds.push( StreamlineSeedPoint( rate, cellIdx, faceIdx ) );
+            m_seeds.emplace( rate, cellIdx, faceIdx );
         }
     }
 
+    // add any nnc connections from this cell, too
+    auto nncSeeds = nncCandidates( cell.mainGridCellIndex(), m_phases, direction );
+    for ( auto& seed : nncSeeds )
+    {
+        m_seeds.push( seed );
+    }
+
+    // keep adding streamlines until we have no more seeds left in our prioritized (by flow rate) list
     while ( !m_seeds.empty() )
     {
-        const size_t                             cellIdx = m_seeds.top().m_cellIdx;
+        auto&                                    seed    = m_seeds.top();
-        const cvf::StructGridInterface::FaceType faceIdx = m_seeds.top().m_faceIdx;
+        const size_t                             index   = seed.index();
+        const cvf::StructGridInterface::FaceType faceIdx = seed.face();
         m_seeds.pop();
 
         RimStreamline* streamline = new RimStreamline( simWellName );
-
+        growStreamline( streamline, index, faceIdx, direction );
-        growStreamline( streamline, cellIdx, faceIdx, direction );
 
         if ( direction < 0.0 ) streamline->reverse();
 
@@ -105,42 +119,86 @@ void RimStreamlineGenerator::generateTracer( RigCell cell, double direction, QSt
 //--------------------------------------------------------------------------------------------------
 ///
 //--------------------------------------------------------------------------------------------------
-void RimStreamlineGenerator::growStreamline( RimStreamline* streamline, size_t cellIdx, cvf::StructGridInterface::FaceType faceIdx, double direction )
+void RimStreamlineGenerator::growStreamline( RimStreamline* streamline, size_t index, CellFaceType faceIdx, double direction )
 {
-    // get the cell
+    double                rate = 0.0;
-    RigCell cell = m_dataAccess->grid()->cell( cellIdx );
-
-    // get rate
     RiaDefines::PhaseType dominantPhaseOut;
+    size_t                currentCellIdx = 0;
+    RigCell               currentCell;
 
-    double rate = m_dataAccess->combinedFaceRate( cell, faceIdx, m_phases, direction, dominantPhaseOut );
+    while ( true )
-
-    // if we go backwards from a producer, the rate needs to be flipped
-    rate *= direction;
-
-    // grow from start cell center to face center, exiting if we reach the max length
-    if ( !growStreamlineFromTo( streamline, cell.center(), cell.faceCenter( faceIdx ), rate, dominantPhaseOut ) ) return;
-
-    while ( rate >= m_flowThreshold )
     {
-        // find next cell and entry face
+        // no face? -> nnc connection
-        cell = cell.neighborCell( faceIdx );
+        if ( faceIdx == CellFaceType::NO_FACE )
-        if ( cell.isInvalid() ) break;
+        {
-        faceIdx = cvf::StructGridInterface::oppositeFace( faceIdx );
+            const auto    nncIdx = index;
+            RigConnection nnc    = m_dataAccess->nncConnection( nncIdx );
 
-        // grow from given face center to cell center, exiting if we reach the max length
+            size_t localCellIdx1 = 0;
-        if ( !growStreamlineFromTo( streamline, cell.faceCenter( faceIdx ), cell.center(), rate, dominantPhaseOut ) ) break;
+            auto   localGrid1    = m_dataAccess->grid()->gridAndGridLocalIdxFromGlobalCellIdx( nnc.c1GlobIdx(), &localCellIdx1 );
+            size_t localCellIdx2 = 0;
+            auto   localGrid2    = m_dataAccess->grid()->gridAndGridLocalIdxFromGlobalCellIdx( nnc.c2GlobIdx(), &localCellIdx2 );
 
-        const size_t cellIdx = cell.gridLocalCellIndex();
+            // no support for LGRs, yet
+            if ( localGrid1->gridId() != 0 ) return;
+            if ( localGrid2->gridId() != 0 ) return;
 
-        if ( m_visitedCells.count( cellIdx ) > 0 ) break;
+            RigCell cell1 = m_dataAccess->grid()->cell( localCellIdx1 );
-        if ( m_wellCells.count( cellIdx ) > 0 ) break;
+            RigCell cell2 = m_dataAccess->grid()->cell( localCellIdx2 );
 
-        m_visitedCells.insert( cellIdx );
+            // get rate
+            rate = m_dataAccess->combinedNNCRate( nncIdx, m_phases, direction, dominantPhaseOut );
+            // if we go backwards from a producer, the rate needs to be flipped
+            rate *= direction;
+
+            // grow from cell1 center to cell2 center, exiting if we reach max length
+            if ( !growStreamlineFromTo( streamline, cell1.center(), cell2.center(), rate, dominantPhaseOut ) ) return;
+
+            // give up if too low rate
+            if ( rate < m_flowThreshold ) return;
+
+            currentCellIdx = localCellIdx2;
+            currentCell    = cell2;
+        }
+        else
+        {
+            // get the current cell
+            const auto cellIdx = index;
+            RigCell    cell    = m_dataAccess->grid()->cell( cellIdx );
+
+            // get rate
+            rate = m_dataAccess->combinedFaceRate( cell, faceIdx, m_phases, direction, dominantPhaseOut );
+
+            // if we go backwards from a producer, the rate needs to be flipped
+            rate *= direction;
+
+            // grow from start cell center to face center, exiting if we reach the max length
+            if ( !growStreamlineFromTo( streamline, cell.center(), cell.faceCenter( faceIdx ), rate, dominantPhaseOut ) ) return;
+
+            // find next cell and entry face
+            currentCell = cell.neighborCell( faceIdx );
+            if ( currentCell.isInvalid() ) return;
+            faceIdx = cvf::StructGridInterface::oppositeFace( faceIdx );
+
+            // grow from given face center to neighbour cell center, exiting if we reach the max length
+            if ( !growStreamlineFromTo( streamline, currentCell.faceCenter( faceIdx ), currentCell.center(), rate, dominantPhaseOut ) )
+                break;
+
+            // give up if too low rate
+            if ( rate < m_flowThreshold ) return;
+
+            currentCellIdx = currentCell.gridLocalCellIndex();
+        }
+
+        // have we been here before?
+        if ( m_visitedCells.count( currentCellIdx ) > 0 ) return;
+        if ( m_wellCells.count( currentCellIdx ) > 0 ) return;
+
+        m_visitedCells.insert( currentCellIdx );
 
         // find the face with max flow where we should exit the cell
-        cvf::StructGridInterface::FaceType                   exitFace = cvf::StructGridInterface::FaceType::NO_FACE;
+        CellFaceType                   exitFace = cvf::StructGridInterface::FaceType::NO_FACE;
-        std::map<cvf::StructGridInterface::FaceType, double> rateMap;
+        std::map<CellFaceType, double> rateMap;
 
         double maxRate = 0.0;
 
@@ -150,7 +208,7 @@ void RimStreamlineGenerator::growStreamline( RimStreamline* streamline, size_t c
             if ( face == faceIdx ) continue;
 
             RiaDefines::PhaseType tempDominantFace;
-            double                faceRate = m_dataAccess->combinedFaceRate( cell, face, m_phases, direction, tempDominantFace );
+            double                faceRate = m_dataAccess->combinedFaceRate( currentCell, face, m_phases, direction, tempDominantFace );
 
             // if we go backwards from a producer, the rate needs to be flipped
             faceRate *= direction;
@@ -160,27 +218,51 @@ void RimStreamlineGenerator::growStreamline( RimStreamline* streamline, size_t c
                 exitFace         = face;
                 maxRate          = faceRate;
                 dominantPhaseOut = tempDominantFace;
+                index            = currentCellIdx;
             }
 
             rateMap[face] = faceRate;
         }
 
+        // check if we have any NNC connections here
+        int  nncId          = -1;
+        auto nncConnections = nncCandidates( currentCellIdx, m_phases, direction );
+        for ( auto& nncPoint : nncConnections )
+        {
+            RiaDefines::PhaseType tempDominantFace;
+            double nncRate = m_dataAccess->combinedNNCRate( nncPoint.index(), m_phases, direction, tempDominantFace ) * direction;
+            if ( nncRate > maxRate )
+            {
+                nncId            = (int)nncPoint.index();
+                maxRate          = nncRate;
+                dominantPhaseOut = tempDominantFace;
+                exitFace         = CellFaceType::NO_FACE;
+            }
+        }
+
         // did we find an exit?
-        if ( exitFace == cvf::StructGridInterface::FaceType::NO_FACE ) break;
+        if ( exitFace == cvf::StructGridInterface::FaceType::NO_FACE )
+        {
+            if ( nncId < 0 ) return;
+            index = (size_t)nncId;
+
+            for ( auto& nncPoint : nncConnections )
+            {
+                if ( nncPoint.index() == index ) continue;
+                m_seeds.push( nncPoint );
+            }
+        }
 
         // add seeds for other faces with flow > threshold
         for ( auto& kvp : rateMap )
         {
             if ( kvp.first == exitFace ) continue;
 
-            if ( kvp.second >= m_flowThreshold ) m_seeds.push( StreamlineSeedPoint( kvp.second, cell.gridLocalCellIndex(), kvp.first ) );
+            if ( kvp.second >= m_flowThreshold )
+                m_seeds.push( StreamlineSeedPoint( kvp.second, currentCell.gridLocalCellIndex(), kvp.first ) );
         }
 
-        rate = maxRate;
+        rate    = maxRate;
-
-        // grow from cell center to exit face center, stopping if we reach the max point limit
-        if ( !growStreamlineFromTo( streamline, cell.center(), cell.faceCenter( exitFace ), rate, dominantPhaseOut ) ) break;
-
         faceIdx = exitFace;
     }
 }
@@ -218,3 +300,31 @@ bool RimStreamlineGenerator::growStreamlineFromTo( RimStreamline*        streaml
 
     return true;
 }
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+std::list<StreamlineSeedPoint> RimStreamlineGenerator::nncCandidates( size_t cellIdx, std::list<RiaDefines::PhaseType> phases, double direction )
+{
+    std::list<StreamlineSeedPoint> foundCells;
+
+    auto mainGrid = m_dataAccess->grid();
+
+    if ( mainGrid->nncData() == nullptr ) return foundCells;
+
+    auto& connections = mainGrid->nncData()->allConnections();
+    for ( size_t i = 0; i < connections.size(); i++ )
+    {
+        if ( connections[i].c1GlobIdx() == cellIdx )
+        {
+            RiaDefines::PhaseType domPhase;
+            double                rate = m_dataAccess->combinedNNCRate( i, phases, direction, domPhase ) * direction;
+            if ( rate > m_flowThreshold )
+            {
+                foundCells.emplace_back( rate, i, cvf::StructGridInterface::FaceType::NO_FACE );
+            }
+        }
+    }
+
+    return foundCells;
+}

ApplicationLibCode/ProjectDataModel/Streamlines/RimStreamlineGenerator.h

@@ -27,6 +27,7 @@
 #include <list>
 #include <queue>
 #include <set>
+#include <utility>
 
 class RigCell;
 class RimStreamline;
@@ -34,6 +35,9 @@ class StreamlineSeedPoint;
 
 class RimStreamlineGenerator : public RimStreamlineGeneratorBase
 {
+public:
+    using CellFaceType = cvf::StructGridInterface::FaceType;
+
 public:
     RimStreamlineGenerator( std::set<size_t>& wellCells );
     ~RimStreamlineGenerator();
@@ -41,9 +45,11 @@ public:
     void generateTracer( RigCell cell, double direction, QString simWellName, std::list<RimStreamline*>& outStreamlines ) override;
 
 protected:
-    void growStreamline( RimStreamline* streamline, size_t cellIdx, cvf::StructGridInterface::FaceType faceIdx, double direction );
+    void growStreamline( RimStreamline* streamline, size_t index, CellFaceType faceIdx, double direction );
 
     bool growStreamlineFromTo( RimStreamline* streamline, cvf::Vec3d startPos, cvf::Vec3d endpos, double rate, RiaDefines::PhaseType dominantPhase );
 
+    std::list<StreamlineSeedPoint> nncCandidates( size_t cellIdx, std::list<RiaDefines::PhaseType> phases, double direction );
+
     std::priority_queue<StreamlineSeedPoint> m_seeds;
 };

ApplicationLibCode/ProjectDataModel/Streamlines/RimStreamlineInViewCollection.cpp

@@ -436,6 +436,9 @@ void RimStreamlineInViewCollection::updateStreamlines()
     auto eclView = firstAncestorOrThisOfType<RimEclipseView>();
     if ( !eclView ) return;
 
+    auto mainGrid = eclView->mainGrid();
+    if ( !mainGrid ) return;
+
     if ( m_shouldGenerateTracers && isActive() )
     {
         // reset generated streamlines
@@ -452,8 +455,8 @@ void RimStreamlineInViewCollection::updateStreamlines()
 
         // set up the data access helper
         RimStreamlineDataAccess dataAccess;
-        bool                    accessOk =
+
-            dataAccess.setupDataAccess( eclipseCase()->eclipseCaseData()->mainGrid(), eclipseCase()->eclipseCaseData(), phases(), timeIdx );
+        bool accessOk = dataAccess.setupDataAccess( mainGrid, eclipseCase()->eclipseCaseData(), phases(), timeIdx );
 
         // did we find the data we needed?
         if ( accessOk )

ApplicationLibCode/ReservoirDataModel/Well/RigWellTargetCandidatesGenerator.cpp

@@ -25,8 +25,10 @@
 
 #include "RigActiveCellInfo.h"
 #include "RigCaseCellResultsData.h"
+#include "RigEclipseCaseData.h"
 #include "RigEclipseResultAddress.h"
 #include "RigMainGrid.h"
+#include "RigNNCData.h"
 #include "RigStatisticsMath.h"
 
 #include "RimEclipseCase.h"
@@ -69,8 +71,10 @@ void RigWellTargetCandidatesGenerator::generateCandidates( RimEclipseCase*
     auto resultsData = eclipseCase->results( RiaDefines::PorosityModelType::MATRIX_MODEL );
     if ( !resultsData ) return;
 
-    std::vector<double> volume = getVolumeVector( *resultsData, volumeType, volumesType, volumeResultType, timeStepIdx );
+    DataContainer data;
-    if ( volume.empty() )
+
+    data.volume = getVolumeVector( *resultsData, volumeType, volumesType, volumeResultType, timeStepIdx );
+    if ( data.volume.empty() )
     {
         RiaLogging::error( "Unable to produce volume vector." );
         return;
@@ -78,31 +82,41 @@ void RigWellTargetCandidatesGenerator::generateCandidates( RimEclipseCase*
 
     RigEclipseResultAddress pressureAddress( RiaDefines::ResultCatType::DYNAMIC_NATIVE, "PRESSURE" );
     resultsData->ensureKnownResultLoaded( pressureAddress );
-    const std::vector<double>& pressure = resultsData->cellScalarResults( pressureAddress, timeStepIdx );
+    data.pressure = resultsData->cellScalarResults( pressureAddress, timeStepIdx );
 
     RigEclipseResultAddress permeabilityXAddress( RiaDefines::ResultCatType::STATIC_NATIVE, "PERMX" );
     resultsData->ensureKnownResultLoaded( permeabilityXAddress );
-    const std::vector<double>& permeabilityX = resultsData->cellScalarResults( permeabilityXAddress, 0 );
+    data.permeabilityX = resultsData->cellScalarResults( permeabilityXAddress, 0 );
 
     RigEclipseResultAddress permeabilityYAddress( RiaDefines::ResultCatType::STATIC_NATIVE, "PERMY" );
     resultsData->ensureKnownResultLoaded( permeabilityYAddress );
-    const std::vector<double>& permeabilityY = resultsData->cellScalarResults( permeabilityYAddress, 0 );
+    data.permeabilityY = resultsData->cellScalarResults( permeabilityYAddress, 0 );
 
     RigEclipseResultAddress permeabilityZAddress( RiaDefines::ResultCatType::STATIC_NATIVE, "PERMZ" );
     resultsData->ensureKnownResultLoaded( permeabilityZAddress );
-    const std::vector<double>& permeabilityZ = resultsData->cellScalarResults( permeabilityZAddress, 0 );
+    data.permeabilityZ = resultsData->cellScalarResults( permeabilityZAddress, 0 );
 
     RigEclipseResultAddress transmissibilityXAddress( RiaDefines::ResultCatType::STATIC_NATIVE, "TRANX" );
     resultsData->ensureKnownResultLoaded( transmissibilityXAddress );
-    const std::vector<double>& transmissibilityX = resultsData->cellScalarResults( transmissibilityXAddress, 0 );
+    data.transmissibilityX = resultsData->cellScalarResults( transmissibilityXAddress, 0 );
 
     RigEclipseResultAddress transmissibilityYAddress( RiaDefines::ResultCatType::STATIC_NATIVE, "TRANY" );
     resultsData->ensureKnownResultLoaded( transmissibilityYAddress );
-    const std::vector<double>& transmissibilityY = resultsData->cellScalarResults( transmissibilityYAddress, 0 );
+    data.transmissibilityY = resultsData->cellScalarResults( transmissibilityYAddress, 0 );
 
     RigEclipseResultAddress transmissibilityZAddress( RiaDefines::ResultCatType::STATIC_NATIVE, "TRANZ" );
     resultsData->ensureKnownResultLoaded( transmissibilityZAddress );
-    const std::vector<double>& transmissibilityZ = resultsData->cellScalarResults( transmissibilityZAddress, 0 );
+    data.transmissibilityZ = resultsData->cellScalarResults( transmissibilityZAddress, 0 );
+
+    auto mainGrid = eclipseCase->eclipseCaseData()->mainGrid();
+    if ( mainGrid != nullptr )
+    {
+        data.transmissibilityNNC = mainGrid->nncData()->staticConnectionScalarResultByName( RiaDefines::propertyNameCombTrans() );
+    }
+    else
+    {
+        data.transmissibilityNNC = nullptr;
+    }
 
     std::vector<double> filterVector;
     if ( limits.filterAddress.isValid() )
@@ -112,7 +126,7 @@ void RigWellTargetCandidatesGenerator::generateCandidates( RimEclipseCase*
     }
     else
     {
-        filterVector.resize( pressure.size(), std::numeric_limits<double>::infinity() );
+        filterVector.resize( data.pressure.size(), std::numeric_limits<double>::infinity() );
     }
 
     std::vector<int> clusters( activeCellCount.value(), 0 );
@@ -122,19 +136,7 @@ void RigWellTargetCandidatesGenerator::generateCandidates( RimEclipseCase*
     int              numClustersFound = 0;
     for ( int clusterId = 1; clusterId <= numClusters; clusterId++ )
     {
-        std::optional<caf::VecIjk> startCell = findStartCell( eclipseCase,
+        std::optional<caf::VecIjk> startCell = findStartCell( eclipseCase, timeStepIdx, limits, data, filterVector, clusters );
-                                                              timeStepIdx,
-                                                              limits,
-                                                              volume,
-                                                              pressure,
-                                                              permeabilityX,
-                                                              permeabilityY,
-                                                              permeabilityZ,
-                                                              transmissibilityX,
-                                                              transmissibilityY,
-                                                              transmissibilityZ,
-                                                              filterVector,
-                                                              clusters );
 
         if ( startCell.has_value() )
         {
@@ -144,22 +146,7 @@ void RigWellTargetCandidatesGenerator::generateCandidates( RimEclipseCase*
                                   .arg( startCell->j() + 1 )
                                   .arg( startCell->k() + 1 ) );
 
-            growCluster( eclipseCase,
+            growCluster( eclipseCase, startCell.value(), limits, data, filterVector, clusters, clusterId, timeStepIdx, maxIterations );
-                         startCell.value(),
-                         limits,
-                         volume,
-                         pressure,
-                         permeabilityX,
-                         permeabilityY,
-                         permeabilityZ,
-                         transmissibilityX,
-                         transmissibilityY,
-                         transmissibilityZ,
-                         filterVector,
-                         clusters,
-                         clusterId,
-                         timeStepIdx,
-                         maxIterations );
             numClustersFound++;
         }
         else
@@ -191,13 +178,19 @@ void RigWellTargetCandidatesGenerator::generateCandidates( RimEclipseCase*
         }
     }
 
-    std::vector<ClusterStatistics> statistics =
+    std::vector<ClusterStatistics> statistics = generateStatistics( eclipseCase,
-        generateStatistics( eclipseCase, pressure, permeabilityX, permeabilityY, permeabilityZ, numClustersFound, timeStepIdx, resultName );
+                                                                    data.pressure,
-    std::vector<double> totalPorvSoil( clusters.size(), std::numeric_limits<double>::infinity() );
+                                                                    data.permeabilityX,
-    std::vector<double> totalPorvSgas( clusters.size(), std::numeric_limits<double>::infinity() );
+                                                                    data.permeabilityY,
-    std::vector<double> totalPorvSoilAndSgas( clusters.size(), std::numeric_limits<double>::infinity() );
+                                                                    data.permeabilityZ,
-    std::vector<double> totalFipOil( clusters.size(), std::numeric_limits<double>::infinity() );
+                                                                    numClustersFound,
-    std::vector<double> totalFipGas( clusters.size(), std::numeric_limits<double>::infinity() );
+                                                                    timeStepIdx,
+                                                                    resultName );
+    std::vector<double>            totalPorvSoil( clusters.size(), std::numeric_limits<double>::infinity() );
+    std::vector<double>            totalPorvSgas( clusters.size(), std::numeric_limits<double>::infinity() );
+    std::vector<double>            totalPorvSoilAndSgas( clusters.size(), std::numeric_limits<double>::infinity() );
+    std::vector<double>            totalFipOil( clusters.size(), std::numeric_limits<double>::infinity() );
+    std::vector<double>            totalFipGas( clusters.size(), std::numeric_limits<double>::infinity() );
 
     auto addValuesForClusterId = []( std::vector<double>& values, const std::vector<int>& clusters, int clusterId, double value )
     {
@@ -209,7 +202,7 @@ void RigWellTargetCandidatesGenerator::generateCandidates( RimEclipseCase*
     };
 
     int clusterId = 1;
-    for ( auto s : statistics )
+    for ( const auto& s : statistics )
     {
         RiaLogging::info( QString( "Cluster #%1 Statistics" ).arg( s.id ) );
         RiaLogging::info( QString( "Number of cells: %1" ).arg( s.numCells ) );
@@ -252,14 +245,7 @@ void RigWellTargetCandidatesGenerator::generateCandidates( RimEclipseCase*
 std::optional<caf::VecIjk> RigWellTargetCandidatesGenerator::findStartCell( RimEclipseCase*            eclipseCase,
                                                                             size_t                     timeStepIdx,
                                                                             const ClusteringLimits&    limits,
-                                                                            const std::vector<double>& volume,
+                                                                            const DataContainer&       data,
-                                                                            const std::vector<double>& pressure,
-                                                                            const std::vector<double>& permeabilityX,
-                                                                            const std::vector<double>& permeabilityY,
-                                                                            const std::vector<double>& permeabilityZ,
-                                                                            const std::vector<double>& transmissibilityX,
-                                                                            const std::vector<double>& transmissibilityY,
-                                                                            const std::vector<double>& transmissibilityZ,
                                                                             const std::vector<double>& filterVector,
                                                                             const std::vector<int>&    clusters )
 {
@@ -278,12 +264,12 @@ std::optional<caf::VecIjk> RigWellTargetCandidatesGenerator::findStartCell( RimE
         size_t resultIndex = resultsData->activeCellInfo()->cellResultIndex( reservoirCellIdx );
         if ( resultIndex != cvf::UNDEFINED_SIZE_T && clusters[resultIndex] == 0 )
         {
-            const double cellVolume   = volume[resultIndex];
+            const double cellVolume   = data.volume[resultIndex];
-            const double cellPressure = pressure[resultIndex];
+            const double cellPressure = data.pressure[resultIndex];
 
-            const double cellPermeabiltyX = permeabilityX[resultIndex];
+            const double cellPermeabiltyX = data.permeabilityX[resultIndex];
-            const double cellPermeabiltyY = permeabilityY[resultIndex];
+            const double cellPermeabiltyY = data.permeabilityY[resultIndex];
-            const double cellPermeabiltyZ = permeabilityZ[resultIndex];
+            const double cellPermeabiltyZ = data.permeabilityZ[resultIndex];
             const bool permeabilityValidInAnyDirection = ( cellPermeabiltyX >= limits.permeability || cellPermeabiltyY >= limits.permeability ||
                                                            cellPermeabiltyZ >= limits.permeability );
 
@@ -309,14 +295,7 @@ std::optional<caf::VecIjk> RigWellTargetCandidatesGenerator::findStartCell( RimE
 void RigWellTargetCandidatesGenerator::growCluster( RimEclipseCase*            eclipseCase,
                                                     const caf::VecIjk&         startCell,
                                                     const ClusteringLimits&    limits,
-                                                    const std::vector<double>& volume,
+                                                    const DataContainer&       data,
-                                                    const std::vector<double>& pressure,
-                                                    const std::vector<double>& permeabilityX,
-                                                    const std::vector<double>& permeabilityY,
-                                                    const std::vector<double>& permeabilityZ,
-                                                    const std::vector<double>& transmissibilityX,
-                                                    const std::vector<double>& transmissibilityY,
-                                                    const std::vector<double>& transmissibilityZ,
                                                     const std::vector<double>& filterVector,
                                                     std::vector<int>&          clusters,
                                                     int                        clusterId,
@@ -332,19 +311,7 @@ void RigWellTargetCandidatesGenerator::growCluster( RimEclipseCase*            e
 
     for ( int i = 0; i < maxIterations; i++ )
     {
-        foundCells = findCandidates( *eclipseCase,
+        foundCells = findCandidates( eclipseCase, foundCells, limits, data, filterVector, clusters );
-                                     foundCells,
-                                     limits,
-                                     volume,
-                                     pressure,
-                                     permeabilityX,
-                                     permeabilityY,
-                                     permeabilityZ,
-                                     transmissibilityX,
-                                     transmissibilityY,
-                                     transmissibilityZ,
-                                     filterVector,
-                                     clusters );
         if ( foundCells.empty() ) break;
         assignClusterIdToCells( *resultsData->activeCellInfo(), foundCells, clusters, clusterId );
     }
@@ -353,24 +320,18 @@ void RigWellTargetCandidatesGenerator::growCluster( RimEclipseCase*            e
 //--------------------------------------------------------------------------------------------------
 ///
 //--------------------------------------------------------------------------------------------------
-std::vector<size_t> RigWellTargetCandidatesGenerator::findCandidates( const RimEclipseCase&      eclipseCase,
+std::vector<size_t> RigWellTargetCandidatesGenerator::findCandidates( RimEclipseCase*            eclipseCase,
                                                                       const std::vector<size_t>& previousCells,
                                                                       const ClusteringLimits&    limits,
-                                                                      const std::vector<double>& volume,
+                                                                      const DataContainer&       data,
-                                                                      const std::vector<double>& pressure,
-                                                                      const std::vector<double>& permeabilityX,
-                                                                      const std::vector<double>& permeabilityY,
-                                                                      const std::vector<double>& permeabilityZ,
-                                                                      const std::vector<double>& transmissibilityX,
-                                                                      const std::vector<double>& transmissibilityY,
-                                                                      const std::vector<double>& transmissibilityZ,
                                                                       const std::vector<double>& filterVector,
                                                                       std::vector<int>&          clusters )
 {
     std::vector<size_t> candidates;
-    auto                resultsData = eclipseCase.results( RiaDefines::PorosityModelType::MATRIX_MODEL );
+    auto                resultsData = eclipseCase->results( RiaDefines::PorosityModelType::MATRIX_MODEL );
+    auto                mainGrid    = eclipseCase->eclipseCaseData()->mainGrid();
 
-    const std::vector<cvf::StructGridInterface::FaceType> faces = {
+    const std::vector<CellFaceType> faces = {
         cvf::StructGridInterface::FaceType::POS_I,
         cvf::StructGridInterface::FaceType::NEG_I,
         cvf::StructGridInterface::FaceType::POS_J,
@@ -381,7 +342,7 @@ std::vector<size_t> RigWellTargetCandidatesGenerator::findCandidates( const RimE
 
     for ( size_t cellIdx : previousCells )
     {
-        const RigCell& cell = eclipseCase.mainGrid()->cell( cellIdx );
+        const RigCell& cell = mainGrid->cell( cellIdx );
         if ( cell.isInvalid() ) continue;
 
         RigGridBase* grid               = cell.hostGrid();
@@ -392,7 +353,7 @@ std::vector<size_t> RigWellTargetCandidatesGenerator::findCandidates( const RimE
 
         grid->ijkFromCellIndex( gridLocalCellIndex, &i, &j, &k );
 
-        for ( cvf::StructGridInterface::FaceType face : faces )
+        for ( CellFaceType face : faces )
         {
             size_t gridLocalNeighborCellIdx;
             if ( grid->cellIJKNeighbor( i, j, k, face, &gridLocalNeighborCellIdx ) )
@@ -401,16 +362,17 @@ std::vector<size_t> RigWellTargetCandidatesGenerator::findCandidates( const RimE
                 size_t neighborResultIndex = resultsData->activeCellInfo()->cellResultIndex( neighborResvCellIdx );
                 if ( neighborResultIndex != cvf::UNDEFINED_SIZE_T && clusters[neighborResultIndex] == 0 )
                 {
-                    double permeability     = getValueForFace( permeabilityX, permeabilityY, permeabilityZ, face, neighborResultIndex );
+                    double permeability =
-                    double transmissibility = getTransmissibilityValueForFace( transmissibilityX,
+                        getValueForFace( data.permeabilityX, data.permeabilityY, data.permeabilityZ, face, neighborResultIndex );
-                                                                               transmissibilityY,
+                    double transmissibility = getTransmissibilityValueForFace( data.transmissibilityX,
-                                                                               transmissibilityZ,
+                                                                               data.transmissibilityY,
+                                                                               data.transmissibilityZ,
                                                                                face,
                                                                                resultIndex,
                                                                                neighborResultIndex );
                     bool   filterValue      = !std::isinf( filterVector[neighborResultIndex] ) && filterVector[neighborResultIndex] > 0.0;
 
-                    if ( volume[neighborResultIndex] > limits.volume && pressure[neighborResultIndex] > limits.pressure &&
+                    if ( data.volume[neighborResultIndex] > limits.volume && data.pressure[neighborResultIndex] > limits.pressure &&
                          permeability > limits.permeability && transmissibility > limits.transmissibility && filterValue )
                     {
                         candidates.push_back( neighborResvCellIdx );
@@ -419,6 +381,29 @@ std::vector<size_t> RigWellTargetCandidatesGenerator::findCandidates( const RimE
                 }
             }
         }
+
+        if ( data.transmissibilityNNC != nullptr )
+        {
+            auto nncCells = nncConnectionCellAndResult( cellIdx, mainGrid );
+            for ( auto& [cellInfo, nncResultIdx] : nncCells )
+            {
+                auto& [otherCellIdx, face] = cellInfo;
+                double transmissibility    = data.transmissibilityNNC->at( nncResultIdx );
+
+                size_t otherResultIndex = resultsData->activeCellInfo()->cellResultIndex( otherCellIdx );
+
+                double permeability = getValueForFace( data.permeabilityX, data.permeabilityY, data.permeabilityZ, face, otherResultIndex );
+
+                bool filterValue = !std::isinf( filterVector[otherResultIndex] ) && filterVector[otherResultIndex] > 0.0;
+
+                if ( data.volume[otherResultIndex] > limits.volume && data.pressure[otherResultIndex] > limits.pressure &&
+                     permeability > limits.permeability && transmissibility > limits.transmissibility && filterValue )
+                {
+                    candidates.push_back( otherCellIdx );
+                    clusters[otherResultIndex] = -1;
+                }
+            }
+        }
     }
 
     return candidates;
@@ -531,11 +516,11 @@ std::optional<size_t> RigWellTargetCandidatesGenerator::getActiveCellCount( RimE
 //--------------------------------------------------------------------------------------------------
 ///
 //--------------------------------------------------------------------------------------------------
-double RigWellTargetCandidatesGenerator::getValueForFace( const std::vector<double>&         x,
+double RigWellTargetCandidatesGenerator::getValueForFace( const std::vector<double>& x,
-                                                          const std::vector<double>&         y,
+                                                          const std::vector<double>& y,
-                                                          const std::vector<double>&         z,
+                                                          const std::vector<double>& z,
-                                                          cvf::StructGridInterface::FaceType face,
+                                                          CellFaceType               face,
-                                                          size_t                             resultIndex )
+                                                          size_t                     resultIndex )
 {
     if ( face == cvf::StructGridInterface::FaceType::POS_I || face == cvf::StructGridInterface::FaceType::NEG_I ) return x[resultIndex];
     if ( face == cvf::StructGridInterface::FaceType::POS_J || face == cvf::StructGridInterface::FaceType::NEG_J ) return y[resultIndex];
@@ -546,12 +531,12 @@ double RigWellTargetCandidatesGenerator::getValueForFace( const std::vector<doub
 //--------------------------------------------------------------------------------------------------
 ///
 //--------------------------------------------------------------------------------------------------
-double RigWellTargetCandidatesGenerator::getTransmissibilityValueForFace( const std::vector<double>&         x,
+double RigWellTargetCandidatesGenerator::getTransmissibilityValueForFace( const std::vector<double>& x,
-                                                                          const std::vector<double>&         y,
+                                                                          const std::vector<double>& y,
-                                                                          const std::vector<double>&         z,
+                                                                          const std::vector<double>& z,
-                                                                          cvf::StructGridInterface::FaceType face,
+                                                                          CellFaceType               face,
-                                                                          size_t                             resultIndex,
+                                                                          size_t                     resultIndex,
-                                                                          size_t                             neighborResultIndex )
+                                                                          size_t                     neighborResultIndex )
 {
     // For negative directions (NEG_I, NEG_J, NEG_K) use the value from the neighbor cell
     bool isPos = face == cvf::StructGridInterface::FaceType::POS_I || face == cvf::StructGridInterface::FaceType::POS_J ||
@@ -959,3 +944,25 @@ cvf::BoundingBox RigWellTargetCandidatesGenerator::computeBoundingBoxForResult(
 
     return boundingBox;
 }
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+std::list<std::pair<std::pair<size_t, RigWellTargetCandidatesGenerator::CellFaceType>, size_t>>
+    RigWellTargetCandidatesGenerator::nncConnectionCellAndResult( size_t cellIdx, RigMainGrid* mainGrid )
+{
+    std::list<std::pair<std::pair<size_t, CellFaceType>, size_t>> foundCells;
+
+    if ( mainGrid->nncData() == nullptr ) return foundCells;
+
+    auto& connections = mainGrid->nncData()->allConnections();
+    for ( size_t i = 0; i < connections.size(); i++ )
+    {
+        if ( connections[i].c1GlobIdx() == cellIdx )
+        {
+            foundCells.push_back( std::make_pair( std::make_pair( connections[i].c2GlobIdx(), connections[i].face() ), i ) );
+        }
+    }
+
+    return foundCells;
+}

ApplicationLibCode/ReservoirDataModel/Well/RigWellTargetCandidatesGenerator.h

@@ -25,11 +25,14 @@
 
 #include "RigEclipseResultAddress.h"
 
+#include <list>
 #include <map>
 #include <optional>
+#include <utility>
 
 class RigActiveCellInfo;
 class RigCaseCellResultsData;
+class RigMainGrid;
 class RimEclipseCase;
 class RimEclipseCaseEnsemble;
 class RimRegularGridCase;
@@ -41,6 +44,8 @@ class RimRegularGridCase;
 class RigWellTargetCandidatesGenerator
 {
 public:
+    using CellFaceType = cvf::StructGridInterface::FaceType;
+
     enum class VolumeType
     {
         OIL,
@@ -72,6 +77,20 @@ public:
         RigEclipseResultAddress filterAddress;
     };
 
+    struct DataContainer
+    {
+        std::vector<double>        volume;
+        std::vector<double>        pressure;
+        std::vector<double>        permeabilityX;
+        std::vector<double>        permeabilityY;
+        std::vector<double>        permeabilityZ;
+        std::vector<double>        permeabilityNNC;
+        std::vector<double>        transmissibilityX;
+        std::vector<double>        transmissibilityY;
+        std::vector<double>        transmissibilityZ;
+        const std::vector<double>* transmissibilityNNC;
+    };
+
     static void generateCandidates( RimEclipseCase*         eclipseCase,
                                     size_t                  timeStepIdx,
                                     VolumeType              volumeType,
@@ -124,45 +143,24 @@ private:
     static std::optional<caf::VecIjk> findStartCell( RimEclipseCase*            eclipseCase,
                                                      size_t                     timeStepIdx,
                                                      const ClusteringLimits&    limits,
-                                                     const std::vector<double>& volume,
+                                                     const DataContainer&       data,
-                                                     const std::vector<double>& pressure,
-                                                     const std::vector<double>& permeabilityX,
-                                                     const std::vector<double>& permeabilityY,
-                                                     const std::vector<double>& permeabilityZ,
-                                                     const std::vector<double>& transmissibilityX,
-                                                     const std::vector<double>& transmissibilityY,
-                                                     const std::vector<double>& transmissibilityZ,
                                                      const std::vector<double>& filterVector,
                                                      const std::vector<int>&    clusters );
 
     static void growCluster( RimEclipseCase*            eclipseCase,
                              const caf::VecIjk&         startCell,
                              const ClusteringLimits&    limits,
-                             const std::vector<double>& volume,
+                             const DataContainer&       data,
-                             const std::vector<double>& pressure,
-                             const std::vector<double>& permeabilityX,
-                             const std::vector<double>& permeabilityY,
-                             const std::vector<double>& permeabilityZ,
-                             const std::vector<double>& transmissibilityX,
-                             const std::vector<double>& transmissibilityY,
-                             const std::vector<double>& transmissibilityZ,
                              const std::vector<double>& filterVector,
                              std::vector<int>&          clusters,
                              int                        clusterId,
                              size_t                     timeStepIdx,
                              int                        maxIterations );
 
-    static std::vector<size_t> findCandidates( const RimEclipseCase&      eclipseCase,
+    static std::vector<size_t> findCandidates( RimEclipseCase*            eclipseCase,
                                                const std::vector<size_t>& previousCells,
                                                const ClusteringLimits&    limits,
-                                               const std::vector<double>& volume,
+                                               const DataContainer&       data,
-                                               const std::vector<double>& pressure,
-                                               const std::vector<double>& permeabilityX,
-                                               const std::vector<double>& permeabilityY,
-                                               const std::vector<double>& permeabilityZ,
-                                               const std::vector<double>& transmissibilityX,
-                                               const std::vector<double>& transmissibilityY,
-                                               const std::vector<double>& transmissibilityZ,
                                                const std::vector<double>& filterVector,
                                                std::vector<int>&          clusters );
 
@@ -179,18 +177,18 @@ private:
 
     static void createResultVector( RimEclipseCase& eclipseCase, const QString& resultName, const std::vector<int>& clusterIds, double value );
 
-    static double getValueForFace( const std::vector<double>&         x,
+    static double getValueForFace( const std::vector<double>& x,
-                                   const std::vector<double>&         y,
+                                   const std::vector<double>& y,
-                                   const std::vector<double>&         z,
+                                   const std::vector<double>& z,
-                                   cvf::StructGridInterface::FaceType face,
+                                   CellFaceType               face,
-                                   size_t                             resultIndex );
+                                   size_t                     resultIndex );
 
-    static double getTransmissibilityValueForFace( const std::vector<double>&         x,
+    static double getTransmissibilityValueForFace( const std::vector<double>& x,
-                                                   const std::vector<double>&         y,
+                                                   const std::vector<double>& y,
-                                                   const std::vector<double>&         z,
+                                                   const std::vector<double>& z,
-                                                   cvf::StructGridInterface::FaceType face,
+                                                   CellFaceType               face,
-                                                   size_t                             resultIndex,
+                                                   size_t                     resultIndex,
-                                                   size_t                             neighborResultIndex );
+                                                   size_t                     neighborResultIndex );
 
     static std::vector<RigWellTargetCandidatesGenerator::ClusterStatistics> generateStatistics( RimEclipseCase*            eclipseCase,
                                                                                                 const std::vector<double>& pressure,
@@ -211,4 +209,6 @@ private:
                                                 std::vector<int>&                                    occupancy );
 
     static cvf::BoundingBox computeBoundingBoxForResult( RimEclipseCase& eclipseCase, const QString& resultName, size_t timeStepIndex );
+
+    static std::list<std::pair<std::pair<size_t, CellFaceType>, size_t>> nncConnectionCellAndResult( size_t cellIdx, RigMainGrid* mainGrid );
 };