[LGR] Aquifer cells/connections refinement is not supported or ignored for CpGrid

CMakeLists_files.cmake

@@ -112,6 +112,7 @@ list(APPEND TEST_SOURCE_FILES
   tests/cpgrid/lgr/adapt_cpgrid_test.cpp
   tests/cpgrid/lgr/addLgrs_in_allActiveCartesianGrid_test.cpp
   tests/cpgrid/lgr/addLgrsOnDistributedGrid_test.cpp
+  tests/cpgrid/lgr/aquifer_cells_and_conn_not_refined_test.cpp
   tests/cpgrid/lgr/autoRefine_test.cpp
   tests/cpgrid/lgr/global_refine_test.cpp
   tests/cpgrid/lgr/id_entity_entityrep_test.cpp

opm/grid/cpgrid/CpGrid.cpp

@@ -1773,18 +1773,8 @@ template cpgrid::Entity<1> createEntity(const CpGrid&, int, bool); // needed in
 
 bool CpGrid::mark(int refCount, const cpgrid::Entity<0>& element, bool throwOnFailure)
 {
-    // Throw if element has a neighboring cell from a different level.
-    // E.g., a coarse cell touching the boundary of an LGR, or
-    // a refined cell with a coarser/finner neighboring cell.
-    for (const auto& intersection : Dune::intersections(leafGridView(), element)){
-        if (intersection.neighbor() && (intersection.outside().level() != element.level()) && (element.level()==0)) {
-            // Refinement of cells at LGR boundaries is not supported, yet.
-            if (throwOnFailure)
-                OPM_THROW(std::invalid_argument, "Refinement of cells at LGR boundaries is not supported, yet.");
-            else
-                return false;
-        }
-    }
+    Opm::Lgr::throwIfIsCoarseAtLgrBoundary(leafGridView(), element, throwOnFailure);
+
     // For serial run, mark elements also in the level they were born.
     std::optional<bool> mark0;
     if(currentData().size()>1) {
@@ -1807,6 +1797,19 @@ int CpGrid::getMark(const cpgrid::Entity<0>& element) const
 
 bool CpGrid::preAdapt()
 {
+    // Aquifer data is only in level zero grid
+    auto levelZeroMarkGetter = [&](const cpgrid::Entity<0>& element)
+    {
+        return current_data_->front()->getMark(element); // The mark gets propagated to all lower levels.
+    };
+
+    // Marked aquifer cells or connections will be ignored in the refinement process
+    Opm::Lgr::throwIfAquCellOrConnHasBeenMarked(/* levelZeroData = */ *currentData().front(),
+                                                /* levelZeroView = */ levelGridView(0),
+                                                /* levelZeroAquiferCells = */ currentData().front()->sortedNumAquiferCells(),
+                                                levelZeroMarkGetter,
+                                                /* throwOnFailure = */ false);
+
     // Check if elements in pre-adapt existing grids have been marked for refinment.
     // Serial run: currentData() = data_. Parallel run: currentData() = distributed_data_.
     bool isPreAdapted = false; // 0
@@ -2724,6 +2727,19 @@ void CpGrid::addLgrsUpdateLeafView(const std::vector<std::array<int,3>>& cells_p
         }
         tmp_maxLevel = tmp_maxLevel + startIJK_vec_parent_grid.size(); // Update the maxLevel
 
+        // Aquifer data is only in level zero grid
+        auto levelZeroMarkGetter = [&](const cpgrid::Entity<0>& element)
+        {
+            return current_data_->front()->getMark(element); // The mark gets propagated to all lower levels.
+        };
+
+        // If there are marked aquifer cells or connections, throw
+        Opm::Lgr::throwIfAquCellOrConnHasBeenMarked(/* levelZeroData = */ *currentData().front(),
+                                                    /* levelZeroView = */ levelGridView(0),
+                                                    /* levelZeroAquiferCells = */ currentData().front()->sortedNumAquiferCells(),
+                                                    levelZeroMarkGetter,
+                                                    /* throwOnFailure = */ true);
+
         //   2. Create the corresponding LGRs. and  3. Update the leaf grid view.
         refineAndUpdateGrid(cells_per_dim_vec_parent_grid,
                             assignRefinedLevel,

opm/grid/cpgrid/LgrHelpers.cpp

@@ -2273,5 +2273,22 @@ void containsEightDifferentCorners(const std::array<int,8>& cell_to_point)
     }
 }
 
+bool throwIfAquiferCell(const std::vector<int>& levelAquiferCells,
+                        const Dune::cpgrid::Entity<0>& element,
+                        bool throwOnFailure)
+{
+    if (!levelAquiferCells.empty()) {
+        bool isAquifer = std::ranges::find(levelAquiferCells, element.getLevelElem().index()) != levelAquiferCells.end();
+        if (isAquifer) {
+            if (throwOnFailure)
+                OPM_THROW(std::invalid_argument, "Refinement of aquifer cells is not supported, yet.");
+            else 
+                return false;
+        }
+    }
+    return true;
+}
+
+
 } // namespace Lgr
 } // namespace Opm

opm/grid/cpgrid/LgrHelpers.hpp

@@ -775,6 +775,54 @@ bool compatibleSubdivisions(const std::vector<std::array<int,3>>& cells_per_dim_
 
 void containsEightDifferentCorners(const std::array<int,8>& cell_to_point);
 
+template <class LevelZeroData, class LevelZeroView, typename MarkGetter>
+bool throwIfAquCellOrConnHasBeenMarked(const LevelZeroData& levelZeroData,
+                                       const LevelZeroView& levelZeroView,
+                                       const std::vector<int>& levelZeroAquiferCells,
+                                       const MarkGetter& markGetter,
+                                       bool throwOnFailure)
+{
+    for (const auto& aquCellIdx : levelZeroAquiferCells) {
+        const auto& aquElem = Dune::cpgrid::Entity<0>(levelZeroData, aquCellIdx, true);
+
+        if (markGetter(aquElem) == 1) {
+            if (throwOnFailure)
+                OPM_THROW(std::invalid_argument, "Refinement of cells connected to aquifers is not supported, yet.");
+            else
+                return false;
+        }
+
+        for (const auto& intersection : Dune::intersections(levelZeroView, aquElem)){
+            if (intersection.neighbor()) {
+                if (markGetter(intersection.outside()) == 1) {
+                    if (throwOnFailure)
+                        OPM_THROW(std::invalid_argument, "Refinement of cells connected to aquifers is not supported, yet.");
+                    else
+                        return false;
+                }
+            }
+        }
+    }
+    return true;
+}
+
+template <class LeafView>
+bool throwIfIsCoarseAtLgrBoundary(const LeafView& leafView,
+                                  const Dune::cpgrid::Entity<0>& element,
+                                  bool throwOnFailure)
+{
+    for (const auto& intersection : Dune::intersections(leafView, element)){
+        if (intersection.neighbor() && (intersection.outside().level() != element.level()) && (element.level() == 0)) {
+            // Refinement of cells at LGR boundaries is not supported, yet.
+            if (throwOnFailure)
+                OPM_THROW(std::invalid_argument, "Refinement of cells at LGR boundaries is not supported, yet.");
+            else
+                return false;
+        }
+    }
+    return true;
+}
+
 } // namespace Lgr
 } // namespace Opm
 

tests/cpgrid/lgr/LgrChecks.hpp

@@ -95,6 +95,11 @@ void checkVertexAndFaceIndexAreNonNegative(const Dune::CpGrid& grid);
 void checkFaceHas4VerticesAndMax2NeighboringCells(const Dune::CpGrid& grid,
                                                   const std::vector<std::shared_ptr<Dune::cpgrid::CpGridData>>& data);
 
+#if HAVE_OPM_COMMON
+void createGridFromDeckString(Dune::CpGrid& grid,
+                              const std::string& deck_string);
+#endif
+
 #if HAVE_OPM_COMMON
 void createGridAndAddLgrs(Dune::CpGrid& grid,
                           const std::string& deck_string,
@@ -534,23 +539,28 @@ void Opm::checkFaceHas4VerticesAndMax2NeighboringCells(const Dune::CpGrid& grid,
 }
 
 #if HAVE_OPM_COMMON
-void Opm::createGridAndAddLgrs(Dune::CpGrid& grid,
-                               const std::string& deck_string,
-                               const std::vector<std::array<int, 3>>& cells_per_dim_vec,
-                               const std::vector<std::array<int, 3>>& startIJK_vec,
-                               const std::vector<std::array<int, 3>>& endIJK_vec,
-                               const std::vector<std::string>& lgr_name_vec)
+void Opm::createGridFromDeckString(Dune::CpGrid& grid,
+                                   const std::string& deck_string)
 {
     Parser parser;
     const auto deck = parser.parseString(deck_string);
     EclipseState ecl_state(deck);
     EclipseGrid eclipse_grid = ecl_state.getInputGrid();
 
     grid.processEclipseFormat(&eclipse_grid, &ecl_state, false, false, false);
+}
+#endif
 
+void Opm::createGridAndAddLgrs(Dune::CpGrid& grid,
+                               const std::string& deck_string,
+                               const std::vector<std::array<int, 3>>& cells_per_dim_vec,
+                               const std::vector<std::array<int, 3>>& startIJK_vec,
+                               const std::vector<std::array<int, 3>>& endIJK_vec,
+                               const std::vector<std::string>& lgr_name_vec)
+{
+    Opm::createGridFromDeckString(grid, deck_string);
     grid.addLgrsUpdateLeafView(cells_per_dim_vec, startIJK_vec, endIJK_vec, lgr_name_vec);
 }
-#endif
 
 void Opm::createGridAndAddLgrs(Dune::CpGrid& grid,
                                const std::array<double, 3>& cell_sizes,