mirror of
https://github.com/OPM/ResInsight.git
synced 2025-02-25 18:55:39 -06:00
#5915 improve performance of NNC computation and limit to active cells
This commit is contained in:
Gaute Lindkvist
@@ -18,6 +18,7 @@
|
||||
|
||||
#include "RigCellFaceGeometryTools.h"
|
||||
|
||||
#include "RigActiveCellInfo.h"
|
||||
#include "RigCell.h"
|
||||
#include "RigMainGrid.h"
|
||||
#include "RigNncConnection.h"
|
||||
@@ -28,6 +29,8 @@
|
||||
|
||||
#include <QDebug>
|
||||
|
||||
#include <omp.h>
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
@@ -40,7 +43,7 @@ cvf::StructGridInterface::FaceType
|
||||
{
|
||||
// Try to find the shared face
|
||||
|
||||
bool isPossibleNeighborInDirection[6] = {true, true, true, true, true, true};
|
||||
bool isPossibleNeighborInDirection[6] = { true, true, true, true, true, true };
|
||||
|
||||
if ( c1.hostGrid() == c2.hostGrid() )
|
||||
{
|
||||
@@ -116,55 +119,42 @@ cvf::StructGridInterface::FaceType
|
||||
return cvf::StructGridInterface::NO_FACE;
|
||||
}
|
||||
|
||||
void assignThreadConnections( std::set<std::pair<size_t, size_t>>& existingPairs,
|
||||
RigConnectionContainer& allConnections,
|
||||
RigConnectionContainer& threadConnections )
|
||||
{
|
||||
for ( size_t i = 0; i < threadConnections.size(); ++i )
|
||||
{
|
||||
#pragma omp critical
|
||||
{
|
||||
RigConnection connection = threadConnections[i];
|
||||
auto it = existingPairs.emplace( connection.m_c1GlobIdx, connection.m_c2GlobIdx );
|
||||
if ( it.second )
|
||||
{
|
||||
allConnections.push_back( connection );
|
||||
}
|
||||
}
|
||||
}
|
||||
threadConnections.clear();
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
std::vector<RigConnection> RigCellFaceGeometryTools::computeOtherNncs( const RigMainGrid* mainGrid,
|
||||
const std::vector<RigConnection>& nativeConnections )
|
||||
RigConnectionContainer RigCellFaceGeometryTools::computeOtherNncs( const RigMainGrid* mainGrid,
|
||||
const RigConnectionContainer& nativeConnections,
|
||||
const RigActiveCellInfo* activeCellInfo )
|
||||
{
|
||||
// Compute Non-Neighbor Connections (NNC) not reported by Eclipse. NNCs with zero transmissibility are not reported
|
||||
// by Eclipse. Use faults as basis for subset of cells to find NNC connection for. The imported connections from
|
||||
// Eclipse are located at the beginning of the connections vector.
|
||||
|
||||
std::vector<RigConnection> otherConnections;
|
||||
std::set<std::pair<size_t, size_t>> nativeCellPairs;
|
||||
|
||||
class CellPair
|
||||
for ( size_t i = 0; i < nativeConnections.size(); ++i )
|
||||
{
|
||||
public:
|
||||
CellPair( size_t globalIdx1, size_t globalIdx2 )
|
||||
{
|
||||
if ( globalIdx1 < globalIdx2 )
|
||||
{
|
||||
m_globalCellIdx1 = globalIdx1;
|
||||
m_globalCellIdx2 = globalIdx2;
|
||||
}
|
||||
else
|
||||
{
|
||||
m_globalCellIdx1 = globalIdx2;
|
||||
m_globalCellIdx2 = globalIdx1;
|
||||
}
|
||||
}
|
||||
|
||||
bool operator<( const CellPair& other ) const
|
||||
{
|
||||
if ( m_globalCellIdx1 != other.m_globalCellIdx1 )
|
||||
{
|
||||
return m_globalCellIdx1 < other.m_globalCellIdx1;
|
||||
}
|
||||
|
||||
return ( m_globalCellIdx2 < other.m_globalCellIdx2 );
|
||||
}
|
||||
|
||||
private:
|
||||
size_t m_globalCellIdx1;
|
||||
size_t m_globalCellIdx2;
|
||||
};
|
||||
|
||||
std::set<CellPair> nativeCellPairs;
|
||||
|
||||
for ( const auto& c : nativeConnections )
|
||||
{
|
||||
nativeCellPairs.emplace( CellPair( c.m_c1GlobIdx, c.m_c2GlobIdx ) );
|
||||
RigConnection c = nativeConnections[i];
|
||||
nativeCellPairs.emplace( c.m_c1GlobIdx, c.m_c2GlobIdx );
|
||||
}
|
||||
|
||||
if ( nativeConnections.size() != nativeCellPairs.size() )
|
||||
@@ -175,164 +165,181 @@ std::vector<RigConnection> RigCellFaceGeometryTools::computeOtherNncs( const Rig
|
||||
qDebug() << message;
|
||||
}
|
||||
|
||||
std::set<CellPair> otherCellPairs;
|
||||
|
||||
const cvf::Collection<RigFault>& faults = mainGrid->faults();
|
||||
|
||||
for ( size_t faultIdx = 0; faultIdx < faults.size(); faultIdx++ )
|
||||
std::set<std::pair<size_t, size_t>> existingPairs;
|
||||
RigConnectionContainer otherConnections;
|
||||
|
||||
for ( int faultIdx = 0; faultIdx < (int)faults.size(); faultIdx++ )
|
||||
{
|
||||
const RigFault* fault = faults.at( faultIdx );
|
||||
|
||||
const std::vector<RigFault::FaultFace>& faultFaces = fault->faultFaces();
|
||||
|
||||
#pragma omp parallel for
|
||||
for ( int faceIdx = 0; faceIdx < static_cast<int>( faultFaces.size() ); faceIdx++ )
|
||||
#pragma omp parallel
|
||||
{
|
||||
const RigFault::FaultFace& f = faultFaces[faceIdx];
|
||||
|
||||
size_t sourceReservoirCellIndex = f.m_nativeReservoirCellIndex;
|
||||
cvf::StructGridInterface::FaceType sourceCellFace = f.m_nativeFace;
|
||||
|
||||
if ( sourceReservoirCellIndex >= mainGrid->cellCount() )
|
||||
RigConnectionContainer threadConnections;
|
||||
#pragma omp for
|
||||
for ( int faceIdx = 0; faceIdx < static_cast<int>( faultFaces.size() ); faceIdx++ )
|
||||
{
|
||||
continue;
|
||||
}
|
||||
const RigFault::FaultFace& f = faultFaces[faceIdx];
|
||||
|
||||
const std::vector<cvf::Vec3d>& mainGridNodes = mainGrid->nodes();
|
||||
|
||||
cvf::BoundingBox bb;
|
||||
std::array<size_t, 4> sourceFaceIndices;
|
||||
mainGrid->globalCellArray()[sourceReservoirCellIndex].faceIndices( sourceCellFace, &sourceFaceIndices );
|
||||
|
||||
bb.add( mainGridNodes[sourceFaceIndices[0]] );
|
||||
bb.add( mainGridNodes[sourceFaceIndices[1]] );
|
||||
bb.add( mainGridNodes[sourceFaceIndices[2]] );
|
||||
bb.add( mainGridNodes[sourceFaceIndices[3]] );
|
||||
|
||||
std::vector<size_t> closeCells;
|
||||
mainGrid->findIntersectingCells( bb, &closeCells );
|
||||
|
||||
cvf::StructGridInterface::FaceType candidateFace = cvf::StructGridInterface::oppositeFace( sourceCellFace );
|
||||
|
||||
size_t neighborCellIndex = std::numeric_limits<size_t>::max();
|
||||
size_t ni = std::numeric_limits<size_t>::max();
|
||||
size_t nj = std::numeric_limits<size_t>::max();
|
||||
size_t nk = std::numeric_limits<size_t>::max();
|
||||
|
||||
{
|
||||
size_t i;
|
||||
size_t j;
|
||||
size_t k;
|
||||
mainGrid->ijkFromCellIndex( sourceReservoirCellIndex, &i, &j, &k );
|
||||
|
||||
mainGrid->neighborIJKAtCellFace( i, j, k, sourceCellFace, &ni, &nj, &nk );
|
||||
|
||||
if ( mainGrid->isCellValid( ni, nj, nk ) )
|
||||
bool atLeastOneCellActive = !activeCellInfo ||
|
||||
( activeCellInfo->isActive( f.m_nativeReservoirCellIndex ) ||
|
||||
activeCellInfo->isActive( f.m_oppositeReservoirCellIndex ) );
|
||||
if ( atLeastOneCellActive )
|
||||
{
|
||||
neighborCellIndex = mainGrid->cellIndexFromIJK( ni, nj, nk );
|
||||
RigConnectionContainer faceConnections = extractConnectionsForFace( f, mainGrid, nativeCellPairs );
|
||||
threadConnections.insert( faceConnections );
|
||||
}
|
||||
}
|
||||
|
||||
for ( size_t candidateCellIndex : closeCells )
|
||||
{
|
||||
if ( candidateCellIndex == sourceReservoirCellIndex )
|
||||
{
|
||||
// Exclude cellIndex for source cell
|
||||
continue;
|
||||
}
|
||||
|
||||
if ( candidateCellIndex == neighborCellIndex )
|
||||
{
|
||||
// Exclude direct neighbor
|
||||
continue;
|
||||
}
|
||||
|
||||
if ( candidateCellIndex >= mainGrid->cellCount() )
|
||||
{
|
||||
continue;
|
||||
}
|
||||
|
||||
if ( neighborCellIndex != std::numeric_limits<size_t>::max() )
|
||||
{
|
||||
// Find target IJK index based on source cell and cell face
|
||||
// Exclude cells not matching destination target index
|
||||
|
||||
size_t ci = std::numeric_limits<size_t>::max();
|
||||
size_t cj = std::numeric_limits<size_t>::max();
|
||||
size_t ck = std::numeric_limits<size_t>::max();
|
||||
mainGrid->ijkFromCellIndex( candidateCellIndex, &ci, &cj, &ck );
|
||||
|
||||
auto gridAxis = cvf::StructGridInterface::gridAxisFromFace( sourceCellFace );
|
||||
if ( gridAxis == cvf::StructGridInterface::GridAxisType::AXIS_I )
|
||||
{
|
||||
if ( ni != ci )
|
||||
{
|
||||
continue;
|
||||
}
|
||||
}
|
||||
else if ( gridAxis == cvf::StructGridInterface::GridAxisType::AXIS_J )
|
||||
{
|
||||
if ( nj != cj )
|
||||
{
|
||||
continue;
|
||||
}
|
||||
}
|
||||
else if ( gridAxis == cvf::StructGridInterface::GridAxisType::AXIS_K )
|
||||
{
|
||||
if ( nk != ck )
|
||||
{
|
||||
continue;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
CellPair candidate( sourceReservoirCellIndex, candidateCellIndex );
|
||||
|
||||
if ( nativeCellPairs.count( candidate ) > 0 )
|
||||
{
|
||||
continue;
|
||||
}
|
||||
|
||||
std::vector<size_t> polygon;
|
||||
std::vector<cvf::Vec3d> intersections;
|
||||
|
||||
std::array<size_t, 4> candidateFaceIndices;
|
||||
mainGrid->globalCellArray()[candidateCellIndex].faceIndices( candidateFace, &candidateFaceIndices );
|
||||
|
||||
bool foundOverlap =
|
||||
cvf::GeometryTools::calculateOverlapPolygonOfTwoQuads( &polygon,
|
||||
&intersections,
|
||||
(cvf::EdgeIntersectStorage<size_t>*)nullptr,
|
||||
cvf::wrapArrayConst( &mainGridNodes ),
|
||||
sourceFaceIndices.data(),
|
||||
candidateFaceIndices.data(),
|
||||
1e-6 );
|
||||
|
||||
if ( foundOverlap )
|
||||
{
|
||||
RigConnection conn;
|
||||
conn.m_c1GlobIdx = sourceReservoirCellIndex;
|
||||
conn.m_c1Face = sourceCellFace;
|
||||
conn.m_c2GlobIdx = candidateCellIndex;
|
||||
|
||||
conn.m_polygon = RigCellFaceGeometryTools::extractPolygon( mainGridNodes, polygon, intersections );
|
||||
|
||||
#pragma omp critical( critical_section_nnc_computations )
|
||||
{
|
||||
auto itBoolPair = otherCellPairs.insert( candidate );
|
||||
if ( itBoolPair.second )
|
||||
{
|
||||
otherConnections.emplace_back( conn );
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
// Merge together connections per thread
|
||||
assignThreadConnections( existingPairs, otherConnections, threadConnections );
|
||||
} // end parallel region
|
||||
}
|
||||
|
||||
return otherConnections;
|
||||
}
|
||||
|
||||
RigConnectionContainer
|
||||
RigCellFaceGeometryTools::extractConnectionsForFace( const RigFault::FaultFace& face,
|
||||
const RigMainGrid* mainGrid,
|
||||
const std::set<std::pair<size_t, size_t>>& nativeCellPairs )
|
||||
{
|
||||
RigConnectionContainer faceConnections;
|
||||
size_t sourceReservoirCellIndex = face.m_nativeReservoirCellIndex;
|
||||
cvf::StructGridInterface::FaceType sourceCellFace = face.m_nativeFace;
|
||||
|
||||
if ( sourceReservoirCellIndex >= mainGrid->cellCount() )
|
||||
{
|
||||
return {};
|
||||
}
|
||||
|
||||
const std::vector<cvf::Vec3d>& mainGridNodes = mainGrid->nodes();
|
||||
|
||||
cvf::BoundingBox bb;
|
||||
std::array<size_t, 4> sourceFaceIndices;
|
||||
mainGrid->globalCellArray()[sourceReservoirCellIndex].faceIndices( sourceCellFace, &sourceFaceIndices );
|
||||
|
||||
bb.add( mainGridNodes[sourceFaceIndices[0]] );
|
||||
bb.add( mainGridNodes[sourceFaceIndices[1]] );
|
||||
bb.add( mainGridNodes[sourceFaceIndices[2]] );
|
||||
bb.add( mainGridNodes[sourceFaceIndices[3]] );
|
||||
|
||||
std::vector<size_t> closeCells;
|
||||
mainGrid->findIntersectingCells( bb, &closeCells );
|
||||
|
||||
cvf::StructGridInterface::FaceType candidateFace = cvf::StructGridInterface::oppositeFace( sourceCellFace );
|
||||
|
||||
size_t neighborCellIndex = std::numeric_limits<size_t>::max();
|
||||
size_t ni = std::numeric_limits<size_t>::max();
|
||||
size_t nj = std::numeric_limits<size_t>::max();
|
||||
size_t nk = std::numeric_limits<size_t>::max();
|
||||
|
||||
{
|
||||
size_t i;
|
||||
size_t j;
|
||||
size_t k;
|
||||
mainGrid->ijkFromCellIndex( sourceReservoirCellIndex, &i, &j, &k );
|
||||
|
||||
mainGrid->neighborIJKAtCellFace( i, j, k, sourceCellFace, &ni, &nj, &nk );
|
||||
|
||||
if ( mainGrid->isCellValid( ni, nj, nk ) )
|
||||
{
|
||||
neighborCellIndex = mainGrid->cellIndexFromIJK( ni, nj, nk );
|
||||
}
|
||||
}
|
||||
|
||||
for ( size_t candidateCellIndex : closeCells )
|
||||
{
|
||||
if ( candidateCellIndex == sourceReservoirCellIndex )
|
||||
{
|
||||
// Exclude cellIndex for source cell
|
||||
continue;
|
||||
}
|
||||
|
||||
if ( candidateCellIndex == neighborCellIndex )
|
||||
{
|
||||
// Exclude direct neighbor
|
||||
continue;
|
||||
}
|
||||
|
||||
if ( candidateCellIndex >= mainGrid->cellCount() )
|
||||
{
|
||||
continue;
|
||||
}
|
||||
|
||||
if ( neighborCellIndex != std::numeric_limits<size_t>::max() )
|
||||
{
|
||||
// Find target IJK index based on source cell and cell face
|
||||
// Exclude cells not matching destination target index
|
||||
|
||||
size_t ci = std::numeric_limits<size_t>::max();
|
||||
size_t cj = std::numeric_limits<size_t>::max();
|
||||
size_t ck = std::numeric_limits<size_t>::max();
|
||||
mainGrid->ijkFromCellIndex( candidateCellIndex, &ci, &cj, &ck );
|
||||
|
||||
auto gridAxis = cvf::StructGridInterface::gridAxisFromFace( sourceCellFace );
|
||||
if ( gridAxis == cvf::StructGridInterface::GridAxisType::AXIS_I )
|
||||
{
|
||||
if ( ni != ci )
|
||||
{
|
||||
continue;
|
||||
}
|
||||
}
|
||||
else if ( gridAxis == cvf::StructGridInterface::GridAxisType::AXIS_J )
|
||||
{
|
||||
if ( nj != cj )
|
||||
{
|
||||
continue;
|
||||
}
|
||||
}
|
||||
else if ( gridAxis == cvf::StructGridInterface::GridAxisType::AXIS_K )
|
||||
{
|
||||
if ( nk != ck )
|
||||
{
|
||||
continue;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
std::pair<size_t, size_t> candidate( sourceReservoirCellIndex, candidateCellIndex );
|
||||
|
||||
if ( nativeCellPairs.count( candidate ) > 0 )
|
||||
{
|
||||
continue;
|
||||
}
|
||||
|
||||
std::vector<size_t> polygon;
|
||||
std::vector<cvf::Vec3d> intersections;
|
||||
|
||||
std::array<size_t, 4> candidateFaceIndices;
|
||||
mainGrid->globalCellArray()[candidateCellIndex].faceIndices( candidateFace, &candidateFaceIndices );
|
||||
|
||||
bool foundOverlap =
|
||||
cvf::GeometryTools::calculateOverlapPolygonOfTwoQuads( &polygon,
|
||||
&intersections,
|
||||
(cvf::EdgeIntersectStorage<size_t>*)nullptr,
|
||||
cvf::wrapArrayConst( &mainGridNodes ),
|
||||
sourceFaceIndices.data(),
|
||||
candidateFaceIndices.data(),
|
||||
1e-6 );
|
||||
|
||||
if ( foundOverlap )
|
||||
{
|
||||
RigConnection conn;
|
||||
conn.m_c1GlobIdx = sourceReservoirCellIndex;
|
||||
conn.m_c2GlobIdx = candidateCellIndex;
|
||||
conn.m_c1Face = sourceCellFace;
|
||||
|
||||
conn.m_polygon = RigCellFaceGeometryTools::extractPolygon( mainGridNodes, polygon, intersections );
|
||||
|
||||
faceConnections.push_back( conn );
|
||||
}
|
||||
}
|
||||
return faceConnections;
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
|
||||
Reference in New Issue
Block a user