diff --git a/ApplicationLibCode/FileInterface/RifReaderOpmCommon.cpp b/ApplicationLibCode/FileInterface/RifReaderOpmCommon.cpp index a9175888b0..9e59cc60d0 100644 --- a/ApplicationLibCode/FileInterface/RifReaderOpmCommon.cpp +++ b/ApplicationLibCode/FileInterface/RifReaderOpmCommon.cpp @@ -19,6 +19,7 @@ #include "RifReaderOpmCommon.h" #include "RiaEclipseFileNameTools.h" +#include "RiaEclipseUnitTools.h" #include "RiaLogging.h" #include "RiaOpmParserTools.h" #include "RiaPreferencesSystem.h" @@ -341,35 +342,80 @@ bool RifReaderOpmCommon::importGrid( RigMainGrid* mainGrid, RigEclipseCaseData* transferStaticNNCData( opmGrid, lgrGrids, mainGrid ); } - auto opmMapAxes = opmGrid.get_mapaxes(); - if ( opmMapAxes.size() == 6 ) + applyMapAxes( opmGrid, mainGrid ); + + return true; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +double RifReaderOpmCommon::mapAxesScaleFactor( const std::string& mapUnits, int gridUnit ) +{ + // opm-common scales MAPAXES to meter based on the MAPUNITS keyword, see EGrid.cpp. The string comparisons must + // mirror opm-common exactly, as the intention is to detect what opm-common did and not what the file says. + double mapUnitInMeter = 0.0; + if ( mapUnits == "METRES" ) + mapUnitInMeter = 1.0; + else if ( mapUnits == "FEET" ) + mapUnitInMeter = RiaEclipseUnitTools::meterPerFeet(); + else if ( mapUnits == "CM" ) + mapUnitInMeter = 0.01; + + // MAPUNITS is missing or not recognized by opm-common, the map axes are left unscaled + if ( mapUnitInMeter == 0.0 ) return 1.0; + + // 1 = Metric, 2 = Field, 3 = Lab + double gridUnitInMeter = 0.0; + if ( gridUnit == 1 ) + gridUnitInMeter = 1.0; + else if ( gridUnit == 2 ) + gridUnitInMeter = RiaEclipseUnitTools::meterPerFeet(); + else if ( gridUnit == 3 ) + gridUnitInMeter = 0.01; + + // Unknown grid unit, undo the scaling applied by opm-common + if ( gridUnitInMeter == 0.0 ) return 1.0 / mapUnitInMeter; + + // Convert from meter to grid unit + return 1.0 / gridUnitInMeter; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +void RifReaderOpmCommon::applyMapAxes( Opm::EclIO::EGrid& opmGrid, RigMainGrid* mainGrid ) +{ + // get_mapaxes() returns a fixed size array that is only assigned when the MAPAXES keyword is present + if ( !opmGrid.with_mapaxes() ) return; + + // The cell corner coordinates are in grid units, scale the map axes to match before they are used as a translation + const double scaleFactor = mapAxesScaleFactor( opmGrid.get_mapunits(), m_gridUnit ); + + const auto& opmMapAxes = opmGrid.get_mapaxes(); + std::array mapAxes; + for ( size_t i = 0; i < opmMapAxes.size(); ++i ) { - std::array mapAxes; - for ( size_t i = 0; i < opmMapAxes.size(); ++i ) - { - mapAxes[i] = opmMapAxes[i]; - } + mapAxes[i] = opmMapAxes[i] * scaleFactor; + } - double norm_denominator = mapAxes[2] * mapAxes[5] - mapAxes[4] * mapAxes[3]; + double norm_denominator = mapAxes[2] * mapAxes[5] - mapAxes[4] * mapAxes[3]; - // Set the map axes transformation matrix on the main grid - mainGrid->setMapAxes( mapAxes ); - mainGrid->setUseMapAxes( norm_denominator != 0.0 ); + // Set the map axes transformation matrix on the main grid + mainGrid->setMapAxes( mapAxes ); + mainGrid->setUseMapAxes( norm_denominator != 0.0 ); - auto transform = mainGrid->mapAxisTransform(); + auto transform = mainGrid->mapAxisTransform(); - // Invert the transformation matrix to convert from file coordinates to domain coordinates - transform.invert(); + // Invert the transformation matrix to convert from file coordinates to domain coordinates + transform.invert(); #pragma omp parallel for - for ( long i = 0; i < static_cast( mainGrid->nodes().size() ); i++ ) - { - auto& n = mainGrid->nodes()[i]; - n.transformPoint( transform ); - } + for ( long i = 0; i < static_cast( mainGrid->nodes().size() ); i++ ) + { + auto& n = mainGrid->nodes()[i]; + n.transformPoint( transform ); } - - return true; } //-------------------------------------------------------------------------------------------------- diff --git a/ApplicationLibCode/FileInterface/RifReaderOpmCommon.h b/ApplicationLibCode/FileInterface/RifReaderOpmCommon.h index fa2690a12a..dc88fd961a 100644 --- a/ApplicationLibCode/FileInterface/RifReaderOpmCommon.h +++ b/ApplicationLibCode/FileInterface/RifReaderOpmCommon.h @@ -74,9 +74,15 @@ class RifReaderOpmCommon : public RifReaderInterface static GridDimensions readGridDimensions( const QString& gridFileName ); + // opm-common scales MAPAXES to meter based on the MAPUNITS keyword, while the cell corner coordinates are left in + // the units given by GRIDUNIT. Returns the factor required to bring the map axes into grid units. + static double mapAxesScaleFactor( const std::string& mapUnits, int gridUnit ); + protected: virtual bool importGrid( RigMainGrid* mainGrid, RigEclipseCaseData* caseData ); + void applyMapAxes( Opm::EclIO::EGrid& opmGrid, RigMainGrid* mainGrid ); + void transferActiveCells( Opm::EclIO::EGrid& opmGrid, size_t cellStartIndex, RigEclipseCaseData* eclipseCaseData, diff --git a/ApplicationLibCode/FileInterface/RifReaderOpmCommonActive.cpp b/ApplicationLibCode/FileInterface/RifReaderOpmCommonActive.cpp index 429180e676..33ba39a995 100644 --- a/ApplicationLibCode/FileInterface/RifReaderOpmCommonActive.cpp +++ b/ApplicationLibCode/FileInterface/RifReaderOpmCommonActive.cpp @@ -222,33 +222,7 @@ bool RifReaderOpmCommonActive::importGrid( RigMainGrid* /* mainGrid*/, RigEclips transferStaticNNCData( opmGrid, lgrGrids, activeGrid ); } - auto opmMapAxes = opmGrid.get_mapaxes(); - if ( opmMapAxes.size() == 6 ) - { - std::array mapAxes; - for ( size_t i = 0; i < opmMapAxes.size(); ++i ) - { - mapAxes[i] = opmMapAxes[i]; - } - - double norm_denominator = mapAxes[2] * mapAxes[5] - mapAxes[4] * mapAxes[3]; - - // Set the map axes transformation matrix on the main grid - activeGrid->setMapAxes( mapAxes ); - activeGrid->setUseMapAxes( norm_denominator != 0.0 ); - - auto transform = activeGrid->mapAxisTransform(); - - // Invert the transformation matrix to convert from file coordinates to domain coordinates - transform.invert(); - -#pragma omp parallel for - for ( long i = 0; i < static_cast( activeGrid->nodes().size() ); i++ ) - { - auto& n = activeGrid->nodes()[i]; - n.transformPoint( transform ); - } - } + applyMapAxes( opmGrid, activeGrid ); return true; } diff --git a/ApplicationLibCode/UnitTests/CMakeLists.txt b/ApplicationLibCode/UnitTests/CMakeLists.txt index 922bab36c9..d5b12cd68d 100644 --- a/ApplicationLibCode/UnitTests/CMakeLists.txt +++ b/ApplicationLibCode/UnitTests/CMakeLists.txt @@ -17,6 +17,7 @@ set(SOURCE_UNITTEST_FILES ${CMAKE_CURRENT_LIST_DIR}/RifEclipseOutputFileTools-Test.cpp ${CMAKE_CURRENT_LIST_DIR}/RifOpmFlowDeckFile-Test.cpp ${CMAKE_CURRENT_LIST_DIR}/RifReaderEclipseOutput-Test.cpp + ${CMAKE_CURRENT_LIST_DIR}/RifReaderOpmCommon-Test.cpp ${CMAKE_CURRENT_LIST_DIR}/RifReaderEclipseSummary-Test.cpp ${CMAKE_CURRENT_LIST_DIR}/RigActiveCellInfo-Test.cpp ${CMAKE_CURRENT_LIST_DIR}/RigEclipseCaseDataTools-Test.cpp diff --git a/ApplicationLibCode/UnitTests/RifReaderOpmCommon-Test.cpp b/ApplicationLibCode/UnitTests/RifReaderOpmCommon-Test.cpp new file mode 100644 index 0000000000..5d60353c74 --- /dev/null +++ b/ApplicationLibCode/UnitTests/RifReaderOpmCommon-Test.cpp @@ -0,0 +1,57 @@ +///////////////////////////////////////////////////////////////////////////////// +// +// 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 "RifReaderOpmCommon.h" + +#include "RiaEclipseUnitTools.h" + +//-------------------------------------------------------------------------------------------------- +/// opm-common scales MAPAXES to meter based on MAPUNITS, while the cell corner coordinates are left in the units +/// given by GRIDUNIT. Verify that the map axes are scaled back into grid units. +//-------------------------------------------------------------------------------------------------- +TEST( RifReaderOpmCommon, MapAxesScaleFactor ) +{ + const int gridUnitMetric = 1; + const int gridUnitField = 2; + const int gridUnitLab = 3; + const int gridUnitUnknown = -1; + + const double feetToMeter = RiaEclipseUnitTools::meterPerFeet(); + + // MAPUNITS is missing, opm-common has not scaled the map axes + EXPECT_DOUBLE_EQ( 1.0, RifReaderOpmCommon::mapAxesScaleFactor( "", gridUnitField ) ); + EXPECT_DOUBLE_EQ( 1.0, RifReaderOpmCommon::mapAxesScaleFactor( "", gridUnitMetric ) ); + + // MAPUNITS is not recognized by opm-common, the map axes are left unscaled + EXPECT_DOUBLE_EQ( 1.0, RifReaderOpmCommon::mapAxesScaleFactor( "FT", gridUnitField ) ); + + // Same unit for map axes and grid, the scaling applied by opm-common must be undone + EXPECT_DOUBLE_EQ( 1.0, RifReaderOpmCommon::mapAxesScaleFactor( "METRES", gridUnitMetric ) ); + EXPECT_DOUBLE_EQ( 1.0 / feetToMeter, RifReaderOpmCommon::mapAxesScaleFactor( "FEET", gridUnitField ) ); + EXPECT_DOUBLE_EQ( 100.0, RifReaderOpmCommon::mapAxesScaleFactor( "CM", gridUnitLab ) ); + + // Map axes and grid have different units, convert from meter to grid unit + EXPECT_DOUBLE_EQ( 1.0, RifReaderOpmCommon::mapAxesScaleFactor( "FEET", gridUnitMetric ) ); + EXPECT_DOUBLE_EQ( 1.0 / feetToMeter, RifReaderOpmCommon::mapAxesScaleFactor( "METRES", gridUnitField ) ); + + // Unknown grid unit, undo the scaling applied by opm-common + EXPECT_DOUBLE_EQ( 1.0 / feetToMeter, RifReaderOpmCommon::mapAxesScaleFactor( "FEET", gridUnitUnknown ) ); + EXPECT_DOUBLE_EQ( 1.0, RifReaderOpmCommon::mapAxesScaleFactor( "METRES", gridUnitUnknown ) ); +} diff --git a/ApplicationLibCode/UnitTests/RivSingleCellPartGenerator-Test.cpp b/ApplicationLibCode/UnitTests/RivSingleCellPartGenerator-Test.cpp index d783789d2f..cdff505ba9 100644 --- a/ApplicationLibCode/UnitTests/RivSingleCellPartGenerator-Test.cpp +++ b/ApplicationLibCode/UnitTests/RivSingleCellPartGenerator-Test.cpp @@ -32,7 +32,7 @@ namespace //-------------------------------------------------------------------------------------------------- /// Build a regular ni x nj x nk box grid in memory (no file, no view) //-------------------------------------------------------------------------------------------------- -cvf::ref buildBoxGrid( int ni, int nj, int nk ) +cvf::ref buildBoxGridForPartGenerator( int ni, int nj, int nk ) { RigReservoirBuilder builder; builder.setIJKCount( cvf::Vec3st( ni, nj, nk ) ); @@ -51,7 +51,7 @@ cvf::ref buildBoxGrid( int ni, int nj, int nk ) //-------------------------------------------------------------------------------------------------- TEST( RivSingleCellPartGeneratorTest, ValidCellIndexCreatesDrawable ) { - cvf::ref caseData = buildBoxGrid( 2, 3, 4 ); + cvf::ref caseData = buildBoxGridForPartGenerator( 2, 3, 4 ); RivSingleCellPartGenerator partGen( caseData.p(), 0, 5, cvf::Vec3d::ZERO ); @@ -65,7 +65,7 @@ TEST( RivSingleCellPartGeneratorTest, ValidCellIndexCreatesDrawable ) //-------------------------------------------------------------------------------------------------- TEST( RivSingleCellPartGeneratorTest, CellIndexOutOfBoundsIsIgnored ) { - cvf::ref caseData = buildBoxGrid( 2, 3, 4 ); + cvf::ref caseData = buildBoxGridForPartGenerator( 2, 3, 4 ); const size_t cellCount = caseData->mainGrid()->cellCount(); @@ -84,7 +84,7 @@ TEST( RivSingleCellPartGeneratorTest, CellIndexOutOfBoundsIsIgnored ) //-------------------------------------------------------------------------------------------------- TEST( RivSingleCellPartGeneratorTest, GridIndexOutOfBoundsIsIgnored ) { - cvf::ref caseData = buildBoxGrid( 2, 3, 4 ); + cvf::ref caseData = buildBoxGridForPartGenerator( 2, 3, 4 ); RivSingleCellPartGenerator partGen( caseData.p(), 1, 0, cvf::Vec3d::ZERO ); @@ -98,7 +98,7 @@ TEST( RivSingleCellPartGeneratorTest, GridIndexOutOfBoundsIsIgnored ) //-------------------------------------------------------------------------------------------------- TEST( RivSingleCellPartGeneratorTest, UndefinedCellIndexIsIgnored ) { - cvf::ref caseData = buildBoxGrid( 2, 3, 4 ); + cvf::ref caseData = buildBoxGridForPartGenerator( 2, 3, 4 ); RivSingleCellPartGenerator partGen( caseData.p(), 0, cvf::UNDEFINED_SIZE_T, cvf::Vec3d::ZERO );