using pinched grid for numerical aquifers.

examples/grdecl2vtu.cpp

@@ -143,7 +143,7 @@ try
     {
         const int* actnum = deck.hasKeyword("ACTNUM") ? deck.getKeyword("ACTNUM").getIntData().data() : nullptr;
         Opm::EclipseGrid ecl_grid(deck , actnum);
-        grid.processEclipseFormat(&ecl_grid, false);
+        grid.processEclipseFormat(&ecl_grid, nullptr, false);
     }
 
     VTKWriter<CpGrid::LeafGridView> vtkwriter(grid.leafGridView());

opm/grid/CpGrid.hpp

@@ -259,16 +259,17 @@ namespace Dune
         ///         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 poreVolume pore volumes for use in MINPV processing, if asked for in deck
-        std::vector<std::size_t> processEclipseFormat(const Opm::EclipseGrid* ecl_grid, bool periodic_extension, bool turn_normals = false, bool clip_z = false,
-                                                      const std::vector<double>& poreVolume = std::vector<double>(),
-                                                      const Opm::NNC& = Opm::NNC(),
-                                                      const std::unordered_map<size_t, double>& aquifer_cell_volumes = std::unordered_map<size_t, double>());
+        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
 
         /// Read the Eclipse grid format ('grdecl').

opm/grid/cpgrid/CpGrid.cpp

@@ -493,7 +493,10 @@ CpGrid::scatterGrid(EdgeWeightMethod method,
         g.coord = &coord[0];
         g.zcorn = &zcorn[0];
         g.actnum = &actnum[0];
-        current_view_data_->processEclipseFormat(g, {}, 0.0, false, false);
+        using NNCMap = std::set<std::pair<int, int>>;
+        using NNCMaps = std::array<NNCMap, 2>;
+        NNCMaps nnc;
+        current_view_data_->processEclipseFormat(g, nullptr, nnc, 0.0, false, false);
         // global grid only on rank 0
         current_view_data_->ccobj_.broadcast(current_view_data_->logical_cartesian_size_.data(),
                                              current_view_data_->logical_cartesian_size_.size(),
@@ -522,15 +525,12 @@ CpGrid::scatterGrid(EdgeWeightMethod method,
 
 #if HAVE_ECL_INPUT
     std::vector<std::size_t> CpGrid::processEclipseFormat(const Opm::EclipseGrid* ecl_grid,
+                                                          Opm::EclipseState* ecl_state,
                                                           bool periodic_extension,
-                                                          bool turn_normals, bool clip_z,
-                                                          const std::vector<double>& poreVolume,
-                                                          const Opm::NNC& nncs,
-                                                          const std::unordered_map<size_t, double>& aquifer_cell_volumes)
+                                                          bool turn_normals, bool clip_z)
     {
-        auto removed_cells = current_view_data_->processEclipseFormat(ecl_grid, periodic_extension,
-                                                                      turn_normals, clip_z,
-                                                                      poreVolume, nncs, aquifer_cell_volumes);
+        auto removed_cells = current_view_data_->processEclipseFormat(ecl_grid, ecl_state, periodic_extension,
+                                                                      turn_normals, clip_z);
         current_view_data_->ccobj_.broadcast(current_view_data_->logical_cartesian_size_.data(),
                                              current_view_data_->logical_cartesian_size_.size(),
                                              0);
@@ -541,7 +541,10 @@ CpGrid::scatterGrid(EdgeWeightMethod method,
     void CpGrid::processEclipseFormat(const grdecl& input_data, double z_tolerance,
                                       bool remove_ij_boundary, bool turn_normals)
     {
-        current_view_data_->processEclipseFormat(input_data, {}, z_tolerance, remove_ij_boundary, turn_normals);
+        using NNCMap = std::set<std::pair<int, int>>;
+        using NNCMaps = std::array<NNCMap, 2>;
+        NNCMaps nnc;
+        current_view_data_->processEclipseFormat(input_data, nullptr, nnc, z_tolerance, remove_ij_boundary, turn_normals);
         current_view_data_->ccobj_.broadcast(current_view_data_->logical_cartesian_size_.data(),
                                              current_view_data_->logical_cartesian_size_.size(),
                                              0);

opm/grid/cpgrid/CpGridData.hpp

@@ -88,7 +88,10 @@
 #include "Entity2IndexDataHandle.hpp"
 #include "DataHandleWrappers.hpp"
 #include "GlobalIdMapping.hpp"
-
+namespace Opm
+{
+class EclipseState;
+}
 namespace Dune
 {
 class CpGrid;
@@ -197,24 +200,29 @@ class CpGridData
     /// Read the Eclipse grid format ('grdecl').
     /// \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 poreVolume pore volumes for use in MINPV processing, if asked for in deck
-    std::vector<std::size_t> processEclipseFormat(const Opm::EclipseGrid* ecl_grid, bool periodic_extension, bool turn_normals = false, bool clip_z = false,
-                                                  const std::vector<double>& poreVolume = std::vector<double>(), const Opm::NNC& nncs = Opm::NNC(),
-                                                  const std::unordered_map<size_t, double>& aquifer_cell_volumes = std::unordered_map<size_t, double>());
+    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
 
     /// Read the Eclipse grid format ('grdecl').
     /// \param input_data the data in grdecl format, declared in preprocess.h.
+    /// \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 nnc is the non-neighboring connections
     /// \param z_tolerance points along a pillar that are closer together in z
     ///        coordinate than this parameter, will be replaced by a single point.
     /// \param remove_ij_boundary if true, will remove (i, j) boundaries. Used internally.
-    void processEclipseFormat(const grdecl& input_data, const std::array<std::set<std::pair<int, int>>, 2>& nnc, double z_tolerance, bool remove_ij_boundary, bool turn_normals = false,
-                              const std::unordered_map<size_t, double>& aquifer_cell_volumes = std::unordered_map<size_t, double>());
+    void processEclipseFormat(const grdecl& input_data, Opm::EclipseState* ecl_state,
+                              std::array<std::set<std::pair<int, int>>, 2>& nnc, double z_tolerance, bool remove_ij_boundary, bool turn_normals = false);
 
 
     /// @brief

opm/grid/cpgrid/processEclipseFormat.cpp

@@ -40,6 +40,8 @@
 
 #include "CpGridData.hpp"
 
+#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
+
 #include <opm/grid/common/GeometryHelpers.hpp>
 #include <opm/grid/cpgrid/EntityRep.hpp>
 
@@ -114,19 +116,16 @@ namespace cpgrid
 {
 
 #if HAVE_ECL_INPUT
-    std::vector<std::size_t> CpGridData::processEclipseFormat(const Opm::EclipseGrid* ecl_grid_ptr, bool periodic_extension, bool turn_normals, bool clip_z,
-                                                              const std::vector<double>& poreVolume, const Opm::NNC& nncs,
-                                                              const std::unordered_map<size_t, double>& aquifer_cell_volumes)
+    std::vector<std::size_t> CpGridData::processEclipseFormat(const Opm::EclipseGrid* ecl_grid_ptr,
+                                                              Opm::EclipseState* ecl_state,
+                                                              bool periodic_extension, bool turn_normals, bool clip_z)
     {
         std::vector<std::size_t> removed_cells;
         if (ccobj_.rank() != 0 ) {
             // Store global grid only on rank 0
             return removed_cells;
         }
 
-        if (!ecl_grid_ptr)
-            OPM_THROW(std::logic_error, "We need a valid pointer to an eclipse grid on rank 0!");
-
         const Opm::EclipseGrid& ecl_grid = *ecl_grid_ptr;
         std::vector<double> coordData = ecl_grid.getCOORD();
         std::vector<int> actnumData = ecl_grid.getACTNUM();
@@ -146,7 +145,7 @@ namespace cpgrid
         Opm::MinpvProcessor::Result minpv_result;
 
         // Possibly process MINPV
-        if (!poreVolume.empty() && (ecl_grid.getMinpvMode() != Opm::MinpvMode::ModeEnum::Inactive)) {
+        if (ecl_state && (ecl_grid.getMinpvMode() != Opm::MinpvMode::ModeEnum::Inactive)) {
             Opm::MinpvProcessor mp(g.dims[0], g.dims[1], g.dims[2]);
             // Currently PINCH is always assumed to be active
             const size_t cartGridSize = g.dims[0] * g.dims[1] * g.dims[2];
@@ -156,6 +155,7 @@ namespace cpgrid
             }
             const double z_tolerance = ecl_grid.isPinchActive() ?  ecl_grid.getPinchThresholdThickness() : 0.0;
             const bool nogap = ecl_grid.getPinchGapMode() ==  Opm::PinchMode::ModeEnum::NOGAP;
+            const auto& poreVolume = ecl_state->fieldProps().porv(true);
             minpv_result = mp.process(thickness, z_tolerance, poreVolume, ecl_grid.getMinpvVector(), actnumData, false, zcornData.data(), nogap);
             if (minpv_result.nnc.size() > 0) {
                 this->zcorn = zcornData;
@@ -168,12 +168,15 @@ namespace cpgrid
             nnc_cells[PinchNNC].insert({cell1, cell2});
 
         // Add explicit NNCs.
-        for (const auto single_nnc : nncs.input()) {
-            // Repeated NNCs will only exist in the map once (repeated
-            // insertions have no effect). The code that computes the
-            // transmissibilities is responsible for ensuring repeated NNC
-            // transmissibilities are added.
-            nnc_cells[ExplicitNNC].insert({single_nnc.cell1, single_nnc.cell2});
+        if (ecl_state) {
+            const auto& nncs = ecl_state->getInputNNC();
+            for (const auto single_nnc : nncs.input()) {
+                // Repeated NNCs will only exist in the map once (repeated
+                // insertions have no effect). The code that computes the
+                // transmissibilities is responsible for ensuring repeated NNC
+                // transmissibilities are added.
+                nnc_cells[ExplicitNNC].insert({single_nnc.cell1, single_nnc.cell2});
+            }
         }
 
         // this variable is only required because getCellZvals() needs
@@ -231,10 +234,10 @@ namespace cpgrid
             grdecl new_g;
             addOuterCellLayer(g, new_coord, new_zcorn, new_actnum, new_g);
             // Make the grid.
-            processEclipseFormat(new_g, nnc_cells, z_tolerance, true, turn_normals, aquifer_cell_volumes);
+            processEclipseFormat(new_g, ecl_state, nnc_cells, z_tolerance, true, turn_normals);
         } else {
             // Make the grid.
-            processEclipseFormat(g, nnc_cells, z_tolerance, false, turn_normals, aquifer_cell_volumes);
+            processEclipseFormat(g, ecl_state, nnc_cells, z_tolerance, false, turn_normals);
         }
 
         return minpv_result.removed_cells;
@@ -246,8 +249,8 @@ namespace cpgrid
 
 
     /// Read the Eclipse grid format ('.grdecl').
-    void CpGridData::processEclipseFormat(const grdecl& input_data, const NNCMaps& nnc, double z_tolerance, bool remove_ij_boundary, bool turn_normals,
-                                          const std::unordered_map<size_t, double>& aquifer_cell_volumes)
+    void CpGridData::processEclipseFormat(const grdecl& input_data, Opm::EclipseState* ecl_state,
+                                          NNCMaps& nnc, double z_tolerance, bool remove_ij_boundary, bool turn_normals)
     {
         if( ccobj_.rank() != 0 )
         {
@@ -258,17 +261,43 @@ namespace cpgrid
         std::cout << "Processing eclipse data." << std::endl;
 #endif
         processed_grid output;
-        const size_t global_nc = input_data.dims[0] * input_data.dims[1] * input_data.dims[2];
-        std::vector<int> is_aquifer_cell(global_nc, 0);
-        for ([[maybe_unused]]const auto& [global_index, volume] : aquifer_cell_volumes) {
-            is_aquifer_cell[global_index] = 1;
+        if (ecl_state) {
+            const auto& aquifer = ecl_state->aquifer();
+            const auto aquifer_cell_volumes = aquifer.numericalAquifers().aquiferCellVolumes();
+            const size_t global_nc = input_data.dims[0] * input_data.dims[1] * input_data.dims[2];
+            std::vector<int> is_aquifer_cell(global_nc, 0);
+            for ([[maybe_unused]]const auto&[global_index, volume] : aquifer_cell_volumes) {
+                is_aquifer_cell[global_index] = 1;
+            }
+            process_grdecl(&input_data, z_tolerance, is_aquifer_cell.data(), &output);
+        } else {
+            process_grdecl(&input_data, z_tolerance, nullptr, &output);
         }
-        process_grdecl(&input_data, z_tolerance, is_aquifer_cell.data(), &output);
         if (remove_ij_boundary) {
             removeOuterCellLayer(output);
             // removeUnusedNodes(output);
         }
 
+        if (ecl_state) {
+            const auto& aquifer = ecl_state->aquifer();
+            if (aquifer.hasNumericalAquifer()) {
+                const size_t global_nc = input_data.dims[0] * input_data.dims[1] * input_data.dims[2];
+                std::vector<int> new_actnum(global_nc, 0);
+                for (int i = 0; i < output.number_of_cells; ++i) {
+                    new_actnum[output.local_cell_index[i]] = 1;
+                }
+                const auto& ecl_grid = ecl_state->getInputGrid();
+                const auto& fp = ecl_state->fieldProps();
+                aquifer.mutableNumericalAquifers().postProcessConnections(ecl_grid, new_actnum);
+                const auto& aquifer_nnc = aquifer.numericalAquifers().aquiferConnectionNNCs(ecl_grid, fp);
+                // We need to update the nnc in the ecl_state
+                ecl_state->appendInputNNC(aquifer_nnc);
+                for (const auto single_nnc : aquifer_nnc) {
+                    nnc[ExplicitNNC].insert({single_nnc.cell1, single_nnc.cell2});
+                }
+            }
+        }
+
         // Move data into the grid's structures.
 #ifdef VERBOSE
         std::cout << "Building topology." << std::endl;
@@ -282,10 +311,13 @@ namespace cpgrid
 #endif
         // here we need the cell volumes based on the active index order
         std::unordered_map<size_t, double> aquifer_cell_volumes_local;
-        for (auto nc = this->global_cell_.size(), i = 0*nc; i < nc; ++i) {
-            auto aquCellPos = aquifer_cell_volumes.find(this->global_cell_[i]);
-            if (aquCellPos != aquifer_cell_volumes.end()) {
-                aquifer_cell_volumes_local.emplace(i, aquCellPos->second);
+        if (ecl_state && ecl_state->aquifer().hasNumericalAquifer()) {
+            const auto& aquifer_cell_volumes = ecl_state->aquifer().numericalAquifers().aquiferCellVolumes();
+            for (auto nc = this->global_cell_.size(), i = 0 * nc; i < nc; ++i) {
+                auto aquCellPos = aquifer_cell_volumes.find(this->global_cell_[i]);
+                if (aquCellPos != aquifer_cell_volumes.end()) {
+                    aquifer_cell_volumes_local.emplace(i, aquCellPos->second);
+                }
             }
         }
         buildGeom(output, cell_to_face_, cell_to_point_, face_to_output_face, aquifer_cell_volumes_local, geometry_.geomVector(std::integral_constant<int,0>()),

tests/cpgrid/grid_nnc.cpp

@@ -64,15 +64,17 @@ struct Fixture
                   const std::vector<std::pair<int, int>>& ex_nb,
                   const bool use_deck_porv = false)
     {
-        Opm::EclipseState es(parser.parseFile(filename));
+        Opm::Deck deck = parser.parseFile(filename);
+        Opm::EclipseState es(deck);
+        es.appendInputNNC(nnc.input());
         std::vector<double> porv;
 
         if (use_deck_porv) {
             porv = es.fieldProps().porv(true);
         }
 
         Dune::CpGrid grid;
-        grid.processEclipseFormat(&es.getInputGrid(), false, false, false, porv, nnc);
+        grid.processEclipseFormat(&es.getInputGrid(), &es, false, false, false);
         const auto& gv = grid.leafGridView();
         ElementMapper elmap(gv, Dune::mcmgElementLayout());
         int elemcount = 0;

tests/test_cpgrid.cpp

@@ -158,13 +158,11 @@ int main(int argc, char** argv )
 
     Opm::Parser parser;
     const auto deck = parser.parseString(deckString);
-    std::vector<double> porv;
 
     Grid grid;
-    const int* actnum = deck.hasKeyword("ACTNUM") ? deck.getKeyword("ACTNUM").getIntData().data() : nullptr;
-    Opm::EclipseGrid ecl_grid(deck , actnum);
+    Opm::EclipseGrid ecl_grid(deck);
 
-    grid.processEclipseFormat(&ecl_grid, false, false, false, porv);
+    grid.processEclipseFormat(&ecl_grid, nullptr, false, false, false);
     testGrid( grid, "CpGrid_ecl", 8, 27 );
 #endif