Merge pull request #5505 from lisajulia/fix/COMPDAT-NLDD

Fix the usage of the COMPDAT keyword in an ACTIONX for the NLDD solver for sequential runs

opm/simulators/flow/BlackoilModelNldd.hpp

@@ -953,9 +953,10 @@ private:
             ? param.num_local_domains_
             : num_cells / default_cells_per_domain;
 
+        const auto& possibleFutureConnectionSet = this->model_.simulator().vanguard().schedule().getPossibleFutureConnections();
         return ::Opm::partitionCells(param.local_domain_partition_method_,
                                      num_domains,
-                                     grid.leafGridView(), wells, zoltan_ctrl);
+                                     grid.leafGridView(), wells, possibleFutureConnectionSet, zoltan_ctrl);
     }
 
     std::vector<int> reconstitutePartitionVector() const

opm/simulators/flow/GenericCpGridVanguard.cpp

@@ -216,7 +216,7 @@ doLoadBalance_(const Dune::EdgeWeightMethod             edgeWeightsMethod,
 #if HAVE_OPENCL || HAVE_ROCSPARSE || HAVE_CUDA
         if (partitionJacobiBlocks) {
             this->cell_part_ = this->grid_->
-                zoltanPartitionWithoutScatter(&wells, &possibleFutureConnections, faceTrans.data(),
+                zoltanPartitionWithoutScatter(&wells, possibleFutureConnections, faceTrans.data(),
                                               numJacobiBlocks,
                                               imbalanceTol);
         }
@@ -292,7 +292,7 @@ distributeGrid(const Dune::EdgeWeightMethod
                const bool                                            loadBalancerSet,
                const std::vector<double>&                            faceTrans,
                const std::vector<Well>&                              wells,
-               const std::unordered_map<std::string, std::set<std::array<int,3>>>& possibleFutureConnections,
+               const std::unordered_map<std::string, std::set<int>>& possibleFutureConnections,
                EclipseState&                                         eclState1,
                FlowGenericVanguard::ParallelWellStruct&              parallelWells)
 {
@@ -330,7 +330,7 @@ distributeGrid(const Dune::EdgeWeightMethod
                const bool                                            loadBalancerSet,
                const std::vector<double>&                            faceTrans,
                const std::vector<Well>&                              wells,
-               const std::unordered_map<std::string, std::set<std::array<int,3>>>& possibleFutureConnections,
+               const std::unordered_map<std::string, std::set<int>>& possibleFutureConnections,
                ParallelEclipseState*                                 eclState,
                FlowGenericVanguard::ParallelWellStruct&              parallelWells)
 {
@@ -350,13 +350,13 @@ distributeGrid(const Dune::EdgeWeightMethod
             : std::vector<int>{};
         //For this case, simple partitioning is selected automatically
         parallelWells =
-            std::get<1>(this->grid_->loadBalance(handle, parts, &wells, &possibleFutureConnections, ownersFirst,
+            std::get<1>(this->grid_->loadBalance(handle, parts, &wells, possibleFutureConnections, ownersFirst,
                                                  addCornerCells, overlapLayers));
     }
     else {
         parallelWells =
             std::get<1>(this->grid_->loadBalance(handle, edgeWeightsMethod,
-                                                 &wells, &possibleFutureConnections,
+                                                 &wells, possibleFutureConnections,
                                                  serialPartitioning,
                                                  faceTrans.data(), ownersFirst,
                                                  addCornerCells, overlapLayers,

opm/simulators/flow/GenericCpGridVanguard.hpp

@@ -174,7 +174,7 @@ private:
                         const bool                                            loadBalancerSet,
                         const std::vector<double>&                            faceTrans,
                         const std::vector<Well>&                              wells,
-                        const std::unordered_map<std::string, std::set<std::array<int,3>>>& possibleFutureConnections,
+                        const std::unordered_map<std::string, std::set<int>>& possibleFutureConnections,
                         EclipseState&                                         eclState,
                         FlowGenericVanguard::ParallelWellStruct&              parallelWells);
 
@@ -187,7 +187,7 @@ private:
                         const bool                                            loadBalancerSet,
                         const std::vector<double>&                            faceTrans,
                         const std::vector<Well>&                              wells,
-                        const std::unordered_map<std::string, std::set<std::array<int,3>>>& possibleFutureConnections,
+                        const std::unordered_map<std::string, std::set<int>>& possibleFutureConnections,
                         ParallelEclipseState*                                 eclState,
                         FlowGenericVanguard::ParallelWellStruct&              parallelWells);
 

opm/simulators/flow/partitionCells.cpp

@@ -108,6 +108,7 @@ public:
     template <class GridView, class Element>
     void buildLocalGraph(const GridView&                                       grid_view,
                          const std::vector<Opm::Well>&                         wells,
+                         const std::unordered_map<std::string, std::set<int>>& possibleFutureConnections,
                          const Opm::ZoltanPartitioningControl<Element>&        zoltan_ctrl);
 
     /// Partition rank's interior cells into non-overlapping domains using
@@ -180,6 +181,7 @@ private:
     template <typename Comm>
     void connectWells(const Comm                                            comm,
                       const std::vector<Opm::Well>&                         wells,
+                      const std::unordered_map<std::string, std::set<int>>& possibleFutureConnections,
                       const std::unordered_map<int, int>&                   g2l);
 };
 
@@ -196,9 +198,10 @@ ZoltanPartitioner::ZoltanPartitioner(const GridView&
 template <class GridView, class Element>
 void ZoltanPartitioner::buildLocalGraph(const GridView&                                       grid_view,
                                         const std::vector<Opm::Well>&                         wells,
+                                        const std::unordered_map<std::string, std::set<int>>& possibleFutureConnections,
                                         const Opm::ZoltanPartitioningControl<Element>&        zoltan_ctrl)
 {
-    this->connectWells(grid_view.comm(), wells, this->connectElements(grid_view, zoltan_ctrl));
+    this->connectWells(grid_view.comm(), wells, possibleFutureConnections, this->connectElements(grid_view, zoltan_ctrl));
 }
 
 template <class GridView, class Element>
@@ -284,6 +287,7 @@ ZoltanPartitioner::connectElements(const GridView&
 template <typename Comm>
 void ZoltanPartitioner::connectWells(const Comm                                            comm,
                                      const std::vector<Opm::Well>&                         wells,
+                                     const std::unordered_map<std::string, std::set<int>>& possibleFutureConnections,
                                      const std::unordered_map<int, int>&                   g2l)
 {
     auto distributedWells = 0;
@@ -301,6 +305,17 @@ void ZoltanPartitioner::connectWells(const Comm                          comm,
 
             cellIx.push_back(locPos->second);
         }
+        const auto possibleFutureConnectionSetIt = possibleFutureConnections.find(well.name());
+        if (possibleFutureConnectionSetIt != possibleFutureConnections.end()) {
+            for (auto& global_index : possibleFutureConnectionSetIt->second) {
+                auto locPos = g2l.find(global_index);
+                if (locPos == g2l.end()) {
+                    ++otherProc;
+                    continue;
+                }
+                cellIx.push_back(locPos->second);
+            }
+        }
 
         if ((otherProc > 0) && !cellIx.empty()) {
             ++distributedWells;
@@ -341,6 +356,7 @@ std::pair<std::vector<int>, int>
 partitionCellsZoltan(const int                                             num_domains,
                      const GridView&                                       grid_view,
                      const std::vector<Opm::Well>&                         wells,
+                     const std::unordered_map<std::string, std::set<int>>& possibleFutureConnections,
                      const Opm::ZoltanPartitioningControl<Element>&        zoltan_ctrl)
 {
     if (num_domains <= 1) {     // No partitioning => every cell in domain zero.
@@ -355,7 +371,7 @@ partitionCellsZoltan(const int                                      num_domains,
     }
 
     auto partitioner = ZoltanPartitioner { grid_view, zoltan_ctrl.local_to_global };
-    partitioner.buildLocalGraph(grid_view, wells, zoltan_ctrl);
+    partitioner.buildLocalGraph(grid_view, wells, possibleFutureConnections, zoltan_ctrl);
 
     return partitioner.partition(num_domains, grid_view, zoltan_ctrl);
 }
@@ -575,12 +591,13 @@ Opm::partitionCells(const std::string& method,
                     const int          num_local_domains,
                     const GridView&    grid_view,
                     [[maybe_unused]] const std::vector<Well>&                              wells,
+                    [[maybe_unused]] const std::unordered_map<std::string, std::set<int>>& possibleFutureConnections,
                     [[maybe_unused]] const ZoltanPartitioningControl<Element>&             zoltan_ctrl)
 {
     if (method == "zoltan") {
 #if HAVE_MPI && HAVE_ZOLTAN
 
-        return partitionCellsZoltan(num_local_domains, grid_view, wells, zoltan_ctrl);
+        return partitionCellsZoltan(num_local_domains, grid_view, wells, possibleFutureConnections, zoltan_ctrl);
 
 #else // !HAVE_MPI || !HAVE_ZOLTAN
 
@@ -672,6 +689,7 @@ Opm::partitionCellsSimple(const int num_cells, const int num_domains)
                         const std::remove_cv_t<std::remove_reference_t<        \
                         decltype(std::declval<Grid>().leafGridView())>>&,      \
                         const std::vector<Opm::Well>&,                         \
+                        const std::unordered_map<std::string, std::set<int>>&, \
                         const Opm::ZoltanPartitioningControl<                  \
                         typename std::remove_cv_t<std::remove_reference_t<     \
                         decltype(std::declval<Grid>().leafGridView())>>::template Codim<0>::Entity>&)

opm/simulators/flow/partitionCells.hpp

@@ -21,6 +21,7 @@
 #define OPM_ASPINPARTITION_HEADER_INCLUDED
 
 #include <functional>
+#include <set>
 #include <string>
 #include <utility>
 #include <vector>
@@ -92,6 +93,7 @@ partitionCells(const std::string&                        method,
                const int                                             num_local_domains,
                const GridView&                                       grid_view,
                const std::vector<Well>&                              wells,
+               const std::unordered_map<std::string, std::set<int>>& possibleFutureConnections,
                const ZoltanPartitioningControl<Element>&             zoltan_ctrl);
 
 /// Read a partitioning from file, assumed to contain one number per cell, its partition number.

opm/simulators/linalg/ISTLSolverBda.cpp

@@ -78,7 +78,7 @@ void BdaSolverInfo<Matrix,Vector>::
 prepare(const Grid& grid,
         const Dune::CartesianIndexMapper<Grid>& cartMapper,
         const std::vector<Well>& wellsForConn,
-        const std::unordered_map<std::string, std::set<std::array<int,3>>>& possibleFutureConnections,
+        const std::unordered_map<std::string, std::set<int>>& possibleFutureConnections,
         const std::vector<int>& cellPartition,
         const std::size_t nonzeroes,
         const bool useWellConn)
@@ -246,7 +246,7 @@ using BV = Dune::BlockVector<Dune::FieldVector<Scalar,Dim>>;
     prepare(const Grid&,                                           \
             const Dune::CartesianIndexMapper<Grid>&,               \
             const std::vector<Well>&,                              \
-            const std::unordered_map<std::string, std::set<std::array<int,3>>>&, \
+            const std::unordered_map<std::string, std::set<int>>&, \
             const std::vector<int>&,                               \
             const std::size_t, const bool);
 using PolyHedralGrid3D = Dune::PolyhedralGrid<3, 3>;

opm/simulators/linalg/ISTLSolverBda.hpp

@@ -60,7 +60,7 @@ struct BdaSolverInfo
     void prepare(const Grid& grid,
                  const Dune::CartesianIndexMapper<Grid>& cartMapper,
                  const std::vector<Well>& wellsForConn,
-                 const std::unordered_map<std::string, std::set<std::array<int,3>>>& possibleFutureConnections,
+                 const std::unordered_map<std::string, std::set<int>>& possibleFutureConnections,
                  const std::vector<int>& cellPartition,
                  const std::size_t nonzeroes,
                  const bool useWellConn);

opm/simulators/linalg/findOverlapRowsAndColumns.hpp

@@ -46,7 +46,7 @@ namespace detail
     /// \param useWellConn Boolean that is true when UseWellContribusion is true
     /// \param wellGraph Cell IDs of well cells stored in a graph.
     template<class Grid, class CartMapper, class W>
-    void setWellConnections(const Grid& grid, const CartMapper& cartMapper, const W& wells, const std::unordered_map<std::string, std::set<std::array<int,3>>>& possibleFutureConnections, bool useWellConn, std::vector<std::set<int>>& wellGraph, int numJacobiBlocks)
+    void setWellConnections(const Grid& grid, const CartMapper& cartMapper, const W& wells, const std::unordered_map<std::string, std::set<int>>& possibleFutureConnections, bool useWellConn, std::vector<std::set<int>>& wellGraph, int numJacobiBlocks)
     {
         if ( grid.comm().size() > 1 || numJacobiBlocks > 1)
         {
@@ -62,7 +62,7 @@ namespace detail
                     cart[ cartMapper.cartesianIndex( i ) ] = i;
 
                 Dune::cpgrid::WellConnections well_indices;
-                well_indices.init(wells, &possibleFutureConnections, cpgdim, cart);
+                well_indices.init(wells, possibleFutureConnections, cpgdim, cart);
 
                 for (auto& well : well_indices)
                 {

opm/simulators/wells/BlackoilWellModel_impl.hpp

@@ -194,11 +194,22 @@ namespace Opm {
 
         // Create cartesian to compressed mapping
         const auto& schedule_wells = this->schedule().getWellsatEnd();
+        const auto& possibleFutureConnections = this->schedule().getPossibleFutureConnections();
 
         // initialize the additional cell connections introduced by wells.
         for (const auto& well : schedule_wells)
         {
             std::vector<int> wellCells = this->getCellsForConnections(well);
+            // Now add the cells of the possible future connections
+            const auto possibleFutureConnectionSetIt = possibleFutureConnections.find(well.name());
+            if (possibleFutureConnectionSetIt != possibleFutureConnections.end()) {
+                for (auto& global_index : possibleFutureConnectionSetIt->second) {
+                    int compressed_idx = compressedIndexForInterior(global_index);
+                    if (compressed_idx >= 0) { // Ignore connections in inactive/remote cells.
+                        wellCells.push_back(compressed_idx);
+                    }
+                }
+            }
             for (int cellIdx : wellCells) {
                 neighbors[cellIdx].insert(wellCells.begin(),
                                           wellCells.end());

parallelTests.cmake

@@ -236,3 +236,15 @@ add_test_compareSeparateECLFiles(CASENAME actionx_compdat_8_procs
                                  REL_TOL ${rel_tol}
                                  IGNORE_EXTRA_KW BOTH
                                  MPI_PROCS 8)
+
+add_test_compareSeparateECLFiles(CASENAME actionx_compdat_nldd_1_proc
+                                 DIR1 actionx
+                                 FILENAME1 COMPDAT_SHORT
+                                 DIR2 actionx
+                                 FILENAME2 ACTIONX_COMPDAT_SHORT
+                                 SIMULATOR flow
+                                 ABS_TOL ${abs_tol}
+                                 REL_TOL ${rel_tol}
+                                 IGNORE_EXTRA_KW BOTH
+                                 MPI_PROCS 1
+                                 TEST_ARGS --nonlinear-solver=nldd --matrix-add-well-contributions=true)