diff --git a/ApplicationLibCode/ReservoirDataModel/Well/RigEclipseWellLogExtractor.cpp b/ApplicationLibCode/ReservoirDataModel/Well/RigEclipseWellLogExtractor.cpp index 851388a7b5..d278970963 100644 --- a/ApplicationLibCode/ReservoirDataModel/Well/RigEclipseWellLogExtractor.cpp +++ b/ApplicationLibCode/ReservoirDataModel/Well/RigEclipseWellLogExtractor.cpp @@ -28,6 +28,7 @@ #include "RigResultAccessor.h" #include "RigWellLogExtractionTools.h" #include "RigWellPath.h" +#include "RigWellPathGeometryTools.h" #include "RigWellPathIntersectionTools.h" #include "cvfBoundingBox.h" @@ -35,8 +36,12 @@ #include "cafAssert.h" +#include #include +#include +#include #include +#include //================================================================================================== /// @@ -107,6 +112,11 @@ void RigEclipseWellLogExtractor::calculateIntersection() insertIntersectionsInMap( intersections, p1, md1, p2, md2, tolerance, &uniqueIntersections ); } + if ( uniqueIntersections.empty() && m_caseData->mainGrid()->isRadial() ) + { + uniqueIntersections = radialGridAxisIntersections( tolerance ); + } + if ( uniqueIntersections.empty() && m_wellPathGeometry->wellPathPoints().size() > 1 ) { // When entering this function, all well path points are either completely outside the grid @@ -171,6 +181,86 @@ void RigEclipseWellLogExtractor::calculateIntersection() populateReturnArrays( uniqueIntersections ); } +//-------------------------------------------------------------------------------------------------- +/// A well in a radial grid is located on the grid axis, inside the inner radius or on a face between angular cells, and is not +/// intersected by any cell. Use the innermost cell column and place the cells along the well path by depth. Requires the innermost +/// ring to consist of a single 360 degree cell for every layer, as picking one of several angular sectors would be arbitrary. +//-------------------------------------------------------------------------------------------------- +std::map RigEclipseWellLogExtractor::radialGridAxisIntersections( double tolerance ) +{ + std::map uniqueIntersections; + + const RigMainGrid* mainGrid = m_caseData->mainGrid(); + + // Only the innermost ring can safely represent the well's position when its cells span the full 360 degrees for every layer. + // If the ring is split into several angular sectors, picking a single sector would be an arbitrary and potentially wrong choice. + if ( mainGrid->cellCountJ() != 1 ) return uniqueIntersections; + + const std::vector& wellPathMds = m_wellPathGeometry->measuredDepths(); + const std::vector wellPathTvds = m_wellPathGeometry->trueVerticalDepths(); + if ( wellPathMds.size() < 2 ) return uniqueIntersections; + + const auto [minTvdIt, maxTvdIt] = std::minmax_element( wellPathTvds.begin(), wellPathTvds.end() ); + const double minTvd = *minTvdIt; + const double maxTvd = *maxTvdIt; + + std::vector> cellDepths; + for ( size_t k = 0; k < mainGrid->cellCountK(); k++ ) + { + const size_t cellIndex = mainGrid->cellIndexFromIJK( 0, 0, k ); + if ( mainGrid->cell( cellIndex ).isInvalid() ) continue; + + double topTvd = std::numeric_limits::infinity(); + double bottomTvd = -std::numeric_limits::infinity(); + for ( const auto& corner : mainGrid->cellCornerVertices( cellIndex ) ) + { + topTvd = std::min( topTvd, -corner.z() ); + bottomTvd = std::max( bottomTvd, -corner.z() ); + } + + if ( bottomTvd <= minTvd || topTvd >= maxTvd ) continue; + + cellDepths.emplace_back( std::max( topTvd, minTvd ), std::min( bottomTvd, maxTvd ), cellIndex ); + } + + if ( cellDepths.empty() ) return uniqueIntersections; + + // MD interpolation requires unique TVD values ordered along the well path. K is not guaranteed to increase with depth, and + // neighbour cells share their boundary TVD. + std::vector boundaryTvds; + for ( const auto& [topTvd, bottomTvd, cellIndex] : cellDepths ) + { + boundaryTvds.push_back( topTvd ); + boundaryTvds.push_back( bottomTvd ); + } + std::sort( boundaryTvds.begin(), boundaryTvds.end() ); + boundaryTvds.erase( std::unique( boundaryTvds.begin(), boundaryTvds.end() ), boundaryTvds.end() ); + + const std::vector boundaryMds = RigWellPathGeometryTools::interpolateMdFromTvd( wellPathMds, wellPathTvds, boundaryTvds ); + if ( boundaryMds.size() != boundaryTvds.size() ) return uniqueIntersections; + + auto mdAtTvd = [&]( double tvd ) + { + const auto it = std::lower_bound( boundaryTvds.begin(), boundaryTvds.end(), tvd ); + return boundaryMds[std::distance( boundaryTvds.begin(), it )]; + }; + + for ( const auto& [topTvd, bottomTvd, cellIndex] : cellDepths ) + { + for ( bool isEntering : { true, false } ) + { + const double md = mdAtTvd( isEntering ? topTvd : bottomTvd ); + const cvf::Vec3d point = m_wellPathGeometry->interpolatedPointAlongWellPath( md ); + + HexIntersectionInfo info( point, isEntering, cvf::StructGridInterface::NO_FACE, cellIndex ); + RigMDCellIdxEnterLeaveKey enterLeaveKey( md, cellIndex, isEntering, tolerance ); + uniqueIntersections.insert( std::make_pair( enterLeaveKey, info ) ); + } + } + + return uniqueIntersections; +} + //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- diff --git a/ApplicationLibCode/ReservoirDataModel/Well/RigEclipseWellLogExtractor.h b/ApplicationLibCode/ReservoirDataModel/Well/RigEclipseWellLogExtractor.h index ac43441027..9b31f228d8 100644 --- a/ApplicationLibCode/ReservoirDataModel/Well/RigEclipseWellLogExtractor.h +++ b/ApplicationLibCode/ReservoirDataModel/Well/RigEclipseWellLogExtractor.h @@ -42,9 +42,10 @@ public: const RigEclipseCaseData* caseData() { return m_caseData.p(); } private: - void calculateIntersection(); - std::vector findCloseCellIndices( const cvf::BoundingBox& bb ); - cvf::Vec3d calculateLengthInCell( size_t cellIndex, const cvf::Vec3d& startPoint, const cvf::Vec3d& endPoint ) const override; + void calculateIntersection(); + std::map radialGridAxisIntersections( double tolerance ); + std::vector findCloseCellIndices( const cvf::BoundingBox& bb ); + cvf::Vec3d calculateLengthInCell( size_t cellIndex, const cvf::Vec3d& startPoint, const cvf::Vec3d& endPoint ) const override; double computeLengthThreshold() const; diff --git a/ApplicationLibCode/UnitTests/RigWellLogExtractor-Test.cpp b/ApplicationLibCode/UnitTests/RigWellLogExtractor-Test.cpp index 569319b525..61349dae3d 100644 --- a/ApplicationLibCode/UnitTests/RigWellLogExtractor-Test.cpp +++ b/ApplicationLibCode/UnitTests/RigWellLogExtractor-Test.cpp @@ -141,6 +141,122 @@ TEST( RigEclipseWellLogExtractor, HorizontalWellAcrossFlippedGrid ) } } +namespace +{ +//-------------------------------------------------------------------------------------------------- +/// 3x1x4 grid with TVD 1000-1004, located at x = [1, 4]. A vertical well at x = 0 is outside all cells, as a well on the axis of a +/// radial grid is. +//-------------------------------------------------------------------------------------------------- +cvf::ref createRadialAxisTestReservoir( bool isRadial, bool kIncreasesUpwards ) +{ + cvf::ref reservoir = new RigEclipseCaseData( nullptr ); + cvf::ref mockFileInterface = new RifReaderMockModel; + + const double topZ = -1000.0; + const double bottomZ = -1004.0; + if ( kIncreasesUpwards ) + { + mockFileInterface->setWorldCoordinates( cvf::Vec3d( 1, 0, bottomZ ), cvf::Vec3d( 4, 1, topZ ) ); + } + else + { + mockFileInterface->setWorldCoordinates( cvf::Vec3d( 1, 0, topZ ), cvf::Vec3d( 4, 1, bottomZ ) ); + } + mockFileInterface->setCellCounts( cvf::Vec3st( 3, 1, 4 ) ); + mockFileInterface->enableWellData( false ); + mockFileInterface->open( "", reservoir.p() ); + + reservoir->mainGrid()->setIsRadial( isRadial ); + reservoir->mainGrid()->computeCachedData(); + + return reservoir; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +std::vector radialAxisIntersections( RigEclipseCaseData* reservoir, double topTvd, double bottomTvd, double topMd ) +{ + cvf::ref wellPathGeometry = new RigWellPath; + wellPathGeometry->setWellPathPoints( { cvf::Vec3d( 0, 0.5, -topTvd ), cvf::Vec3d( 0, 0.5, -bottomTvd ) }, + { topMd, topMd + bottomTvd - topTvd } ); + + cvf::ref e = new RigEclipseWellLogExtractor( reservoir, wellPathGeometry.p(), "" ); + return e->cellIntersectionInfosAlongWellPath(); +} +} // namespace + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +TEST( RigEclipseWellLogExtractor, RadialGridWellOnAxisUsesInnermostColumn ) +{ + for ( bool kIncreasesUpwards : { false, true } ) + { + auto reservoir = createRadialAxisTestReservoir( true, kIncreasesUpwards ); + + auto intersections = radialAxisIntersections( reservoir.p(), 1000.0, 1004.0, 100.0 ); + ASSERT_EQ( 4u, intersections.size() ) << "kIncreasesUpwards=" << kIncreasesUpwards; + + for ( size_t i = 0; i < intersections.size(); i++ ) + { + const size_t k = kIncreasesUpwards ? 3 - i : i; + EXPECT_EQ( reservoir->mainGrid()->cellIndexFromIJK( 0, 0, k ), intersections[i].globCellIndex ); + EXPECT_NEAR( 100.0 + i, intersections[i].startMD, 1e-6 ); + EXPECT_NEAR( 101.0 + i, intersections[i].endMD, 1e-6 ); + } + } +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +TEST( RigEclipseWellLogExtractor, RadialGridWellOnAxisClipsCellsToWellPathDepth ) +{ + auto reservoir = createRadialAxisTestReservoir( true, false ); + + auto intersections = radialAxisIntersections( reservoir.p(), 1001.5, 1003.0, 0.0 ); + ASSERT_EQ( 2u, intersections.size() ); + + EXPECT_EQ( reservoir->mainGrid()->cellIndexFromIJK( 0, 0, 1 ), intersections[0].globCellIndex ); + EXPECT_NEAR( 0.0, intersections[0].startMD, 1e-6 ); + EXPECT_NEAR( 0.5, intersections[0].endMD, 1e-6 ); + + EXPECT_EQ( reservoir->mainGrid()->cellIndexFromIJK( 0, 0, 2 ), intersections[1].globCellIndex ); + EXPECT_NEAR( 0.5, intersections[1].startMD, 1e-6 ); + EXPECT_NEAR( 1.5, intersections[1].endMD, 1e-6 ); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +TEST( RigEclipseWellLogExtractor, NonRadialGridWellOutsideGridHasNoIntersections ) +{ + auto reservoir = createRadialAxisTestReservoir( false, false ); + + EXPECT_TRUE( radialAxisIntersections( reservoir.p(), 1000.0, 1004.0, 100.0 ).empty() ); +} + +//-------------------------------------------------------------------------------------------------- +/// The innermost ring fallback only applies when it is a single 360 degree cell for every layer. When the ring is split into +/// several angular sectors, picking a single sector would be arbitrary, so no fallback intersections should be inserted. +//-------------------------------------------------------------------------------------------------- +TEST( RigEclipseWellLogExtractor, RadialGridWellOnAxisWithMultipleAngularSectorsHasNoFallback ) +{ + cvf::ref reservoir = new RigEclipseCaseData( nullptr ); + cvf::ref mockFileInterface = new RifReaderMockModel; + + mockFileInterface->setWorldCoordinates( cvf::Vec3d( 1, 0, -1004.0 ), cvf::Vec3d( 4, 1, -1000.0 ) ); + mockFileInterface->setCellCounts( cvf::Vec3st( 3, 2, 4 ) ); + mockFileInterface->enableWellData( false ); + mockFileInterface->open( "", reservoir.p() ); + + reservoir->mainGrid()->setIsRadial( true ); + reservoir->mainGrid()->computeCachedData(); + + EXPECT_TRUE( radialAxisIntersections( reservoir.p(), 1000.0, 1004.0, 100.0 ).empty() ); +} + //-------------------------------------------------------------------------------------------------- /// setFlipAxis must actually mirror the grid node coordinates for every combination of Flip X / /// Flip Y. Covers the data-level part of https://github.com/OPM/ResInsight/issues/13967 where a