Scale MAPAXES to grid units in the opm-common grid reader

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<double, 6> mapAxes;
+    for ( size_t i = 0; i < opmMapAxes.size(); ++i )
     {
-        std::array<double, 6> 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<long>( mainGrid->nodes().size() ); i++ )
-        {
-            auto& n = mainGrid->nodes()[i];
-            n.transformPoint( transform );
-        }
+    for ( long i = 0; i < static_cast<long>( mainGrid->nodes().size() ); i++ )
+    {
+        auto& n = mainGrid->nodes()[i];
+        n.transformPoint( transform );
     }
-
-    return true;
 }
 
 //--------------------------------------------------------------------------------------------------

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,

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<double, 6> 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<long>( activeGrid->nodes().size() ); i++ )
-        {
-            auto& n = activeGrid->nodes()[i];
-            n.transformPoint( transform );
-        }
-    }
+    applyMapAxes( opmGrid, activeGrid );
 
     return true;
 }

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

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 <http://www.gnu.org/licenses/gpl.html>
+//  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 ) );
+}

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<RigEclipseCaseData> buildBoxGrid( int ni, int nj, int nk )
+cvf::ref<RigEclipseCaseData> buildBoxGridForPartGenerator( int ni, int nj, int nk )
 {
     RigReservoirBuilder builder;
     builder.setIJKCount( cvf::Vec3st( ni, nj, nk ) );
@@ -51,7 +51,7 @@ cvf::ref<RigEclipseCaseData> buildBoxGrid( int ni, int nj, int nk )
 //--------------------------------------------------------------------------------------------------
 TEST( RivSingleCellPartGeneratorTest, ValidCellIndexCreatesDrawable )
 {
-    cvf::ref<RigEclipseCaseData> caseData = buildBoxGrid( 2, 3, 4 );
+    cvf::ref<RigEclipseCaseData> 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<RigEclipseCaseData> caseData = buildBoxGrid( 2, 3, 4 );
+    cvf::ref<RigEclipseCaseData> caseData = buildBoxGridForPartGenerator( 2, 3, 4 );
 
     const size_t cellCount = caseData->mainGrid()->cellCount();
 
@@ -84,7 +84,7 @@ TEST( RivSingleCellPartGeneratorTest, CellIndexOutOfBoundsIsIgnored )
 //--------------------------------------------------------------------------------------------------
 TEST( RivSingleCellPartGeneratorTest, GridIndexOutOfBoundsIsIgnored )
 {
-    cvf::ref<RigEclipseCaseData> caseData = buildBoxGrid( 2, 3, 4 );
+    cvf::ref<RigEclipseCaseData> 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<RigEclipseCaseData> caseData = buildBoxGrid( 2, 3, 4 );
+    cvf::ref<RigEclipseCaseData> caseData = buildBoxGridForPartGenerator( 2, 3, 4 );
 
     RivSingleCellPartGenerator partGen( caseData.p(), 0, cvf::UNDEFINED_SIZE_T, cvf::Vec3d::ZERO );