Add Experimental Edge Conformal CPG Processing Mode

Initially supported fully only for the 'UnstructuredGrid'
(create_grid_cornerpoint() constructor, PolyhedralGrid wrapper),
this commit introduces amended logic that aims to create edge
conformal cell complexes.  Specifically, the faces of the
'processed_grid' structure from 'process_grdecl()' contains *all*
vertices along the pillars, not just the top and bottom vertices of
the traditional corner-point description.  Moreover, the
create_grid_cornerpoint() will insert additional vertices in the top
and bottom *faces* when the edge conformal mode is activated.  This
aspect is not yet available in the CpGrid wrapper, though the
vertices along the pillars will be included in that case too.

Client code may activate this mode by setting the 'edge_conformal'
flag to 'true'.  The flag is a bool (i.e., _Bool) in both the C and
the C++ code, but transmitted as an 'int' across the language
boundary.

At this time, the mode is intended to support geo-mechanical
workflows and should typically not be enabled in production runs of
regular reservoir simulations.

While here, split some long lines and mark objects 'const' where
possible.

Author: hnil

CMakeLists_files.cmake

@@ -58,6 +58,7 @@ list (APPEND MAIN_SOURCE_FILES
   opm/grid/cornerpoint_grid.c
   opm/grid/cpgpreprocess/facetopology.c
   opm/grid/cpgpreprocess/geometry.c
+  opm/grid/cpgpreprocess/make_edge_conformal.cpp
   opm/grid/cpgpreprocess/preprocess.c
   opm/grid/cpgpreprocess/uniquepoints.c
   opm/grid/UnstructuredGrid.c
@@ -247,6 +248,7 @@ list (APPEND PUBLIC_HEADER_FILES
   opm/grid/cornerpoint_grid.h
   opm/grid/cpgpreprocess/facetopology.h
   opm/grid/cpgpreprocess/geometry.h
+  opm/grid/cpgpreprocess/make_edge_conformal.hpp
   opm/grid/cpgpreprocess/preprocess.h
   opm/grid/cpgpreprocess/uniquepoints.h
   opm/grid/transmissibility/trans_tpfa.h

opm/grid/CpGrid.hpp

@@ -224,58 +224,115 @@ namespace Dune
 #if HAVE_ECL_INPUT
         /// Read the Eclipse grid format ('grdecl').
         ///
-        /// \return Vector of global indices to the cells which have
-        ///         been removed in the grid processing due to small pore volume. Function only returns
-        ///         indices on rank 0, the vector is empty of other ranks.
-        /// \param ecl_grid the high-level object from opm-parser which represents the simulation's grid
-        ///        In a parallel run this may be a nullptr on all ranks but rank zero.
-        /// \param ecl_state the object from opm-parser provide information regarding to pore volume, NNC,
-        ///        aquifer information when ecl_state is available. NNC and aquifer connection
-        ///        information will also be updated during the function call when available and necessary.
-        /// \param periodic_extension if true, the grid will be (possibly) refined, so that
-        ///        intersections/faces along i and j boundaries will match those on the other
-        ///        side. That is, i- faces will match i+ faces etc.
-        /// \param turn_normals if true, all normals will be turned. This is intended for handling inputs with wrong orientations.
-        /// \param clip_z if true, the grid will be clipped so that the top and bottom will be planar.
-        /// \param pinchActive Force specific pinch behaviour. If true a face will connect two vertical cells, that are
-        ///           topological connected, even if there are cells with zero volume between them. If false these
-        ///           cells will not be connected despite their faces coinciding.
+        /// \return Vector of global indices to the cells which have been
+        /// removed in the grid processing due to small pore volume.
+        /// Function only returns indices on rank 0, the vector is empty of
+        /// other ranks.
+        ///
+        /// \param[in] ecl_grid the high-level object from opm-parser which
+        /// represents the simulation's grid.  In a parallel run this may be
+        /// a nullptr on all ranks but rank zero.
+        ///
+        /// \param[in,out] ecl_state the object from opm-parser provide
+        /// information regarding to pore volume, NNC, aquifer information
+        /// when ecl_state is available.  NNC and aquifer connection
+        /// nformation will also be updated during the function call when
+        /// available and necessary.
+        ///
+        /// \param[in] periodic_extension Whether or not to process the
+        /// resulting grid in order to have intersections/faces along i and
+        /// j boundaries match those on the other side. That is, i- faces
+        /// will match i+ faces etc.
+        ///
+        /// \param[in] turn_normals Whether or not to turn all normals.
+        /// This is intended for handling inputs with wrong orientations.
+        ///
+        /// \param[in] clip_z Whether or not to clip the result grid in
+        /// order to have planar top and bottom surfaces.
+        ///
+        /// \param[in] pinchActive Whether or not to force specific pinch
+        /// behaviour.  If set, a face will connect two vertical cells, that
+        /// are topological connected, even if there are cells with zero
+        /// volume between them. If false these cells will not be connected
+        /// despite their faces coinciding.
+        ///
+        /// \param[in] edge_conformal Whether or not to construct an
+        /// edge-conformal grid.  Typically useful in geo-mechanical
+        /// applications.
         std::vector<std::size_t>
         processEclipseFormat(const Opm::EclipseGrid* ecl_grid,
                              Opm::EclipseState* ecl_state,
-                             bool periodic_extension, bool turn_normals, bool clip_z,
-                             bool pinchActive);
+                             bool periodic_extension,
+                             bool turn_normals,
+                             bool clip_z,
+                             bool pinchActive,
+                             bool edge_conformal);
 
         /// Read the Eclipse grid format ('grdecl').
         ///
-        /// Pinch behaviour is determind from the parameter ecl_grid. If ecl_grid is a nullptr or PINCH was specified for
-        /// the grid, then a face will connect two vertical cells, that are topological connected, even if there are
-        /// cells with zero volume between them, Otherwise these cells will not be connected despite their faces coinciding.
-        ///
-        /// \return Vector of global indices to the cells which have
-        ///         been removed in the grid processing due to small pore volume. Function only returns
-        ///         indices on rank 0, the vector is empty of other ranks.
-        /// \param ecl_grid the high-level object from opm-parser which represents the simulation's grid
-        ///        In a parallel run this may be a nullptr on all ranks but rank zero.
-        /// \param ecl_state the object from opm-parser provide information regarding to pore volume, NNC,
-        ///        aquifer information when ecl_state is available. NNC and aquifer connection
-        ///        information will also be updated during the function call when available and necessary.
-        /// \param periodic_extension if true, the grid will be (possibly) refined, so that
-        ///        intersections/faces along i and j boundaries will match those on the other
-        ///        side. That is, i- faces will match i+ faces etc.
-        /// \param turn_normals if true, all normals will be turned. This is intended for handling inputs with wrong orientations.
-        /// \param clip_z if true, the grid will be clipped so that the top and bottom will be planar.
+        /// Pinch behaviour is determind from the parameter ecl_grid. If
+        /// ecl_grid is a nullptr or PINCH was specified for the grid, then
+        /// a face will connect two vertical cells, that are topological
+        /// connected, even if there are cells with zero volume between
+        /// them, Otherwise these cells will not be connected despite their
+        /// faces coinciding.
+        ///
+        /// \return Vector of global indices to the cells which have been
+        /// removed in the grid processing due to small pore
+        /// volume. Function only returns indices on rank 0, the vector is
+        /// empty of other ranks.
+        ///
+        /// \param[in] ecl_grid the high-level object from opm-parser which
+        /// represents the simulation's grid.  In a parallel run this may be
+        /// a nullptr on all ranks but rank zero.
+        ///
+        /// \param[in,out] ecl_state the object from opm-parser provide
+        /// information regarding to pore volume, NNC, aquifer information
+        /// when ecl_state is available.  NNC and aquifer connection
+        /// nformation will also be updated during the function call when
+        /// available and necessary.
+        ///
+        /// \param[in] periodic_extension Whether or not to process the
+        /// resulting grid in order to have intersections/faces along i and
+        /// j boundaries match those on the other side. That is, i- faces
+        /// will match i+ faces etc.
+        ///
+        /// \param[in] turn_normals Whether or not to turn all normals.
+        /// This is intended for handling inputs with wrong orientations.
+        ///
+        /// \param[in] clip_z Whether or not to clip the result grid in
+        /// order to have planar top and bottom surfaces.
+        ///
+        /// \param[in] edge_conformal Whether or not to construct an
+        /// edge-conformal grid.  Typically useful in geo-mechanical
+        /// applications.
         std::vector<std::size_t>
         processEclipseFormat(const Opm::EclipseGrid* ecl_grid,
                              Opm::EclipseState* ecl_state,
-                             bool periodic_extension, bool turn_normals = false, bool clip_z = false);
-
-#endif
+                             bool periodic_extension,
+                             bool turn_normals = false,
+                             bool clip_z = false,
+                             bool edge_conformal = false);
+#endif // HAVE_ECL_INPUT
 
         /// Read the Eclipse grid format ('grdecl').
-        /// \param input_data the data in grdecl format, declared in preprocess.h.
-        /// \param remove_ij_boundary if true, will remove (i, j) boundaries. Used internally.
-        void processEclipseFormat(const grdecl& input_data, bool remove_ij_boundary, bool turn_normals = false);
+        ///
+        /// \param[in] input_data the data in grdecl format, declared in
+        /// preprocess.h.
+        ///
+        /// \param[in] remove_ij_boundary if true, will remove (i, j)
+        /// boundaries. Used internally.
+        ///
+        /// \param[in] turn_normals if true, all normals will be turned. This is
+        /// intended for handling inputs with wrong orientations.
+        ///
+        /// \param[in] edge_conformal Whether or not to construct an
+        /// edge-conformal grid.  Typically useful in geo-mechanical
+        /// applications.
+        void processEclipseFormat(const grdecl& input_data,
+                                  bool remove_ij_boundary,
+                                  bool turn_normals = false,
+                                  bool edge_conformal = false);
 
         //@}
 

opm/grid/GridManager.cpp

@@ -40,18 +40,20 @@ namespace Opm
 
 #if HAVE_ECL_INPUT
     /// Construct a 3d corner-point grid from a deck.
-    GridManager::GridManager(const Opm::EclipseGrid& inputGrid)
+    GridManager::GridManager(const Opm::EclipseGrid& inputGrid,
+                             const bool edge_conformal)
         : ug_(0)
     {
-        initFromEclipseGrid(inputGrid, std::vector<double>());
+        initFromEclipseGrid(inputGrid, std::vector<double>(), edge_conformal);
     }
 
 
     GridManager::GridManager(const Opm::EclipseGrid& inputGrid,
-                             const std::vector<double>& poreVolumes)
+                             const std::vector<double>& poreVolumes,
+                             const bool edge_conformal)
         : ug_(0)
     {
-        initFromEclipseGrid(inputGrid, poreVolumes);
+        initFromEclipseGrid(inputGrid, poreVolumes, edge_conformal);
     }
 #endif
 
@@ -133,9 +135,10 @@ namespace Opm
 #if HAVE_ECL_INPUT
     // Construct corner-point grid from EclipseGrid.
     void GridManager::initFromEclipseGrid(const Opm::EclipseGrid& inputGrid,
-                                          const std::vector<double>& poreVolumes)
+                                          const std::vector<double>& poreVolumes,
+                                          const bool edge_conformal)
     {
-        struct grdecl g;
+        grdecl g{};
 
         g.dims[0] = inputGrid.getNX();
         g.dims[1] = inputGrid.getNY();
@@ -166,7 +169,7 @@ namespace Opm
                        actnum, opmfil, zcorn.data());
         }
 
-        ug_ = create_grid_cornerpoint(&g, z_tolerance);
+        ug_ = create_grid_cornerpoint(&g, z_tolerance, edge_conformal);
         if (!ug_) {
             OPM_THROW(std::runtime_error, "Failed to construct grid.");
         }

opm/grid/GridManager.hpp

@@ -47,15 +47,17 @@ class EclipseGrid;
 
 #if HAVE_ECL_INPUT
         /// Construct a grid from an EclipseState::EclipseGrid instance.
-        explicit GridManager(const EclipseGrid& inputGrid);
+        explicit GridManager(const EclipseGrid& inputGrid,
+                             bool edge_conformal = false);
 
         /// Construct a grid from an EclipseState::EclipseGrid instance,
         /// giving an explicit set of pore volumes to be used for MINPV
         /// considerations.
         /// \input[in] eclipseGrid    encapsulates a corner-point grid given from a deck
         /// \input[in] poreVolumes    one element per logical cartesian grid element
         GridManager(const EclipseGrid& inputGrid,
-                    const std::vector<double>& poreVolumes);
+                    const std::vector<double>& poreVolumes,
+                    bool edge_conformal);
 #endif
 
         /// Construct a 2d cartesian grid with cells of unit size.
@@ -92,7 +94,8 @@ class EclipseGrid;
 #if HAVE_ECL_INPUT
         // Construct corner-point grid from EclipseGrid.
         void initFromEclipseGrid(const EclipseGrid& inputGrid,
-                                 const std::vector<double>& poreVolumes);
+                                 const std::vector<double>& poreVolumes,
+                                 bool edge_conformal);
 #endif
 
         // The managed UnstructuredGrid.