diff --git a/ApplicationLibCode/ReservoirDataModel/RigMainGrid.cpp b/ApplicationLibCode/ReservoirDataModel/RigMainGrid.cpp index 6b19a14d08..bafafa61e6 100644 --- a/ApplicationLibCode/ReservoirDataModel/RigMainGrid.cpp +++ b/ApplicationLibCode/ReservoirDataModel/RigMainGrid.cpp @@ -546,18 +546,39 @@ void RigMainGrid::calculateFaults( const RigActiveCellInfo* activeCellInfo ) const std::vector& vxs = m_mainGrid->nodes(); + // Each cell only writes the fault index of its own positive faces and the opposite faces of its neighbors, so the + // cells can be processed in parallel. Static scheduling assigns increasing cell ranges to increasing thread + // indices, so appending the thread results in thread order keeps the faces ordered by cell index. + const int threadCount = RiaOpenMPTools::availableThreadCount(); + std::vector> threadFaultFaces( threadCount ); + std::vector> threadFaultFacesInactive( threadCount ); + +#pragma omp parallel + { + const int myThread = RiaOpenMPTools::currentThreadIndex(); + +#pragma omp for schedule( static ) + for ( int gcIdx = 0; gcIdx < static_cast( totalCellCount() ); ++gcIdx ) + { + addUnNamedFaultFaces( gcIdx, + activeCellInfo, + vxs, + unNamedFaultIdx, + unNamedFaultWithInactiveIdx, + threadFaultFaces[myThread], + threadFaultFacesInactive[myThread], + m_faultsPrCellAcc.p() ); + } + } + std::vector& unNamedFaultFaces = unNamedFault->faultFaces(); std::vector& unNamedFaultFacesInactive = unNamedFaultWithInactive->faultFaces(); - for ( int gcIdx = 0; gcIdx < static_cast( totalCellCount() ); ++gcIdx ) + for ( int i = 0; i < threadCount; i++ ) { - addUnNamedFaultFaces( gcIdx, - activeCellInfo, - vxs, - unNamedFaultIdx, - unNamedFaultWithInactiveIdx, - unNamedFaultFaces, - unNamedFaultFacesInactive, - m_faultsPrCellAcc.p() ); + unNamedFaultFaces.insert( unNamedFaultFaces.end(), threadFaultFaces[i].begin(), threadFaultFaces[i].end() ); + unNamedFaultFacesInactive.insert( unNamedFaultFacesInactive.end(), + threadFaultFacesInactive[i].begin(), + threadFaultFacesInactive[i].end() ); } } @@ -596,13 +617,10 @@ void RigMainGrid::addUnNamedFaultFaces( int gcIdx, bool firstNO_FAULTFaceForCell = true; bool isCellActive = true; - char upperLimitForFaceType = cvf::StructGridInterface::FaceType::POS_K; - - // Compare only I and J faces - for ( char faceIdx = 0; faceIdx < upperLimitForFaceType; ++faceIdx ) + // Compare only I and J faces. A negative face is the positive face of the neighbor cell with a lower index, so it is + // handled when that neighbor is processed. + for ( auto face : { cvf::StructGridInterface::FaceType::POS_I, cvf::StructGridInterface::FaceType::POS_J } ) { - cvf::StructGridInterface::FaceType face = cvf::StructGridInterface::FaceType( faceIdx ); - // For faces that has no used defined Fault assigned: if ( m_faultsPrCellAcc->faultIdx( gcIdx, face ) == RigFaultsPrCellAccumulator::NO_FAULT ) @@ -661,7 +679,7 @@ void RigMainGrid::addUnNamedFaultFaces( int gcIdx, if ( static_cast( gcIdx ) < neighborReservoirCellIdx ) { - RigFault::FaultFace ff( gcIdx, cvf::StructGridInterface::FaceType( faceIdx ), neighborReservoirCellIdx ); + RigFault::FaultFace ff( gcIdx, face, neighborReservoirCellIdx ); if ( isCellActive && isNeighborCellActive ) { unNamedFaultFaces.push_back( ff ); diff --git a/ApplicationLibCode/UnitTests/CMakeLists.txt b/ApplicationLibCode/UnitTests/CMakeLists.txt index 8da8372a89..f9139b0c5d 100644 --- a/ApplicationLibCode/UnitTests/CMakeLists.txt +++ b/ApplicationLibCode/UnitTests/CMakeLists.txt @@ -31,6 +31,7 @@ set(SOURCE_UNITTEST_FILES ${CMAKE_CURRENT_LIST_DIR}/RigSgasResultCalculator-Test.cpp ${CMAKE_CURRENT_LIST_DIR}/RigPorvSoilSgasResultCalculator-Test.cpp ${CMAKE_CURRENT_LIST_DIR}/RigReservoir-Test.cpp + ${CMAKE_CURRENT_LIST_DIR}/RigMainGridFaults-Test.cpp ${CMAKE_CURRENT_LIST_DIR}/RigMainGridSearchTree-Test.cpp ${CMAKE_CURRENT_LIST_DIR}/RigResdataGridConverter-Test.cpp ${CMAKE_CURRENT_LIST_DIR}/RigGridExportAdapter-Test.cpp diff --git a/ApplicationLibCode/UnitTests/RigMainGridFaults-Test.cpp b/ApplicationLibCode/UnitTests/RigMainGridFaults-Test.cpp new file mode 100644 index 0000000000..4fbc36319b --- /dev/null +++ b/ApplicationLibCode/UnitTests/RigMainGridFaults-Test.cpp @@ -0,0 +1,146 @@ +///////////////////////////////////////////////////////////////////////////////// +// +// Copyright (C) 2026- 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 +// for more details. +// +///////////////////////////////////////////////////////////////////////////////// + +#include "gtest/gtest.h" + +#include "RiaResultNames.h" + +#include "RigActiveCellInfo.h" +#include "RigEclipseCaseData.h" +#include "RigFault.h" +#include "RigMainGrid.h" +#include "RigReservoirBuilder.h" + +namespace +{ +//-------------------------------------------------------------------------------------------------- +/// Build a 1 x 1 x cellCount row of cells along I (no file, no view), and shift every other cell +/// down in Z so that every I interface becomes a geometric (unnamed) fault. +//-------------------------------------------------------------------------------------------------- +cvf::ref buildStaggeredRowGridCase( int cellCount ) +{ + RigReservoirBuilder builder; + builder.setIJKCount( cvf::Vec3st( cellCount, 1, 1 ) ); + builder.setWorldCoordinates( cvf::Vec3d( 0.0, 0.0, 0.0 ), cvf::Vec3d( cellCount, 1.0, -1.0 ) ); + + cvf::ref eclipseCase = new RigEclipseCaseData( nullptr ); + builder.createGridsAndCells( eclipseCase.p() ); + eclipseCase->mainGrid()->computeCachedData(); + + std::vector& nodes = eclipseCase->mainGrid()->nodes(); + for ( int i = 1; i < cellCount; i += 2 ) + { + for ( int corner = 0; corner < 8; corner++ ) + { + nodes[static_cast( i ) * 8 + corner].z() -= 10.0; + } + } + + return eclipseCase; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +void findUnnamedFaults( RigMainGrid* mainGrid, const RigFault*& unnamedFault, const RigFault*& unnamedFaultInactive ) +{ + unnamedFault = nullptr; + unnamedFaultInactive = nullptr; + for ( const auto& fault : mainGrid->faults() ) + { + if ( fault->name() == RiaResultNames::undefinedGridFaultName() ) unnamedFault = fault.p(); + if ( fault->name() == RiaResultNames::undefinedGridFaultWithInactiveName() ) unnamedFaultInactive = fault.p(); + } +} + +} // namespace + +//-------------------------------------------------------------------------------------------------- +/// Every I interface in the staggered grid is a geometric fault. The unnamed fault faces are computed +/// in parallel, one cell range per thread, and the per-thread results are appended in thread order. +/// Verify that all interfaces are found, and that the resulting faces stay ordered by cell index. +//-------------------------------------------------------------------------------------------------- +TEST( RigMainGridFaults, UnnamedFaultFacesAreFoundInCellOrder ) +{ + const int cellCount = 10; + + cvf::ref eclipseCase = buildStaggeredRowGridCase( cellCount ); + RigMainGrid* mainGrid = eclipseCase->mainGrid(); + + RigActiveCellInfo* activeCellInfo = eclipseCase->activeCellInfo( RiaDefines::PorosityModelType::MATRIX_MODEL ); + mainGrid->calculateFaults( activeCellInfo ); + + const RigFault* unnamedFault = nullptr; + const RigFault* unnamedFaultInactive = nullptr; + findUnnamedFaults( mainGrid, unnamedFault, unnamedFaultInactive ); + + ASSERT_NE( unnamedFault, nullptr ); + ASSERT_NE( unnamedFaultInactive, nullptr ); + EXPECT_TRUE( unnamedFaultInactive->faultFaces().empty() ); + + const std::vector& faces = unnamedFault->faultFaces(); + ASSERT_EQ( static_cast( cellCount - 1 ), faces.size() ); + + for ( size_t idx = 0; idx < faces.size(); idx++ ) + { + EXPECT_EQ( idx, faces[idx].m_nativeReservoirCellIndex ) << "Face " << idx << " is out of cell-index order"; + EXPECT_EQ( idx + 1, faces[idx].m_oppositeReservoirCellIndex ); + EXPECT_EQ( cvf::StructGridInterface::POS_I, faces[idx].m_nativeFace ); + } +} + +//-------------------------------------------------------------------------------------------------- +/// Fault faces with an inactive native or neighbor cell are reported in a separate fault face list. +//-------------------------------------------------------------------------------------------------- +TEST( RigMainGridFaults, UnnamedFaultFacesWithInactiveCellGoToSeparateList ) +{ + const int cellCount = 4; + + cvf::ref eclipseCase = buildStaggeredRowGridCase( cellCount ); + RigMainGrid* mainGrid = eclipseCase->mainGrid(); + + // Build an independent active cell info with cell 2 left inactive (its entry is never assigned, so it stays + // at the default "undefined" value). The faces (1,2) and (2,3) should then be reported as inactive fault faces. + cvf::ref activeCellInfo = new RigActiveCellInfo; + activeCellInfo->setReservoirCellCount( static_cast( cellCount ) ); + size_t activeIdx = 0; + for ( int i = 0; i < cellCount; i++ ) + { + if ( i == 2 ) continue; + activeCellInfo->setCellResultIndex( ReservoirCellIndex( static_cast( i ) ), ActiveCellIndex( activeIdx++ ) ); + } + + mainGrid->calculateFaults( activeCellInfo.p() ); + + const RigFault* unnamedFault = nullptr; + const RigFault* unnamedFaultInactive = nullptr; + findUnnamedFaults( mainGrid, unnamedFault, unnamedFaultInactive ); + + ASSERT_NE( unnamedFault, nullptr ); + ASSERT_NE( unnamedFaultInactive, nullptr ); + + ASSERT_EQ( 1u, unnamedFault->faultFaces().size() ); + EXPECT_EQ( 0u, unnamedFault->faultFaces()[0].m_nativeReservoirCellIndex ); + EXPECT_EQ( 1u, unnamedFault->faultFaces()[0].m_oppositeReservoirCellIndex ); + + ASSERT_EQ( 2u, unnamedFaultInactive->faultFaces().size() ); + EXPECT_EQ( 1u, unnamedFaultInactive->faultFaces()[0].m_nativeReservoirCellIndex ); + EXPECT_EQ( 2u, unnamedFaultInactive->faultFaces()[0].m_oppositeReservoirCellIndex ); + EXPECT_EQ( 2u, unnamedFaultInactive->faultFaces()[1].m_nativeReservoirCellIndex ); + EXPECT_EQ( 3u, unnamedFaultInactive->faultFaces()[1].m_oppositeReservoirCellIndex ); +}