New feature: multiply well edge weights for partitioning

CMakeLists_files.cmake

@@ -99,6 +99,7 @@ list(APPEND TEST_SOURCE_FILES
   tests/test_repairzcorn.cpp
   tests/test_sparsetable.cpp
   tests/test_subgridpart.cpp
+  tests/cpgrid/combinedgridwellgraph_test.cpp
   tests/cpgrid/distribution_test.cpp
   tests/cpgrid/entityrep_test.cpp
   tests/cpgrid/entity_test.cpp

opm/grid/GraphOfGrid.cpp

@@ -356,6 +356,37 @@ void GraphOfGrid<Grid>::addNeighboringCellsToWells ()
     }
 }
 
+template<typename Grid>
+void GraphOfGrid<Grid>::multiplyWellConnectivity (const std::set<int>& well, const WeightType& factor)
+{
+    // check neighbors of the vertex and
+    // multiply the edges that lead to other well connections
+    for (const auto& conn : well) {
+        if (graph.contains(conn)) {
+            auto& edges = graph[conn].edges;
+            for (auto& edge : edges) {
+                // process only higher IDs, do not search through well twice
+                if (edge.first < conn) {
+                    if (well.find(edge.first) != well.end()) {
+                        edge.second *= factor;
+                        assert(graph.contains(edge.first));
+                        auto& otherEdges = graph[edge.first].edges;
+                        assert(otherEdges.contains(conn));
+                        otherEdges[conn] *= factor;
+                    }
+                }
+            }
+        } else if (wellID(conn)==-1) {
+            OPM_THROW(std::domain_error, "Vertex is not present in the graph of grid.");
+        } else {
+            // conn got contracted into other vertex
+            // note: if conn is the smallest (=identifying) ID of a well,
+            //       graph contains it and everything works fine...
+            OPM_THROW(std::domain_error, "Mixing vertex contraction (addWell) and multiplyWellConnectivity is not supported.");
+        }
+    }
+}
+
 template class GraphOfGrid<Dune::CpGrid>;
 
 } // namespace Opm

opm/grid/GraphOfGrid.hpp

@@ -179,6 +179,9 @@ class GraphOfGrid{
         }
     }
 
+    /// \brief Multiply edges between well cells by a factor
+    void multiplyWellConnectivity (const std::set<int>& well, const WeightType& factor);
+
 private:
     /// \brief Create a graph representation of the grid
     ///

opm/grid/GraphOfGridWrappers.cpp

@@ -547,12 +547,16 @@ zoltanPartitioningWithGraphOfGrid(const Dune::CpGrid& grid,
     setDefaultZoltanParameters(zz);
     Zoltan_Set_Param(zz, "IMBALANCE_TOL", std::to_string(zoltanImbalanceTol).c_str());
     int layers = 0; // extra layers of cells attached to wells to distance them from boundary
+    float mIWC = 1; // multiply edge weights between a well's cells. Used when allowDistributedWells=true.
     for (const auto& [key, value] : params)
     {
-        if (key=="EnvelopeWellLayers")
+        if (key=="EnvelopeWellLayers") {
             layers = std::stoi(value);
-        else
+        } else if (key=="MultiplyWellConnectivities") {
+            mIWC = std::stof(value);
+        } else {
             Zoltan_Set_Param(zz, key.c_str(), value.c_str());
+        }
     }
 
     // root process has the whole grid, other ranks nothing
@@ -568,6 +572,11 @@ zoltanPartitioningWithGraphOfGrid(const Dune::CpGrid& grid,
         // skip cell contraction if wells can be distributed over multiple processes
         addWellConnections(gog, wellConnections);
         gog.addNeighboringCellsToWells(layers);
+    } else if (mIWC != 1) {
+        // multiply edge weights between connections of a well (not between two wells)
+        for (const auto& well : wellConnections) {
+            gog.multiplyWellConnectivity(well, mIWC);
+        }
     }
 
     // call partitioner

opm/grid/common/MetisPartition.cpp

@@ -276,7 +276,7 @@ metisSerialGraphPartitionGridOnRoot(const CpGrid& cpgrid,
         if( wells )
         {
             // well edge weight for partitioning, big enough that wells should not get split
-            idx_t weWeight = sumOfGridEdges<idx_t>(cpgrid, *gridAndWells);
+            idx_t weWeight = calculateWellEdgeWeight<idx_t>(cpgrid, *gridAndWells);
 
             int neighborCounter = 0;
             for( int cell = 0; cell < n;  cell++ )

opm/grid/common/ZoltanGraphFunctions.cpp

@@ -180,13 +180,29 @@ void getNullEdgeList(void *cpGridPointer, int sizeGID, int sizeLID,
 }
 
 template <typename EdgeWeightType>
-EdgeWeightType sumOfGridEdges(const Dune::CpGrid& grid,
-                              const CombinedGridWellGraph& graph)
+EdgeWeightType calculateWellEdgeWeight(const Dune::CpGrid& grid,
+                                       const CombinedGridWellGraph& graph)
 {
     double total = 0.0;
     for (int edge = 0; edge < grid.numFaces(); ++edge) {
         total += graph.edgeWeight(edge);
     }
+
+    // when multipltWellConnectivities is provided, set the well weight to the average of grid weight times that coefficient
+    double mWC = graph.getMultiplyWellConnectivities();
+    if (mWC >= 0) {
+        if (total != std::numeric_limits<double>::infinity()) {
+            total /= grid.numFaces();
+        } else {
+            // grid is too big, use maximum instead of the average
+            total = 0;
+            for (int edge=0; edge<grid.numFaces(); ++edge) {
+                total = std::max(total, graph.transmissibility(edge));
+            }
+        }
+        total *= mWC;
+    }
+
     const double maxVal = static_cast<double>(std::numeric_limits<EdgeWeightType>::max());
     if (total > maxVal) {
         return std::numeric_limits<EdgeWeightType>::max();
@@ -324,7 +340,7 @@ void getCpGridWellsEdgeList(void *graphPointer, int sizeGID, int sizeLID,
     int neighborCounter = 0;
 
     // well edge weight for partitioning, big enough that wells should not get split
-    float weWeight = sumOfGridEdges<float>(grid, graph);
+    float weWeight = calculateWellEdgeWeight<float>(grid, graph);
 
     for( int cell = 0; cell < numCells;  cell++ )
     {
@@ -436,14 +452,14 @@ void fillNBORGIDForSpecificCellAndIncrementNeighborCounter(const Dune::CpGrid&,
 template
 void fillNBORGIDAndWeightsForSpecificCellAndIncrementNeighborCounterForGridWithWells(const CombinedGridWellGraph&, const int, int*, int&, int*&, int*, const int&);
 template
-int sumOfGridEdges(const Dune::CpGrid& grid, const CombinedGridWellGraph& graph);
+int calculateWellEdgeWeight(const Dune::CpGrid& grid, const CombinedGridWellGraph& graph);
 
 template
 void fillNBORGIDForSpecificCellAndIncrementNeighborCounter(Dune::CpGrid const&, int, long*, int&, long*&);
 template
 void fillNBORGIDAndWeightsForSpecificCellAndIncrementNeighborCounterForGridWithWells(Dune::cpgrid::CombinedGridWellGraph const&, int, long*, int&, long*&, long*, const long&);
 template
-long sumOfGridEdges(const Dune::CpGrid& grid, const CombinedGridWellGraph& graph);
+long calculateWellEdgeWeight(const Dune::CpGrid& grid, const CombinedGridWellGraph& graph);
 
 #endif
 

opm/grid/common/ZoltanGraphFunctions.hpp

@@ -193,6 +193,17 @@ class CombinedGridWellGraph
         else
             return 1.0;
     }
+
+    void setMultiplyWellConnectivities(const double& mWC)
+    {
+        multiplyWellConnectivities = mWC;
+    }
+
+    double getMultiplyWellConnectivities() const
+    {
+        return multiplyWellConnectivities;
+    }
+
 private:
 
     void addCompletionSetToGraph()
@@ -239,17 +250,18 @@ class CombinedGridWellGraph
     int edgeWeightsMethod_;
     WellConnections well_indices_;
     double log_min_;
+    double multiplyWellConnectivities = -1;
 };
 
 /// \brief Get the number of edges of the graph of the grid and the wells for one cell
 int getNumberOfEdgesForSpecificCellForGridWithWells(const CombinedGridWellGraph& graph, int localCellId);
 
-/// \brief Iterate over the grid and get the sum of all edge weights
+/// \brief Iterate over the grid to calculate the edge weights between well connections
 ///
-/// Used as a weight for edges between cells of a well
+/// The result is the sum of all grid edges, or multiplyWellConnectivities*grid_average if the multiplier is positive
 template <typename EdgeWeightType>
-EdgeWeightType sumOfGridEdges(const Dune::CpGrid& grid,
-                              const CombinedGridWellGraph& graph);
+EdgeWeightType calculateWellEdgeWeight(const Dune::CpGrid& grid,
+                                       const CombinedGridWellGraph& graph);
 
 /// \brief Get the list of edges and weights for one cell of a grid with wells
 template<typename ID, typename weightType>

opm/grid/common/ZoltanPartition.cpp

@@ -326,8 +326,14 @@ zoltanGraphPartitionGridOnRoot(const CpGrid& cpgrid,
     }
     setDefaultZoltanParameters(zz);
     Zoltan_Set_Param(zz, "IMBALANCE_TOL", std::to_string(zoltanImbalanceTol).c_str());
-    for (const auto& [key, value] : params)
-        Zoltan_Set_Param(zz, key.c_str(), value.c_str());
+    double mWC = -1;
+    for (const auto& [key, value] : params) {
+        if (key=="MultiplyWellConnectivities") {
+            mWC = std::stod(value);
+        } else {
+            Zoltan_Set_Param(zz, key.c_str(), value.c_str());
+        }
+    }
 
     // For the load balancer one process has the whole grid and
     // all others an empty partition before loadbalancing.
@@ -344,6 +350,7 @@ zoltanGraphPartitionGridOnRoot(const CpGrid& cpgrid,
                                                        transmissibilities,
                                                        partitionIsEmpty,
                                                        edgeWeightsMethod));
+        gridAndWells->setMultiplyWellConnectivities(mWC);
         Dune::cpgrid::setCpGridZoltanGraphFunctions(zz, *gridAndWells,
                                                     partitionIsEmpty);
     }
@@ -523,8 +530,13 @@ class ZoltanSerialPartitioner
         else
             Zoltan_Set_Param(zz, "NUM_GLOBAL_PARTS", std::to_string(cc.size()).c_str());
 
-        for (const auto& [key, value] : params)
-            Zoltan_Set_Param(zz, key.c_str(), value.c_str());
+        for (const auto& [key, value] : params) {
+            if (key=="MultiplyWellConnectivities") {
+                gridAndWells->setMultiplyWellConnectivities(std::stod(value));
+            } else {
+                 Zoltan_Set_Param(zz, key.c_str(), value.c_str());
+            }
+        }
 
         // For the load balancer one process has the whole grid and
         // all others an empty partition before loadbalancing.