API for Exposing Basic Information About Radial Grids

[bska]

This is an API design issue and I welcome any and all comments.

I am currently working on implementing a transmissibility calculation for radial grids (assuming essentially the OLDTRAN method and diagonal tensors in cylinder coordinates, no support for faults). The calculation depends on

1. Being able to identify whether or not a simulation case uses radial coordinates (i.e., essentially whether or not keyword RADIAL is present in RUNSPEC).

2. Extracting the (positive) inner radius of the model geometry (keyword INRAD from the GRID section.

My current approach to obtaining the requisite information changes (mutilates?) the EclipseGrid API as outlined below. I am not too happy about this because it adds data members to the class that have a very specialised purpose and are hardly ever used/useful in a real simulation case. Any other suggestions for where to put these (or similar) query methods? Could/should we key this off the SimulationConfig?

---

Changes to EclipseGrid.hpp

diff --git a/lib/eclipse/include/opm/parser/eclipse/EclipseState/Grid/EclipseGrid.hpp b/lib/eclipse/include/opm/parser/eclipse/EclipseState/Grid/EclipseGrid.hpp
index 3f59f489..193b935c 100644
--- a/lib/eclipse/include/opm/parser/eclipse/EclipseState/Grid/EclipseGrid.hpp
+++ b/lib/eclipse/include/opm/parser/eclipse/EclipseState/Grid/EclipseGrid.hpp
@@ -104,6 +104,8 @@ namespace Opm {
         */
         bool circle( ) const;
         bool isPinchActive( ) const;
+        bool isRadial( ) const;
+        double getInnerRadius( ) const;
         double getPinchThresholdThickness( ) const;
         PinchMode::ModeEnum getPinchOption( ) const;
         PinchMode::ModeEnum getMultzOption( ) const;
@@ -187,6 +189,8 @@ namespace Opm {
         PinchMode::ModeEnum m_multzMode;
         mutable std::vector< int > activeMap;
         bool m_circle = false;
+        bool m_isradial = false;
+        double m_inrad = 0.0;

Changes to EclipseGrid.cpp

diff --git a/lib/eclipse/EclipseState/Grid/EclipseGrid.cpp b/lib/eclipse/EclipseState/Grid/EclipseGrid.cpp
index ab89063e..aef395a2 100644
--- a/lib/eclipse/EclipseState/Grid/EclipseGrid.cpp
+++ b/lib/eclipse/EclipseState/Grid/EclipseGrid.cpp
@@ -187,6 +187,14 @@ namespace Opm {
         return this->m_circle;
     }

+    bool EclipseGrid::isRadial( ) const{
+        return this->m_isradial;
+    }
+
+    double EclipseGrid::getInnerRadius( ) const{
+       return this->m_inrad;
+    }
+
     void EclipseGrid::initGrid( const std::array<int, 3>& dims, const Deck& deck) {
         if (deck.hasKeyword<ParserKeywords::RADIAL>()) {
             initCylindricalGrid( dims, deck );
@@ -402,7 +410,8 @@ namespace Opm {
                 std::vector<double> tj(dims[1] + 1);
                 double z1 = *std::min_element( zcorn.begin() , zcorn.end());
                 double z2 = *std::max_element( zcorn.begin() , zcorn.end());
-                ri[0] = deck.getKeyword<ParserKeywords::INRAD>().getRecord(0).getItem(0).getSIDouble( 0 );
+                this->m_inrad = ri[0] =
+                    deck.getKeyword<ParserKeywords::INRAD>().getRecord(0).getItem(0).getSIDouble( 0 );
                 for (int i = 1; i <= dims[0]; i++)
                     ri[i] = ri[i - 1] + drv[i - 1];

@@ -434,6 +443,7 @@ namespace Opm {
                 }
             }
             initCornerPointGrid( dims , coord, zcorn, nullptr, nullptr);
+            m_isradial = true;
         }
     }

[joakim-hove]

Yes - this would of course mean that some very seldom used properties got disproportionally large room in the api/class - but of course that is a bit the way things are. I would merge this - some possible alternatives:

1. The m_isradial and the m_iscircle could be merged in a bitmask like structure, which could also embody the alternatives cornerpoint and block-centered?
2. The inner radius could be calculated from a i == 0 cell (that would require an api extension to get the coordinates of a cell though. Based on: #1097 this could be implemted as:

const auto& p = grid.getCornerPos(0,0,0,0);
double inrad = sqrt( p[0]*p[0] + p[1]*p[1] );