From eec9a88bd3b048c677824c091a5f58e35e0e6658 Mon Sep 17 00:00:00 2001 From: Magne Sjaastad Date: Mon, 24 Aug 2026 14:26:59 +0200 Subject: [PATCH] #14596 Collect cells per thread in polygon filter instead of using critical sections All four cell filter loops guarded a push_back into a shared container with an unnamed critical section, taking a process wide lock for every cell matching the polygon. On large grids with a permissive polygon this can be slower than running single threaded, and it made the resulting cell order vary between runs. Collect the cells in per thread buffers and append them in thread order after the parallel region, following the pattern used elsewhere in the code base. This removes eight critical sections. --- .../CellFilters/RimPolygonFilter.cpp | 235 +++++++++++------- 1 file changed, 145 insertions(+), 90 deletions(-) diff --git a/ApplicationLibCode/ProjectDataModel/CellFilters/RimPolygonFilter.cpp b/ApplicationLibCode/ProjectDataModel/CellFilters/RimPolygonFilter.cpp index a0b34c8193..fffc2eaceb 100644 --- a/ApplicationLibCode/ProjectDataModel/CellFilters/RimPolygonFilter.cpp +++ b/ApplicationLibCode/ProjectDataModel/CellFilters/RimPolygonFilter.cpp @@ -18,6 +18,8 @@ #include "RimPolygonFilter.h" +#include "RiaOpenMPTools.h" + #include "RigCellGeometryTools.h" #include "RigEclipseCaseData.h" #include "RigFemPartCollection.h" @@ -444,6 +446,22 @@ bool RimPolygonFilter::cellInsidePolygon2D( cvf::Vec3d center, std::array +void appendThreadCells( const std::vector>& threadCells, CellContainer& cells ) +{ + for ( const auto& cellsForThread : threadCells ) + { + cells.insert( cells.end(), cellsForThread.begin(), cellsForThread.end() ); + } +} +} // namespace + //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- @@ -451,29 +469,40 @@ void RimPolygonFilter::updateCellsDepthEclipse( const std::vector& p { // we should look in depth using Z coordinate const int gIdx = static_cast( grid->gridIndex() ); + + const int numberOfThreads = RiaOpenMPTools::availableThreadCount(); + std::vector> threadCells( numberOfThreads ); + // loop over all cells -#pragma omp parallel for - for ( int n = 0; n < (int)grid->cellCount(); n++ ) +#pragma omp parallel { - // valid cell? - RigCell cell = grid->cell( n ); - if ( cell.isInvalid() ) continue; + const int myThread = RiaOpenMPTools::currentThreadIndex(); - // get corner coordinates - std::array hexCorners = grid->cellCornerVertices( n ); - - // get cell ijk for k filter - size_t i, j, k; - grid->ijkFromCellIndex( n, &i, &j, &k ); - if ( !m_intervalTool.isNumberIncluded( k ) ) continue; - - // check if the polygon includes the cell - if ( cellInsidePolygon2D( cell.center(), hexCorners, points ) ) + // NB! We are inside a parallel section, do not use "parallel for" here +#pragma omp for + for ( int n = 0; n < (int)grid->cellCount(); n++ ) { -#pragma omp critical - m_cells[gIdx].push_back( n ); + // valid cell? + RigCell cell = grid->cell( n ); + if ( cell.isInvalid() ) continue; + + // get corner coordinates + std::array hexCorners = grid->cellCornerVertices( n ); + + // get cell ijk for k filter + size_t i, j, k; + grid->ijkFromCellIndex( n, &i, &j, &k ); + if ( !m_intervalTool.isNumberIncluded( k ) ) continue; + + // check if the polygon includes the cell + if ( cellInsidePolygon2D( cell.center(), hexCorners, points ) ) + { + threadCells[myThread].push_back( n ); + } } } + + appendThreadCells( threadCells, m_cells[gIdx] ); } //-------------------------------------------------------------------------------------------------- @@ -492,49 +521,57 @@ void RimPolygonFilter::updateCellsKIndexEclipse( const std::vector& const bool closedPolygon = isPolygonClosed(); const bool singlePoint = ( points.size() == 1 ); + const int numberOfThreads = RiaOpenMPTools::availableThreadCount(); + std::vector> threadCells( numberOfThreads ); + // find all cells in the K layer that matches the polygon -#pragma omp parallel for - for ( int i = 0; i < (int)grid->cellCountI(); i++ ) +#pragma omp parallel { - for ( size_t j = 0; j < grid->cellCountJ(); j++ ) + const int myThread = RiaOpenMPTools::currentThreadIndex(); + + // NB! We are inside a parallel section, do not use "parallel for" here +#pragma omp for + for ( int i = 0; i < (int)grid->cellCountI(); i++ ) { - size_t cellIdx = grid->cellIndexFromIJKUnguarded( i, j, K ); - const RigCell& cell = grid->cell( cellIdx ); - // valid cell? - if ( cell.isInvalid() ) continue; - - // get corner coordinates - std::array hexCorners = grid->cellCornerVertices( cellIdx ); - - if ( closedPolygon ) + for ( size_t j = 0; j < grid->cellCountJ(); j++ ) { - // check if the polygon includes the cell - if ( cellInsidePolygon2D( cell.center(), hexCorners, points ) ) + size_t cellIdx = grid->cellIndexFromIJKUnguarded( i, j, K ); + const RigCell& cell = grid->cell( cellIdx ); + // valid cell? + if ( cell.isInvalid() ) continue; + + // get corner coordinates + std::array hexCorners = grid->cellCornerVertices( cellIdx ); + + if ( closedPolygon ) { -#pragma omp critical - foundCells.push_back( cellIdx ); - } - } - else - { - if ( singlePoint ) - { - if ( RigCellGeometryTools::pointInsideCellNegK2D( points[0], hexCorners ) ) + // check if the polygon includes the cell + if ( cellInsidePolygon2D( cell.center(), hexCorners, points ) ) { -#pragma omp critical - foundCells.push_back( cellIdx ); + threadCells[myThread].push_back( cellIdx ); } } - // check if the polyline touches the top face of the cell - else if ( RigCellGeometryTools::polylineIntersectsCellNegK2D( points, hexCorners ) ) + else { -#pragma omp critical - foundCells.push_back( cellIdx ); + if ( singlePoint ) + { + if ( RigCellGeometryTools::pointInsideCellNegK2D( points[0], hexCorners ) ) + { + threadCells[myThread].push_back( cellIdx ); + } + } + // check if the polyline touches the top face of the cell + else if ( RigCellGeometryTools::polylineIntersectsCellNegK2D( points, hexCorners ) ) + { + threadCells[myThread].push_back( cellIdx ); + } } } } } + appendThreadCells( threadCells, foundCells ); + // now extend all these cells in one K layer to all K layers for ( const size_t cellIdx : foundCells ) { @@ -593,31 +630,41 @@ void RimPolygonFilter::updateCellsForEclipse( const std::vector& poi void RimPolygonFilter::updateCellsDepthGeoMech( const std::vector& points, const RigFemPartGrid* grid, int partId ) { // we should look in depth using Z coordinate + const int numberOfThreads = RiaOpenMPTools::availableThreadCount(); + std::vector> threadCells( numberOfThreads ); + // loop over all cells -#pragma omp parallel for - for ( int i = 0; i < (int)grid->cellCountI(); i++ ) +#pragma omp parallel { - for ( size_t j = 0; j < grid->cellCountJ(); j++ ) + const int myThread = RiaOpenMPTools::currentThreadIndex(); + + // NB! We are inside a parallel section, do not use "parallel for" here +#pragma omp for + for ( int i = 0; i < (int)grid->cellCountI(); i++ ) { - for ( size_t k = 0; k < grid->cellCountK(); k++ ) + for ( size_t j = 0; j < grid->cellCountJ(); j++ ) { - if ( !m_intervalTool.isNumberIncluded( k ) ) continue; - - size_t cellIdx = grid->cellIndexFromIJK( i, j, k ); - if ( cellIdx == cvf::UNDEFINED_SIZE_T ) continue; - - std::array vertices = grid->cellCornerVertices( cellIdx ); - cvf::Vec3d center = grid->cellCentroid( cellIdx ); - - // check if the polygon includes the cell - if ( cellInsidePolygon2D( center, vertices, points ) ) + for ( size_t k = 0; k < grid->cellCountK(); k++ ) { -#pragma omp critical - m_cells[partId].push_back( cellIdx ); + if ( !m_intervalTool.isNumberIncluded( k ) ) continue; + + size_t cellIdx = grid->cellIndexFromIJK( i, j, k ); + if ( cellIdx == cvf::UNDEFINED_SIZE_T ) continue; + + std::array vertices = grid->cellCornerVertices( cellIdx ); + cvf::Vec3d center = grid->cellCentroid( cellIdx ); + + // check if the polygon includes the cell + if ( cellInsidePolygon2D( center, vertices, points ) ) + { + threadCells[myThread].push_back( cellIdx ); + } } } } } + + appendThreadCells( threadCells, m_cells[partId] ); } //-------------------------------------------------------------------------------------------------- @@ -681,47 +728,55 @@ void RimPolygonFilter::updateCellsKIndexGeoMech( const std::vector& // find all cells in this K layer that matches the polygon std::list foundCells; -#pragma omp parallel for - for ( int i = 0; i < (int)grid->cellCountI(); i++ ) + const int numberOfThreads = RiaOpenMPTools::availableThreadCount(); + std::vector> threadCells( numberOfThreads ); + +#pragma omp parallel { - for ( size_t j = 0; j < grid->cellCountJ(); j++ ) + const int myThread = RiaOpenMPTools::currentThreadIndex(); + + // NB! We are inside a parallel section, do not use "parallel for" here +#pragma omp for + for ( int i = 0; i < (int)grid->cellCountI(); i++ ) { - size_t cellIdx = grid->cellIndexFromIJK( i, j, nk ); - if ( cellIdx == cvf::UNDEFINED_SIZE_T ) continue; - - // get corner coordinates - std::array hexCorners = grid->cellCornerVertices( cellIdx ); - cvf::Vec3d center = grid->cellCentroid( cellIdx ); - - if ( closedPolygon ) + for ( size_t j = 0; j < grid->cellCountJ(); j++ ) { - // check if the polygon includes the cell - if ( cellInsidePolygon2D( center, hexCorners, points ) ) + size_t cellIdx = grid->cellIndexFromIJK( i, j, nk ); + if ( cellIdx == cvf::UNDEFINED_SIZE_T ) continue; + + // get corner coordinates + std::array hexCorners = grid->cellCornerVertices( cellIdx ); + cvf::Vec3d center = grid->cellCentroid( cellIdx ); + + if ( closedPolygon ) { -#pragma omp critical - foundCells.push_back( cellIdx ); - } - } - else - { - if ( singlePoint ) - { - if ( RigCellGeometryTools::pointInsideCellNegK2D( points[0], hexCorners ) ) + // check if the polygon includes the cell + if ( cellInsidePolygon2D( center, hexCorners, points ) ) { -#pragma omp critical - foundCells.push_back( cellIdx ); + threadCells[myThread].push_back( cellIdx ); } } - // check if the polyline touches the top face of the cell - else if ( RigCellGeometryTools::polylineIntersectsCellNegK2D( points, hexCorners ) ) + else { -#pragma omp critical - foundCells.push_back( cellIdx ); + if ( singlePoint ) + { + if ( RigCellGeometryTools::pointInsideCellNegK2D( points[0], hexCorners ) ) + { + threadCells[myThread].push_back( cellIdx ); + } + } + // check if the polyline touches the top face of the cell + else if ( RigCellGeometryTools::polylineIntersectsCellNegK2D( points, hexCorners ) ) + { + threadCells[myThread].push_back( cellIdx ); + } } } } } + appendThreadCells( threadCells, foundCells ); + // now extend all these cells in one K layer to all K layers for ( const size_t cellIdx : foundCells ) {