OilfieldToolsNet/ResInsight
4
0
Fork 0
mirror of https://github.com/OPM/ResInsight.git synced 2025-02-25 18:55:39 -06:00
Code Issues Packages Projects Releases Wiki Activity

Parallelize another loop in fault calculation

Browse Source
This commit is contained in:
Gaute Lindkvist 2020-05-14 11:36:20 +02:00
parent c5661f3fce
commit 3a526ab555
3 changed files with 78 additions and 56 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

11
ApplicationCode/ReservoirDataModel/RigGridBase.cpp
View File

@ -351,6 +351,17 @@ bool RigGridBase::cellIJKNeighbor( size_t i, size_t j, size_t k, FaceType face,
return true;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RigGridBase::cellIJKNeighborUnguarded( size_t i, size_t j, size_t k, FaceType face, size_t* neighborCellIndex ) const
{
size_t ni, nj, nk;
neighborIJKAtCellFace( i, j, k, face, &ni, &nj, &nk );
*neighborCellIndex = cellIndexFromIJK( ni, nj, nk );
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------

1
ApplicationCode/ReservoirDataModel/RigGridBase.h
View File

@ -106,6 +106,7 @@ public:
bool isCellValid( size_t i, size_t j, size_t k ) const override;
bool cellIJKNeighbor( size_t i, size_t j, size_t k, FaceType face, size_t* neighborCellIndex ) const override;
void cellIJKNeighborUnguarded( size_t i, size_t j, size_t k, FaceType face, size_t* neighborCellIndex ) const;
private:
std::string m_gridName;

122
ApplicationCode/ReservoirDataModel/RigMainGrid.cpp
View File

@ -450,16 +450,34 @@ void RigMainGrid::calculateFaults( const RigActiveCellInfo* activeCellInfo )
std::vector<RigFault::FaultFace>& unNamedFaultFaces = unNamedFault->faultFaces();
std::vector<RigFault::FaultFace>& unNamedFaultFacesInactive = unNamedFaultWithInactive->faultFaces();
for ( int gcIdx = 0; gcIdx < static_cast<int>( m_cells.size() ); ++gcIdx )
#pragma omp parallel
{
addUnNamedFaultFaces( gcIdx,
activeCellInfo,
vxs,
unNamedFaultIdx,
unNamedFaultWithInactiveIdx,
unNamedFaultFaces,
unNamedFaultFacesInactive,
m_faultsPrCellAcc.p() );
std::vector<RigFault::FaultFace> threadUnNamedFaultFaces;
std::vector<RigFault::FaultFace> threadUnNamedFaultFacesInactive;
#pragma omp for schedule( guided )
for ( int gcIdx = 0; gcIdx < static_cast<int>( m_cells.size() ); ++gcIdx )
{
addUnNamedFaultFaces( gcIdx,
activeCellInfo,
vxs,
unNamedFaultIdx,
unNamedFaultWithInactiveIdx,
threadUnNamedFaultFaces,
threadUnNamedFaultFacesInactive,
m_faultsPrCellAcc.p() );
}
#pragma omp critical
{
unNamedFaultFaces.insert( unNamedFaultFaces.end(),
threadUnNamedFaultFaces.begin(),
threadUnNamedFaultFaces.end() );
unNamedFaultFacesInactive.insert( unNamedFaultFacesInactive.end(),
threadUnNamedFaultFacesInactive.begin(),
threadUnNamedFaultFacesInactive.end() );
}
}
}
@ -480,6 +498,8 @@ void RigMainGrid::addUnNamedFaultFaces( int gcIdx,
return;
}
const double tolerance = 1e-6;
size_t neighborReservoirCellIdx;
size_t neighborGridCellIdx;
size_t i = 0;
@ -490,7 +510,7 @@ void RigMainGrid::addUnNamedFaultFaces( int gcIdx,
bool firstNO_FAULTFaceForCell = true;
bool isCellActive = true;
char upperLimitForFaceType = cvf::StructGridInterface::FaceType::POS_K;
const char upperLimitForFaceType = cvf::StructGridInterface::FaceType::POS_K;
// Compare only I and J faces
for ( char faceIdx = 0; faceIdx < upperLimitForFaceType; ++faceIdx )
@ -513,52 +533,51 @@ void RigMainGrid::addUnNamedFaultFaces( int gcIdx,
firstNO_FAULTFaceForCell = false;
}
if ( !hostGrid->cellIJKNeighbor( i, j, k, face, &neighborGridCellIdx ) )
{
continue;
}
hostGrid->cellIJKNeighborUnguarded( i, j, k, face, &neighborGridCellIdx );
neighborReservoirCellIdx = hostGrid->reservoirCellIndex( neighborGridCellIdx );
if ( m_cells[neighborReservoirCellIdx].isInvalid() )
if ( neighborReservoirCellIdx >= m_cells.size() || m_cells[neighborReservoirCellIdx].isInvalid() )
{
continue;
}
bool isNeighborCellActive = activeCellInfo->isActive( neighborReservoirCellIdx );
double tolerance = 1e-6;
std::array<size_t, 4> faceIdxs;
m_cells[gcIdx].faceIndices( face, &faceIdxs );
std::array<size_t, 4> nbFaceIdxs;
m_cells[neighborReservoirCellIdx].faceIndices( StructGridInterface::oppositeFace( face ), &nbFaceIdxs );
bool sharedFaceVertices = true;
if ( sharedFaceVertices && vxs[faceIdxs[0]].pointDistance( vxs[nbFaceIdxs[0]] ) > tolerance )
sharedFaceVertices = false;
if ( sharedFaceVertices && vxs[faceIdxs[1]].pointDistance( vxs[nbFaceIdxs[3]] ) > tolerance )
sharedFaceVertices = false;
if ( sharedFaceVertices && vxs[faceIdxs[2]].pointDistance( vxs[nbFaceIdxs[2]] ) > tolerance )
sharedFaceVertices = false;
if ( sharedFaceVertices && vxs[faceIdxs[3]].pointDistance( vxs[nbFaceIdxs[1]] ) > tolerance )
sharedFaceVertices = false;
if ( sharedFaceVertices )
{
continue;
}
// To avoid doing this calculation for the opposite face
int faultIdx = unNamedFaultIdx;
if ( !( isCellActive && isNeighborCellActive ) ) faultIdx = unNamedFaultWithInactiveIdx;
faultsPrCellAcc->setFaultIdx( gcIdx, face, faultIdx );
faultsPrCellAcc->setFaultIdx( neighborReservoirCellIdx, StructGridInterface::oppositeFace( face ), faultIdx );
// Add as fault face only if the grid index is less than the neighbors
if ( static_cast<size_t>( gcIdx ) < neighborReservoirCellIdx )
{
bool isNeighborCellActive = activeCellInfo->isActive( neighborReservoirCellIdx );
std::array<size_t, 4> faceIdxs;
m_cells[gcIdx].faceIndices( face, &faceIdxs );
std::array<size_t, 4> nbFaceIdxs;
m_cells[neighborReservoirCellIdx].faceIndices( StructGridInterface::oppositeFace( face ), &nbFaceIdxs );
bool sharedFaceVertices = true;
if ( sharedFaceVertices && vxs[faceIdxs[0]].pointDistance( vxs[nbFaceIdxs[0]] ) > tolerance )
sharedFaceVertices = false;
if ( sharedFaceVertices && vxs[faceIdxs[1]].pointDistance( vxs[nbFaceIdxs[3]] ) > tolerance )
sharedFaceVertices = false;
if ( sharedFaceVertices && vxs[faceIdxs[2]].pointDistance( vxs[nbFaceIdxs[2]] ) > tolerance )
sharedFaceVertices = false;
if ( sharedFaceVertices && vxs[faceIdxs[3]].pointDistance( vxs[nbFaceIdxs[1]] ) > tolerance )
sharedFaceVertices = false;
if ( sharedFaceVertices )
{
continue;
}
// To avoid doing this calculation for the opposite face
int faultIdx = unNamedFaultIdx;
if ( !( isCellActive && isNeighborCellActive ) ) faultIdx = unNamedFaultWithInactiveIdx;
faultsPrCellAcc->setFaultIdx( gcIdx, face, faultIdx );
#pragma omp critical
{
faultsPrCellAcc->setFaultIdx( neighborReservoirCellIdx,
StructGridInterface::oppositeFace( face ),
faultIdx );
}
RigFault::FaultFace ff( gcIdx, cvf::StructGridInterface::FaceType( faceIdx ), neighborReservoirCellIdx );
if ( isCellActive && isNeighborCellActive )
{
@ -569,15 +588,6 @@ void RigMainGrid::addUnNamedFaultFaces( int gcIdx,
unNamedFaultFacesInactive.push_back( ff );
}
}
else
{
CVF_FAIL_MSG( "Found fault with global neighbor index less than the native index. " ); // Should never
// occur. because
// we flag the
// opposite face
// in the
// faultsPrCellAcc
}
}
}
}