mirror of
https://github.com/OPM/ResInsight.git
synced 2026-10-26 07:00:17 -05:00
Magne Sjaastad
c04dbbeec3
Computing additional NNCs spent most of its time running quad overlap tests on candidate cells outside the neighbor pillar column. Most of these touch the fault face along an edge only. Restrict candidates for I and J faces to the neighbor column. This intentionally skips a small number of non-zero area overlaps outside the column.
369 lines
17 KiB
C++
369 lines
17 KiB
C++
/////////////////////////////////////////////////////////////////////////////////
|
|
//
|
|
// 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 "RigCellFaceGeometryTools.h"
|
|
|
|
#include "RiaLogging.h"
|
|
|
|
#include "RigActiveCellInfo.h"
|
|
#include "RigCell.h"
|
|
#include "RigMainGrid.h"
|
|
#include "RigNncConnection.h"
|
|
|
|
#include "cvfGeometryTools.h"
|
|
|
|
#include "cafAssert.h"
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
cvf::StructGridInterface::FaceType RigCellFaceGeometryTools::calculateCellFaceOverlap( const RigCell& c1,
|
|
const RigCell& c2,
|
|
const RigMainGrid& mainGrid,
|
|
std::vector<size_t>* connectionPolygon,
|
|
std::vector<cvf::Vec3d>* connectionIntersections )
|
|
{
|
|
// Try to find the shared face
|
|
|
|
bool isPossibleNeighborInDirection[6] = { true, true, true, true, true, true };
|
|
|
|
if ( c1.hostGrid() == c2.hostGrid() )
|
|
{
|
|
char hasNeighbourInAnyDirection = 0;
|
|
|
|
size_t i1, j1, k1;
|
|
c1.hostGrid()->ijkFromCellIndex( c1.gridLocalCellIndex(), &i1, &j1, &k1 );
|
|
size_t i2, j2, k2;
|
|
c2.hostGrid()->ijkFromCellIndex( c2.gridLocalCellIndex(), &i2, &j2, &k2 );
|
|
|
|
isPossibleNeighborInDirection[cvf::StructGridInterface::POS_I] = ( ( i1 + 1 ) == i2 );
|
|
isPossibleNeighborInDirection[cvf::StructGridInterface::NEG_I] = ( ( i2 + 1 ) == i1 );
|
|
isPossibleNeighborInDirection[cvf::StructGridInterface::POS_J] = ( ( j1 + 1 ) == j2 );
|
|
isPossibleNeighborInDirection[cvf::StructGridInterface::NEG_J] = ( ( j2 + 1 ) == j1 );
|
|
isPossibleNeighborInDirection[cvf::StructGridInterface::POS_K] = ( ( k1 + 1 ) == k2 );
|
|
isPossibleNeighborInDirection[cvf::StructGridInterface::NEG_K] = ( ( k2 + 1 ) == k1 );
|
|
|
|
hasNeighbourInAnyDirection =
|
|
isPossibleNeighborInDirection[cvf::StructGridInterface::POS_I] + isPossibleNeighborInDirection[cvf::StructGridInterface::NEG_I] +
|
|
isPossibleNeighborInDirection[cvf::StructGridInterface::POS_J] + isPossibleNeighborInDirection[cvf::StructGridInterface::NEG_J] +
|
|
isPossibleNeighborInDirection[cvf::StructGridInterface::POS_K] + isPossibleNeighborInDirection[cvf::StructGridInterface::NEG_K];
|
|
|
|
// If cell 2 is not adjacent with respect to any of the six ijk directions,
|
|
// assume that we have no overlapping area.
|
|
|
|
if ( !hasNeighbourInAnyDirection )
|
|
{
|
|
// Add to search map
|
|
// m_cellIdxToFaceToConnectionIdxMap[m_connections[cnIdx].m_c1GlobIdx][cvf::StructGridInterface::NO_FACE].push_back(cnIdx);
|
|
// m_cellIdxToFaceToConnectionIdxMap[m_connections[cnIdx].m_c2GlobIdx][cvf::StructGridInterface::NO_FACE].push_back(cnIdx);
|
|
|
|
// cvf::Trace::show("NNC: No direct neighbors : C1: " + cvf::String((int)m_connections[cnIdx].m_c1GlobIdx) +
|
|
// " C2: " + cvf::String((int)m_connections[cnIdx].m_c2GlobIdx));
|
|
return cvf::StructGridInterface::NO_FACE;
|
|
}
|
|
}
|
|
|
|
for ( unsigned char fIdx = 0; fIdx < 6; ++fIdx )
|
|
{
|
|
if ( !isPossibleNeighborInDirection[fIdx] )
|
|
{
|
|
continue;
|
|
}
|
|
|
|
// Calculate connection polygon
|
|
|
|
std::vector<size_t> polygon;
|
|
std::vector<cvf::Vec3d> intersections;
|
|
std::array<size_t, 4> face1;
|
|
std::array<size_t, 4> face2;
|
|
c1.faceIndices( ( cvf::StructGridInterface::FaceType )( fIdx ), &face1 );
|
|
c2.faceIndices( cvf::StructGridInterface::oppositeFace( ( cvf::StructGridInterface::FaceType )( fIdx ) ), &face2 );
|
|
|
|
bool foundOverlap = cvf::GeometryTools::calculateOverlapPolygonOfTwoQuads( polygon,
|
|
intersections,
|
|
(cvf::EdgeIntersectStorage<size_t>*)nullptr,
|
|
cvf::wrapArrayConst( &mainGrid.nodes() ),
|
|
face1.data(),
|
|
face2.data(),
|
|
1e-6 );
|
|
|
|
if ( foundOverlap )
|
|
{
|
|
if ( connectionPolygon ) ( *connectionPolygon ) = polygon;
|
|
if ( connectionIntersections ) ( *connectionIntersections ) = intersections;
|
|
return ( cvf::StructGridInterface::FaceType )( fIdx );
|
|
}
|
|
}
|
|
|
|
return cvf::StructGridInterface::NO_FACE;
|
|
}
|
|
|
|
void assignThreadConnections( RigConnectionContainer& allConnections, RigConnectionContainer& threadConnections )
|
|
{
|
|
#pragma omp critical( critical_section_RigCellFaceGeometryTools_assignThreadConnections )
|
|
{
|
|
allConnections.push_back( threadConnections );
|
|
}
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
RigConnectionContainer RigCellFaceGeometryTools::computeOtherNncs( const RigMainGrid* mainGrid,
|
|
const RigConnectionContainer& nativeConnections,
|
|
const RigActiveCellInfo* activeCellInfo,
|
|
bool includeInactiveCells )
|
|
{
|
|
// 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::set<std::pair<unsigned, unsigned>> nativeCellPairs;
|
|
|
|
for ( size_t i = 0; i < nativeConnections.size(); ++i )
|
|
{
|
|
RigConnection c = nativeConnections[i];
|
|
nativeCellPairs.emplace( static_cast<unsigned>( c.c1GlobIdx() ), static_cast<unsigned>( c.c2GlobIdx() ) );
|
|
}
|
|
|
|
if ( nativeConnections.size() != nativeCellPairs.size() )
|
|
{
|
|
QString message = QString( "Nnc connection imported from Eclipse are not unique\nNNC count : %1\nUnique : %2" )
|
|
.arg( nativeConnections.size() )
|
|
.arg( nativeCellPairs.size() );
|
|
RiaLogging::warning( message.toStdString(), "RigCellFaceGeometryTools" );
|
|
}
|
|
|
|
const cvf::Collection<RigFault>& faults = mainGrid->faults();
|
|
|
|
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();
|
|
|
|
// Build a vector of active face indices so we don't have to have this check inside the parallel loop.
|
|
// This makes the load balancing much better in the loop.
|
|
std::vector<size_t> activeFaceIndices;
|
|
activeFaceIndices.reserve( faultFaces.size() );
|
|
for ( size_t faceIdx = 0; faceIdx < faultFaces.size(); ++faceIdx )
|
|
{
|
|
const RigFault::FaultFace& f = faultFaces[faceIdx];
|
|
|
|
bool atLeastOneCellActive = true;
|
|
if ( !includeInactiveCells && activeCellInfo && activeCellInfo->reservoirActiveCellCount() > 0u )
|
|
{
|
|
atLeastOneCellActive = activeCellInfo->isActive( ReservoirCellIndex( f.m_nativeReservoirCellIndex ) ) ||
|
|
activeCellInfo->isActive( ReservoirCellIndex( f.m_oppositeReservoirCellIndex ) );
|
|
}
|
|
|
|
if ( atLeastOneCellActive ) activeFaceIndices.push_back( faceIdx );
|
|
}
|
|
|
|
#pragma omp parallel
|
|
{
|
|
RigConnectionContainer threadConnections;
|
|
#pragma omp for schedule( guided )
|
|
for ( int activeFaceIdx = 0; activeFaceIdx < static_cast<int>( activeFaceIndices.size() ); activeFaceIdx++ )
|
|
{
|
|
size_t faceIdx = activeFaceIndices[activeFaceIdx];
|
|
extractConnectionsForFace( faultFaces[faceIdx], mainGrid, nativeCellPairs, threadConnections );
|
|
}
|
|
|
|
// Merge together connections per thread
|
|
assignThreadConnections( otherConnections, threadConnections );
|
|
} // end parallel region
|
|
}
|
|
|
|
size_t nncCountWarningThreshold = 5000000;
|
|
if ( otherConnections.size() > nncCountWarningThreshold )
|
|
{
|
|
auto txt = QString( "Additional NNC count has reached %1, and is above the warning threshold of %2. Faults for "
|
|
"inactive cells can be managed from Preferences->Eclipse Grid->Include Inactive Cells" )
|
|
.arg( otherConnections.size() )
|
|
.arg( nncCountWarningThreshold );
|
|
RiaLogging::warning( txt.toStdString(), "RigCellFaceGeometryTools" );
|
|
}
|
|
|
|
otherConnections.remove_duplicates();
|
|
return otherConnections;
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
void RigCellFaceGeometryTools::extractConnectionsForFace( const RigFault::FaultFace& face,
|
|
const RigMainGrid* mainGrid,
|
|
const std::set<std::pair<unsigned, unsigned>>& nativeCellPairs,
|
|
RigConnectionContainer& connections )
|
|
{
|
|
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->cell( 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 );
|
|
|
|
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->ijkFromCellIndexUnguarded( sourceReservoirCellIndex, &i, &j, &k );
|
|
|
|
RigMainGrid::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 >= mainGrid->cellCount() )
|
|
{
|
|
continue;
|
|
}
|
|
|
|
if ( candidateCellIndex == neighborCellIndex )
|
|
{
|
|
// Exclude direct neighbor
|
|
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->ijkFromCellIndexUnguarded( candidateCellIndex, &ci, &cj, &ck );
|
|
|
|
// For I and J faces, only consider cells in the neighbor pillar column. Cells in other columns mostly
|
|
// touch the face along an edge, and testing them for overlap dominates the computation time.
|
|
auto gridAxis = cvf::StructGridInterface::gridAxisFromFace( sourceCellFace );
|
|
if ( gridAxis == cvf::StructGridInterface::GridAxisType::AXIS_I || gridAxis == cvf::StructGridInterface::GridAxisType::AXIS_J )
|
|
{
|
|
if ( ni != ci || nj != cj )
|
|
{
|
|
continue;
|
|
}
|
|
}
|
|
else if ( gridAxis == cvf::StructGridInterface::GridAxisType::AXIS_K )
|
|
{
|
|
if ( nk != ck )
|
|
{
|
|
continue;
|
|
}
|
|
}
|
|
}
|
|
|
|
// Test if this pair of cells already has a connection. Check both combinations of cell index ordering to avoid
|
|
// duplicate NNC geometry for the same pair of cells
|
|
|
|
auto candidate = std::make_pair( static_cast<unsigned>( sourceReservoirCellIndex ), static_cast<unsigned>( candidateCellIndex ) );
|
|
|
|
if ( nativeCellPairs.count( candidate ) > 0 )
|
|
{
|
|
continue;
|
|
}
|
|
|
|
std::swap( candidate.first, candidate.second );
|
|
|
|
if ( nativeCellPairs.count( candidate ) > 0 )
|
|
{
|
|
continue;
|
|
}
|
|
|
|
std::vector<size_t> polygon;
|
|
std::vector<cvf::Vec3d> intersections;
|
|
|
|
std::array<size_t, 4> candidateFaceIndices;
|
|
mainGrid->cell( 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( sourceReservoirCellIndex,
|
|
candidateCellIndex,
|
|
sourceCellFace,
|
|
RigCellFaceGeometryTools::extractPolygon( mainGridNodes, polygon, intersections ) );
|
|
|
|
connections.push_back( conn );
|
|
}
|
|
}
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
std::vector<cvf::Vec3f> RigCellFaceGeometryTools::extractPolygon( const std::vector<cvf::Vec3d>& nativeNodes,
|
|
const std::vector<size_t>& connectionPolygon,
|
|
const std::vector<cvf::Vec3d>& connectionIntersections )
|
|
{
|
|
std::vector<cvf::Vec3f> allPolygonNodes;
|
|
|
|
for ( size_t polygonIndex : connectionPolygon )
|
|
{
|
|
if ( polygonIndex < nativeNodes.size() )
|
|
allPolygonNodes.push_back( cvf::Vec3f( nativeNodes[polygonIndex] ) );
|
|
else
|
|
allPolygonNodes.push_back( cvf::Vec3f( connectionIntersections[polygonIndex - nativeNodes.size()] ) );
|
|
}
|
|
|
|
return allPolygonNodes;
|
|
}
|