diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake
index 0d31c27e3..178484d85 100644
--- a/CMakeLists_files.cmake
+++ b/CMakeLists_files.cmake
@@ -51,6 +51,7 @@ list(APPEND MAIN_SOURCE_FILES
opm/grid/ColumnExtract.cpp
opm/grid/FaceQuadrature.cpp
opm/grid/GraphOfGrid.cpp
+ opm/grid/CoarseGraphOfGrid.cpp
opm/grid/GraphOfGridWrappers.cpp
opm/grid/GridHelpers.cpp
opm/grid/GridManager.cpp
@@ -91,6 +92,7 @@ list(APPEND TEST_SOURCE_FILES
tests/test_geom2d.cpp
tests/test_graphofgrid.cpp
tests/test_graphofgrid_parallel.cpp
+ tests/test_coarsegraphofgrid.cpp
tests/test_gridutilities.cpp
tests/test_minpvprocessor.cpp
tests/test_polyhedralgrid.cpp
@@ -251,6 +253,7 @@ list(APPEND PUBLIC_HEADER_FILES
opm/grid/ColumnExtract.hpp
opm/grid/FaceQuadrature.hpp
opm/grid/GraphOfGrid.hpp
+ opm/grid/CoarseGraphOfGrid.hpp
opm/grid/GraphOfGridWrappers.hpp
opm/grid/GridHelpers.hpp
opm/grid/GridManager.hpp
diff --git a/opm/grid/CoarseGraphOfGrid.cpp b/opm/grid/CoarseGraphOfGrid.cpp
new file mode 100644
index 000000000..91e34fc92
--- /dev/null
+++ b/opm/grid/CoarseGraphOfGrid.cpp
@@ -0,0 +1,291 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*
+ Copyright 2026 Equinor ASA.
+
+ This file is part of the Open Porous Media project (OPM).
+
+ OPM is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OPM is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OPM. If not, see .
+
+ Consult the COPYING file in the top-level source directory of this
+ module for the precise wording of the license and the list of
+ copyright holders.
+*/
+
+#include
+#include "CoarseGraphOfGrid.hpp"
+
+#include
+
+namespace Opm {
+
+template
+void CoarseGraphOfGrid::mergeWellCellsForCoarseGraph(std::vector& hasWell,
+ std::vector>& wellPerf,
+ const Dune::cpgrid::WellConnections& wellConn)
+{
+ int wellId = 0;
+
+ for (const auto& well : wellConn) {
+
+ bool cellInMultWells = false;
+ std::set otherWell;
+ std::vector perfs;
+
+ // Loop over all connections in the well. Create list perf of connections
+ for (int idx : well) {
+
+ if (!cellInMultWells) {
+ //Check if connection is intersected by other well
+ if (hasWell[idx]!=-1) {
+ cellInMultWells = true;
+ otherWell.insert( hasWell[idx] );
+ } else {
+ hasWell[idx] = wellId;
+ perfs.push_back(idx);
+ }
+ }
+ else {
+ if (hasWell[idx]!=-1) {
+ otherWell.insert( hasWell[idx] );
+ }
+ }
+ }
+ // If current well intersects with other wells. Add the well connections together
+ if (cellInMultWells) {
+ int minOtherWell = *otherWell.begin();
+ for (int idx : well) {
+ if (hasWell[idx]!=minOtherWell) {
+ hasWell[idx] = minOtherWell;
+ wellPerf[minOtherWell].push_back(idx);
+ }
+ }
+ // If more than one intersection add all wells to smallest wellId and delete others.
+ for (int otherWellId : otherWell) {
+ if (otherWellId != minOtherWell) {
+ for (int idx : wellPerf[otherWellId]) {
+ if (hasWell[idx]!=minOtherWell) {
+ hasWell[idx] = minOtherWell;
+ wellPerf[minOtherWell].push_back(idx);
+ }
+ }
+ wellPerf[otherWellId].clear();
+ }
+ }
+ }
+ else {
+ wellPerf.push_back(perfs);
+ wellId++;
+ }
+ }
+}
+
+template
+void CoarseGraphOfGrid::dfsqw(const Row& row, std::priority_queue &q, int v, int master,
+ double w, int maxNode, std::vector& visited,
+ std::vector>>& gEdges,
+ const std::vector& hasWell, const std::vector>& wellPerf)
+{
+
+ auto& current_cnode = coarseNodes.back();
+ auto& current_edges = gEdges.back();
+
+ std::vector wellIdxs;
+ if (hasWell[v] == -1) {
+ visited[v] = true;
+ map_to_coarse_[v] = master;
+ current_cnode.push_back(v);
+
+ // Add all neighboring vertices of v with transmissibility larger than w to the queue.
+ auto col = row.begin();
+ for (; col != row.end(); ++col) {
+ int nab = col.index();
+ double wgt = (*transGraph)[v][nab];
+ if ( wgt > w) {
+ if (!visited[nab]) {
+ q.push({wgt, nab});
+ }
+ }
+ }
+ } else {
+ // If node has a well, merge all connections to master node.
+ int wellId = hasWell[v];
+ const std::vector& perfs = wellPerf[wellId];
+
+ for (const auto& idx : perfs) {
+ visited[idx] = true;
+ map_to_coarse_[idx] = master;
+ current_cnode.push_back(idx);
+ wellIdxs.push_back(idx);
+ }
+ for (const auto& idx : perfs) {
+ auto wrow = (*transGraph)[idx];
+ auto col = wrow.begin();
+ for (; col != wrow.end(); ++col) {
+ int nab = col.index();
+ double wgt = (*transGraph)[idx][nab];
+ if ( wgt > w) {
+ if (!visited[nab]) {
+ q.push({wgt, nab});
+ }
+ }
+ }
+ }
+ }
+
+ // Only merge more vertices if the current coarse node is smaller than maxNode.
+ if ( (int)current_cnode.size() < maxNode ) {
+ if (!q.empty()) {
+
+ // Find strongest connection in queue q not already merged to current_cnode.
+ auto strongCon = q.top();
+ int nab = strongCon.idx;
+ q.pop();
+ while (visited[nab] && !q.empty()) {
+ strongCon = q.top();
+ nab = strongCon.idx;
+ q.pop();
+ }
+ // Call dfsq reflexively on strongest connection in queue
+ if (!visited[nab])
+ dfsqw((*transGraph)[nab],q,nab,master,w,maxNode,visited,gEdges,hasWell,wellPerf);
+ }
+ } else {
+ q = std::priority_queue();
+ }
+
+ // Add connection between v and nab if v and nab are not merged
+ if (wellIdxs.size() > 0) {
+ for (const auto& idx : wellIdxs) {
+ auto wrow = (*transGraph)[idx];
+ auto col = wrow.begin();
+ for (; col != wrow.end(); ++col) {
+ int nab = col.index();
+ if (map_to_coarse_[v]!=map_to_coarse_[nab]) {
+ current_edges.push_back({v,nab,(*transGraph)[idx][nab]});
+ }
+ }
+ }
+ }
+ else {
+ auto col = row.begin();
+ for (; col != row.end(); ++col) {
+ int nab = col.index();
+ if (map_to_coarse_[v]!=map_to_coarse_[nab]) {
+ current_edges.push_back({v,nab,(*transGraph)[v][nab]});
+ }
+ }
+ }
+}
+
+template
+void CoarseGraphOfGrid::createCoarseGraph(const Dune::EdgeWeightMethod edgeWeightMethod,
+ double coarseThreshold,
+ int coarsePartitionMaxNodeSize,
+ bool allowDistributedWells,
+ int root,
+ const Dune::cpgrid::WellConnections& wellConn)
+{
+ int N = grid.size(0);
+ const auto& rank = grid.comm().rank();
+
+ // List to keep track of visited vertices in original fine graph
+ std::vector visited(N, false);
+ // List to know if vertex i is perferated by a well. hasWell[i]==-1 means no, hasWell[i]!=-1 yes
+ std::vector hasWell(N,-1);
+ // For each well i, wellPerf[i] lists vertices perferated by well i.
+ // If wells intersect they are merged.
+ std::vector> wellPerf;
+
+ // Only add data to hasWell and wellPerf if allowDistributedWells==false
+ if (!allowDistributedWells) {
+ mergeWellCellsForCoarseGraph(hasWell, wellPerf, wellConn);
+ }
+ // Map from index of fine graph to index of coarse graph.
+ map_to_coarse_.resize(N, -1);
+
+ // Vector that describes the coarse graph.
+ // Each coarse vertex v has vector<(int,int,double)> = gEdges[v],
+ // where each tuple (idx,nab,wgt) is index idx of fine graph vertex idx contained in v,
+ // (idx,nab) is edge in fine graph, and wgt is transmissibility between idx and nab.
+ std::vector >> gEdges;
+
+ // Counter for coarse graph index
+ int newV = 0;
+ // Keep track of largest coarse vertex.
+ int biggest = 0;
+
+ if (rank == root) {
+ // Loop over all idecies in fine graph
+ for (int v = 0; v < N; ++v) {
+
+ // if idx v is not visited, add it to coarse graph
+ if (!visited[v]) {
+
+ std::priority_queue q;
+
+ // Allocate coarse node v in gEdges and coarseNodes
+ gEdges.emplace_back();
+ coarseNodes.emplace_back();
+
+ // Call depth first search from row transGraph[v]
+ dfsqw((*transGraph)[v],q,v,newV,coarseThreshold,
+ coarsePartitionMaxNodeSize,visited,
+ gEdges,hasWell,wellPerf);
+
+ newV++;
+
+ if ((int)coarseNodes.back().size() > biggest)
+ biggest = coarseNodes.back().size();
+ }
+ }
+ std::cout << "Coarse partitioning graph size: " << coarseNodes.size() <<" Largest node: "<< biggest << std::endl;
+
+ // Construct the coarse graph from gEdges
+ //
+ // Loop over coarse nodes with cIdx=0,...,newV -1
+ for (const std::vector >& es : gEdges ) {
+
+ // ce represents the edges of cIdx. If cIdx has connection with cNab,
+ // wgt=ce[cNab] exists and is greater than 0.
+ std::map ce;
+
+ // Loop over edges of all vertices from fine graph contained in coarse node cIdx
+ for (const std::tuple& fe : es) {
+
+ // Coarse index of the fine graph edge
+ int coarseNab = map_to_coarse_[std::get<1>(fe)];
+ // Transmissibility of fine edge
+ double transVal = std::get<2>(fe);
+ // Besed on edgeWeightMethod, choose weight of coarse edgeWeight
+ double weight = edgeWeightMethod == 0 ? 1.0 : transVal;
+
+ // Only add fine edge to coarse graph if transmissibility is non-zero. Overlap cells
+ // between partitions are not added if connection has non-zero transmissibility.
+ if (transVal > 0) {
+ if ( ce.count(coarseNab) == 1 ) {
+ ce[coarseNab] += weight;
+ } else {
+ ce.insert({coarseNab,weight});
+ }
+ }
+ }
+ cedges.push_back(ce);
+ }
+ }
+}
+
+template class CoarseGraphOfGrid;
+
+} // namespace Opm
diff --git a/opm/grid/CoarseGraphOfGrid.hpp b/opm/grid/CoarseGraphOfGrid.hpp
new file mode 100644
index 000000000..5db2321e8
--- /dev/null
+++ b/opm/grid/CoarseGraphOfGrid.hpp
@@ -0,0 +1,148 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*
+ Copyright 2026 Equinor ASA.
+
+ This file is part of the Open Porous Media project (OPM).
+
+ OPM is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OPM is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OPM. If not, see .
+
+ Consult the COPYING file in the top-level source directory of this
+ module for the precise wording of the license and the list of
+ copyright holders.
+*/
+
+#ifndef OPM_COARSE_GRAPH_OF_GRID_HEADER
+#define OPM_COARSE_GRAPH_OF_GRID_HEADER
+
+#include
+#include
+#include
+
+namespace Opm {
+
+struct WgtIdx {
+
+ double wgt;
+ int idx;
+
+ bool operator<(const WgtIdx& other) const {
+ return wgt < other.wgt;
+ }
+};
+
+/// \brief A class storing a Coarse graph representation of the grid
+///
+/// Similar to GraphOfGrid, but here nodes are merged if
+/// the transmissibility on the edge connecting them is below
+/// a certain threshold.
+template
+class CoarseGraphOfGrid{
+
+ using TransGraph = Dune::BCRSMatrix>;
+ using Row = typename TransGraph::row_type;
+
+public:
+
+ explicit CoarseGraphOfGrid (const Grid& grid_,
+ const Dune::EdgeWeightMethod edgeWeightMethod,
+ const TransGraph* tg,
+ double coarseThreshold,
+ int coarsePartitionMaxNodeSize,
+ bool allowDistributedWells,
+ int root,
+ const Dune::cpgrid::WellConnections& wellConn)
+ : grid(grid_), transGraph(tg)
+ {
+ createCoarseGraph(edgeWeightMethod, coarseThreshold, coarsePartitionMaxNodeSize,
+ allowDistributedWells, root, wellConn);
+ }
+
+ const Grid& getGrid() const
+ {
+ return grid;
+ }
+
+ /// \brief Number of graph vertices
+ int size() const
+ {
+ return coarseNodes.size();
+ }
+
+ const std::vector>& getCoarseNodes() const
+ {
+ return coarseNodes;
+ }
+
+ const std::vector >& getCoarseEdges() const
+ {
+ return cedges;
+ }
+
+ /// \brief returns map of global cell id to coarse vertex id
+ const std::vector& getMapToCoarse() const
+ {
+ return map_to_coarse_;
+ }
+
+ /// \brief Return the list of wells
+ const auto& getWells () const
+ {
+ return wells;
+ }
+private:
+
+ /// \brief Merge vertices that share a common well
+ void mergeWellCellsForCoarseGraph(std::vector& hasWell,
+ std::vector>& wellPerf,
+ const Dune::cpgrid::WellConnections& wells);
+
+ /// \brief Coarsen the graph by doing a (d)epth (f)irst (s)earch with priority (q)ueue
+ /// and merges (w)ells.
+ ///
+ /// Recurseive dfs that merges strongly connected nodes (transmissibility) in the
+ /// partitioning graph.
+ /// A priority queue is used to make sure the largest connections are prioritised in the
+ /// merging of vertices.
+ void dfsqw(const Row& row, std::priority_queue &q, int v, int master,
+ double w, int maxNode, std::vector& visited,
+ std::vector>>& gEdges,
+ const std::vector& hasWell, const std::vector>& wellPerf);
+
+ /// \brief Create the coarse graph merging all vertices connected with a large transmissibility.
+ ///
+ /// Creates a coarsened graph by merging vertices connected with a
+ /// transmissibility larger than coarseThreshold.
+ /// Coarse vertices will not be larger than coarsePartitionMaxNodeSize.
+ /// Also meges wells if allowDistributedWells == true.
+ void createCoarseGraph(const Dune::EdgeWeightMethod edgeWeightMethod,
+ double coarseThreshold,
+ int coarsePartitionMaxNodeSize,
+ bool allowDistributedWells,
+ int root,
+ const Dune::cpgrid::WellConnections& wells);
+
+ const Grid& grid;
+ std::list> wells;
+
+ const TransGraph* transGraph;
+ std::vector map_to_coarse_;
+ std::vector > cedges;
+ std::vector> coarseNodes;
+
+};
+
+} // namespace Opm
+
+#endif // OPM_COARSE_GRAPH_OF_GRID_HEADER
diff --git a/opm/grid/CpGrid.hpp b/opm/grid/CpGrid.hpp
index 7cbb06627..8e44929e4 100644
--- a/opm/grid/CpGrid.hpp
+++ b/opm/grid/CpGrid.hpp
@@ -43,6 +43,7 @@
#include
#include
+#include
#include
@@ -958,12 +959,16 @@ namespace Dune
bool addCornerCells=false, int overlapLayers=1, int partitionMethod = Dune::PartitionMethod::zoltanGoG,
double imbalanceTol = 1.1,
bool allowDistributedWells = false,
- bool useTransToFilterOverlap = true)
+ bool useTransToFilterOverlap = true,
+ const Dune::BCRSMatrix>* graph = nullptr,
+ double coarseThreshold = 1.0,
+ int coarsePartitionMaxNodeSize = -1)
{
auto ret = scatterGrid(method, ownersFirst, wells, possibleFutureConnections, serialPartitioning, transmissibilities,
addCornerCells, overlapLayers, partitionMethod, imbalanceTol, allowDistributedWells,
/* input_cell_parts = */ std::vector{}, /* level = */ 0,
- useTransToFilterOverlap);
+ useTransToFilterOverlap, graph,
+ coarseThreshold, coarsePartitionMaxNodeSize);
using std::get;
if (get<0>(ret))
{
@@ -1481,7 +1486,10 @@ namespace Dune
bool allowDistributedWells = true,
const std::vector& input_cell_part = {},
int level = -1,
- bool useTransToFilterOverlap = true);
+ bool useTransToFilterOverlap = true,
+ const Dune::BCRSMatrix>* transGraph = nullptr,
+ double coarseThreshold = 1.0,
+ int coarsePartitionMaxNodeSize = -1);
/** @brief The data stored in the grid.
*
diff --git a/opm/grid/GraphOfGridWrappers.cpp b/opm/grid/GraphOfGridWrappers.cpp
index 34e540df6..7d8caca77 100644
--- a/opm/grid/GraphOfGridWrappers.cpp
+++ b/opm/grid/GraphOfGridWrappers.cpp
@@ -32,9 +32,10 @@
namespace Opm {
#if HAVE_MPI
+template
int getGraphOfGridNumVertices(void* pGraph, int *err)
{
- const GraphOfGrid& gog = *static_cast*>(pGraph);
+ const GraphType& gog = *static_cast(pGraph);
int size = gog.size();
*err = ZOLTAN_OK;
return size;
@@ -131,12 +132,90 @@ void getGraphOfGridEdgeList(void *pGraph,
*err = ZOLTAN_OK;
}
-template
+void getCoarseGraphVerticesList(void* pGraph,
+ [[maybe_unused]] int dimGlobalID,
+ [[maybe_unused]] int dimLocalID,
+ ZOLTAN_ID_PTR gIDs,
+ [[maybe_unused]] ZOLTAN_ID_PTR lIDs,
+ int weightDim,
+ float *objWeights,
+ int *err)
+{
+ assert(dimGlobalID==1); // ID is a single int
+ assert(weightDim==1); // vertex weight is a single float
+ const CoarseGraphOfGrid& gog = *static_cast*>(pGraph);
+ const std::vector>& cnodes = gog.getCoarseNodes();
+ int i=0;
+ for (const auto& v : cnodes)
+ {
+ gIDs[i] = i;
+ // lIDs are left unused
+ objWeights[i] = v.size();
+ ++i;
+ }
+ *err = ZOLTAN_OK;
+}
+void getCoarseGraphNumEdges(void *pGraph,
+ [[maybe_unused]] int dimGlobalID,
+ [[maybe_unused]] int dimLocalID,
+ [[maybe_unused]] int numCells,
+ ZOLTAN_ID_PTR gIDs,
+ [[maybe_unused]] ZOLTAN_ID_PTR lIDs,
+ int *numEdges,
+ int *err)
+{
+ assert(dimGlobalID==1); // ID is a single int
+ const CoarseGraphOfGrid& gog = *static_cast*>(pGraph);
+
+ const std::vector >& edges = gog.getCoarseEdges();
+
+ for (size_t idx = 0; idx < edges.size(); ++idx)
+ {
+ numEdges[idx] = edges[gIDs[idx]].size();
+ }
+
+ *err = ZOLTAN_OK;
+}
+
+void getCoarseGraphEdgeList(void *pGraph,
+ [[maybe_unused]] int dimGlobalID,
+ [[maybe_unused]] int dimLocalID,
+ [[maybe_unused]] int numCells,
+ [[maybe_unused]] ZOLTAN_ID_PTR gIDs,
+ [[maybe_unused]] ZOLTAN_ID_PTR lIDs,
+ [[maybe_unused]] int *numEdges,
+ ZOLTAN_ID_PTR nborGIDs,
+ int *nborProc,
+ int weightDim,
+ float *edgeWeights,
+ int *err)
+{
+ assert(dimGlobalID==1); // ID is a single int
+ assert(weightDim==1); // edge weight is a single float
+ const CoarseGraphOfGrid& gog = *static_cast*>(pGraph);
+ const std::vector >& edges = gog.getCoarseEdges();
+ int id=0;
+
+ const auto& rank = gog.getGrid().comm().rank();
+ for (const auto& node : edges)
+ {
+ for (const auto& edge : node) {
+ nborGIDs[id] = edge.first;
+ edgeWeights[id] = edge.second;
+ nborProc[id++] = rank;
+ }
+ }
+
+ *err = ZOLTAN_OK;
+}
+
+template
void setGraphOfGridZoltanGraphFunctions(Zoltan_Struct *zz,
- GraphOfGrid& gog,
+ GraphType& gog,
bool pretendNull)
{
- GraphOfGrid* pGraph = &gog;
+ using DecayedGraph = std::decay_t;
+ GraphType* pGraph = &gog;
if (pretendNull)
{
Zoltan_Set_Num_Obj_Fn(zz, Dune::cpgrid::getNullNumCells, pGraph);
@@ -146,12 +225,21 @@ void setGraphOfGridZoltanGraphFunctions(Zoltan_Struct *zz,
}
else
{
- Zoltan_Set_Num_Obj_Fn(zz, getGraphOfGridNumVertices, pGraph);
- Zoltan_Set_Obj_List_Fn(zz, getGraphOfGridVerticesList, pGraph);
- Zoltan_Set_Num_Edges_Multi_Fn(zz, getGraphOfGridNumEdges, pGraph);
- Zoltan_Set_Edge_List_Multi_Fn(zz, getGraphOfGridEdgeList, pGraph);
+ if constexpr (std::is_same_v>) {
+ Zoltan_Set_Num_Obj_Fn(zz, getGraphOfGridNumVertices, pGraph);
+ Zoltan_Set_Obj_List_Fn(zz, getGraphOfGridVerticesList, pGraph);
+ Zoltan_Set_Num_Edges_Multi_Fn(zz, getGraphOfGridNumEdges, pGraph);
+ Zoltan_Set_Edge_List_Multi_Fn(zz, getGraphOfGridEdgeList, pGraph);
+ }
+ else if constexpr (std::is_same_v>) {
+ Zoltan_Set_Num_Obj_Fn(zz, getGraphOfGridNumVertices, pGraph);
+ Zoltan_Set_Obj_List_Fn(zz, getCoarseGraphVerticesList, pGraph);
+ Zoltan_Set_Num_Edges_Multi_Fn(zz, getCoarseGraphNumEdges, pGraph);
+ Zoltan_Set_Edge_List_Multi_Fn(zz, getCoarseGraphEdgeList, pGraph);
+ }
}
}
+
#endif // HAVE_MPI
void addFutureConnectionWells(GraphOfGrid& gog,
@@ -209,9 +297,9 @@ void extendGIDtoRank(const GraphOfGrid& gog,
#if HAVE_MPI
namespace Impl{
-
+template
std::vector>
-extendRootExportList(const GraphOfGrid& gog,
+extendRootExportList(const GOG& gog,
std::vector>& exportList,
int root,
const std::vector& gIDtoRank)
@@ -491,6 +579,128 @@ makeImportAndExportLists(const GraphOfGrid& gog,
std::move(myImportList) );
}
+template
+std::tuple,
+ std::vector>,
+ std::vector >,
+ std::vector > >
+makeImportAndExportLists(const CoarseGraphOfGrid& gog,
+ const Dune::Communication& cc,
+ const std::vector * wells,
+ const Dune::cpgrid::WellConnections& wellConnections,
+ const std::unordered_map>& possibleFutureConnections,
+ int root,
+ int numExport,
+ int numImport,
+ [[maybe_unused]] const Id* exportLocalGids,
+ const Id* exportGlobalGids,
+ const int* exportToPart,
+ [[maybe_unused]] const Id* importGlobalGids,
+ bool allowDistributedWells)
+{
+ int size = gog.getMapToCoarse().size();
+ cc.broadcast(&size, 1, root);
+ int rank = cc.rank();
+ std::vector gIDtoRank(size, rank);
+ std::vector > wellsOnProc;
+
+ // List entry: process to export to, (global) index, process rank, attribute there (not needed?)
+ std::vector> myExportList;
+ // List entry: process to import from, global index, process rank, attribute here, local index (determined later)
+ std::vector> myImportList;
+ float buffer = 1.05; // to allocate extra space for wells in myExportList and myImportList
+ assert(rank==root || numExport==0);
+ assert(rank!=root || numImport==0);
+ // all cells on root are added to its export and its import list
+ std::size_t reserveEx = rank!=root ? 0 : size;
+ std::size_t reserveIm = size*buffer/cc.size();
+ myExportList.reserve(reserveEx);
+ myImportList.reserve(reserveIm);
+ using AttributeSet = Dune::cpgrid::CpGridData::AttributeSet;
+
+ std::vector> importListFromRoot(cc.size());
+ std::vector sizeOfImport(cc.size(), 0);
+ if (rank==root)
+ {
+ std::vector coarsePartRes(gog.size(), root);
+ auto cnodes = gog.getCoarseNodes();
+ for ( int i=0; i < numExport; ++i )
+ {
+ coarsePartRes[exportGlobalGids[i]] = exportToPart[i];
+ sizeOfImport[exportToPart[i]] += cnodes[exportGlobalGids[i]].size();
+ }
+
+ const std::vector m2c = gog.getMapToCoarse();
+ for (int i = 0; i < size; ++i) {
+
+ gIDtoRank[i] = coarsePartRes[m2c[i]];
+ myExportList.emplace_back(i, coarsePartRes[m2c[i]], static_cast(AttributeSet::owner));
+ }
+
+ for ( std::size_t i = 0; i < gIDtoRank.size(); ++i)
+ {
+ if ( gIDtoRank[i] == rank )
+ {
+ myImportList.emplace_back(i, rank, static_cast(AttributeSet::owner), -1 );
+ }
+ else {
+ importListFromRoot[gIDtoRank[i]].emplace_back(i);
+ }
+ }
+ }
+ std::vector newImportList;
+ int newNumImport;
+ if (cc.rank() == root) {
+ std::vector requestSize(2 * (cc.size() - 1));
+
+ for (int i = 0; i < cc.size() - 1; ++i) {
+ int ii = i + (int)(i >= root); // ii takes values {0,...,mpisize-1} but skips root
+ int tag = 15; // a random number
+ MPI_Isend(&sizeOfImport[ii], 1, MPI_INT, ii, tag, cc, &requestSize[2 * i]);
+ MPI_Isend(importListFromRoot[ii].data(), sizeOfImport[ii], MPI_INT,ii, tag + 1, cc, &requestSize[2 * i + 1]);
+ }
+ newNumImport = 0;
+ MPI_Waitall(requestSize.size(), requestSize.data(), MPI_STATUS_IGNORE);
+ } else {
+ int tag = 15; // a random number
+ MPI_Recv(&newNumImport, 1, MPI_INT, root, tag, cc, MPI_STATUS_IGNORE);
+ newImportList.resize(newNumImport);
+ MPI_Recv(newImportList.data(), newNumImport, MPI_INT, root, tag + 1, cc, MPI_STATUS_IGNORE);
+ }
+
+ for ( int i=0; i < newNumImport; ++i )
+ {
+ myImportList.emplace_back(newImportList[i], root, static_cast(AttributeSet::owner), -1);
+ }
+ std::vector> parallel_wells;
+ if( wells )
+ {
+ if (allowDistributedWells) {
+ wellsOnProc = Dune::cpgrid::perforatingWellIndicesOnProc(gIDtoRank, *wells,
+ possibleFutureConnections,
+ gog.getGrid());
+ parallel_wells = Dune::cpgrid::computeParallelWells(wellsOnProc,
+ *wells,
+ cc,
+ root);
+ }
+ else {
+ auto wellRanks = getWellRanks(gIDtoRank, wellConnections);
+ parallel_wells = wellsOnThisRank(*wells, wellRanks, cc, root);
+ }
+ }
+ else
+ {
+ std::ranges::sort(myExportList);
+ std::ranges::sort(myImportList);
+ }
+ return std::make_tuple( std::move(gIDtoRank),
+ std::move(parallel_wells),
+ std::move(myExportList),
+ std::move(myImportList) );
+
+}
+
namespace {
void setDefaultZoltanParameters(Zoltan_Struct* zz)
{
@@ -654,6 +864,115 @@ zoltanPartitioningWithGraphOfGrid(const Dune::CpGrid& grid,
return importExportLists;
}
+std::tuple, std::vector>,
+ std::vector >,
+ std::vector >,
+ Dune::cpgrid::WellConnections>
+zoltanPartitioningWithCoarseGraph(const Dune::CpGrid& grid,
+ const std::vector * wells,
+ const std::unordered_map>& possibleFutureConnections,
+ const Dune::cpgrid::CpGridDataTraits::Communication& cc,
+ Dune::EdgeWeightMethod edgeWeightMethod,
+ int root,
+ const double zoltanImbalanceTol,
+ bool allowDistributedWells,
+ const std::map& params,
+ const Dune::BCRSMatrix>* transGraph,
+ double coarseThreshold,
+ int coarsePartitionMaxNodeSize)
+{
+ float ver = 0;
+ struct Zoltan_Struct *zz;
+ int changes, numGidEntries, numLidEntries, numImport, numExport;
+ ZOLTAN_ID_PTR importGlobalGids, importLocalGids, exportGlobalGids, exportLocalGids;
+ int *importProcs, *importToPart, *exportProcs, *exportToPart;
+ int argc=0;
+ char** argv = 0 ;
+ int rc = Zoltan_Initialize(argc, argv, &ver);
+ if (rc != ZOLTAN_OK) {
+ OPM_THROW(std::runtime_error, "Could not initialize Zoltan!");
+ }
+ zz = Zoltan_Create(cc);
+ if (zz == nullptr) {
+ OPM_THROW(std::runtime_error, "Could not create Zoltan data structures!");
+ }
+ 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());
+ }
+
+ // root process has the whole grid, other ranks nothing
+ bool partitionIsEmpty = cc.rank()!=root;
+
+ auto wellConnections = partitionIsEmpty || !wells ? Dune::cpgrid::WellConnections()
+ : Dune::cpgrid::WellConnections(*wells, possibleFutureConnections, grid);
+
+ // prepare graph and contract well cells
+ // non-root processes have empty grid and no wells
+ CoarseGraphOfGrid cgog(grid, edgeWeightMethod, transGraph,
+ coarseThreshold, coarsePartitionMaxNodeSize,
+ allowDistributedWells, root, wellConnections);
+
+ assert(cgog.size()==0 || !partitionIsEmpty);
+
+ // call partitioner
+ setGraphOfGridZoltanGraphFunctions(zz, cgog, partitionIsEmpty);
+ rc = Zoltan_LB_Partition(zz, /* input (all remaining fields are output) */
+ &changes, /* 1 if partitioning was changed, 0 otherwise */
+ &numGidEntries, /* Number of integers used for a global ID */
+ &numLidEntries, /* Number of integers used for a local ID */
+ &numImport, /* Number of vertices to be sent to me */
+ &importGlobalGids, /* Global IDs of vertices to be sent to me */
+ &importLocalGids, /* Local IDs of vertices to be sent to me */
+ &importProcs, /* Process rank for source of each incoming vertex */
+ &importToPart, /* New partition for each incoming vertex */
+ &numExport, /* Number of vertices I must send to other processes*/
+ &exportGlobalGids, /* Global IDs of the vertices I must send */
+ &exportLocalGids, /* Local IDs of the vertices I must send */
+ &exportProcs, /* Process to which I send each of the vertices */
+ &exportToPart); /* Partition to which each vertex will belong */
+ if (rc == ZOLTAN_WARN) {
+ OpmLog::warning("Zoltan_LB_Partition returned with warning");
+ } else if (rc == ZOLTAN_MEMERR) {
+ OPM_THROW(std::runtime_error, "Memory allocation failure in Zoltan_LB_Partition");
+ } else if (rc == ZOLTAN_FATAL) {
+ OPM_THROW(std::runtime_error, "Error returned from Zoltan_LB_Partition");
+ }
+
+ // arrange output into tuples and add well cells
+ auto prepareIELists = [&]() {
+
+ auto partResult = makeImportAndExportLists(cgog,
+ cc,
+ wells,
+ wellConnections,
+ possibleFutureConnections,
+ root,
+ numExport,
+ numImport,
+ exportLocalGids,
+ exportGlobalGids,
+ exportProcs,
+ importGlobalGids,
+ allowDistributedWells);
+ return std::tuple(std::move(std::get<0>(partResult)),
+ std::move(std::get<1>(partResult)),
+ std::move(std::get<2>(partResult)),
+ std::move(std::get<3>(partResult)),
+ std::move(wellConnections));
+ };
+ auto importExportLists = prepareIELists();
+
+ Zoltan_LB_Free_Part(&exportGlobalGids, &exportLocalGids, &exportProcs, &exportToPart);
+ Zoltan_LB_Free_Part(&importGlobalGids, &importLocalGids, &importProcs, &importToPart);
+ Zoltan_Destroy(&zz);
+
+ return importExportLists;
+}
+
std::vector >
makeExportListsFromGIDtoRank(const std::vector& gIDtoRank, int ccsize)
{
@@ -750,6 +1069,88 @@ applySerialZoltan (const Dune::CpGrid& grid,
Zoltan_Destroy(&zz);
return std::make_tuple(rc, gIDtoRank);
}
+
+std::tuple>
+applySerialZoltanCG (const Dune::CpGrid& grid,
+ const Dune::cpgrid::WellConnections& wellConnections,
+ int numParts,
+ Dune::EdgeWeightMethod edgeWeightMethod,
+ int root,
+ const double zoltanImbalanceTol,
+ bool allowDistributedWells,
+ const std::map& params,
+ const Dune::BCRSMatrix>* transGraph,
+ double coarseThreshold,
+ int coarsePartitionMaxNodeSize)
+{
+ int rc = ZOLTAN_OK;
+ ZOLTAN_ID_PTR importGlobalGids, importLocalGids, exportGlobalGids, exportLocalGids;
+ int numExport = 0, numImport = 0;
+ int *importProcs, *importToPart, *exportProcs, *exportToPart;
+ struct Zoltan_Struct* zz;
+ int changes, numGidEntries, numLidEntries;
+
+ int argc = 0;
+ char** argv = 0;
+ float ver = 0;
+ std::vector gIDtoRank;
+
+ rc = Zoltan_Initialize(argc, argv, &ver);
+ if (rc != ZOLTAN_OK)
+ return std::make_tuple(ZOLTAN_OK + 1, gIDtoRank);
+ zz = Zoltan_Create(MPI_COMM_SELF);
+ if (!zz)
+ return std::make_tuple(ZOLTAN_OK + 2, gIDtoRank);
+ setDefaultZoltanParameters(zz);
+ Zoltan_Set_Param(zz, "IMBALANCE_TOL", std::to_string(zoltanImbalanceTol).c_str());
+ Zoltan_Set_Param(zz, "NUM_GLOBAL_PARTS", std::to_string(numParts).c_str());
+
+ for (const auto& [key, value] : params) {
+ Zoltan_Set_Param(zz, key.c_str(), value.c_str());
+ }
+
+ // prepare graph and contract well cells
+ CoarseGraphOfGrid cgog(grid, edgeWeightMethod, transGraph,
+ coarseThreshold, coarsePartitionMaxNodeSize,
+ allowDistributedWells, root, wellConnections);
+
+
+ // call partitioner
+ setGraphOfGridZoltanGraphFunctions(zz, cgog, false);
+ rc = Zoltan_LB_Partition(zz, /* input (all remaining fields are output) */
+ &changes, /* 1 if partitioning was changed, 0 otherwise */
+ &numGidEntries, /* Number of integers used for a global ID */
+ &numLidEntries, /* Number of integers used for a local ID */
+ &numImport, /* Number of vertices to be sent to me */
+ &importGlobalGids, /* Global IDs of vertices to be sent to me */
+ &importLocalGids, /* Local IDs of vertices to be sent to me */
+ &importProcs, /* Process rank for source of each incoming vertex */
+ &importToPart, /* New partition for each incoming vertex */
+ &numExport, /* Number of vertices I must send to other processes*/
+ &exportGlobalGids, /* Global IDs of the vertices I must send */
+ &exportLocalGids, /* Local IDs of the vertices I must send */
+ &exportProcs, /* Process to which I send each of the vertices */
+ &exportToPart); /* Partition to which each vertex will belong */
+ numImport = 0;
+ if (rc == ZOLTAN_OK) {
+ gIDtoRank.resize(grid.numCells(), root);
+ std::vector coarsePartRes(cgog.size(), root);
+ const std::vector m2c = cgog.getMapToCoarse();
+ for (int i = 0; i < numExport; ++i) {
+ coarsePartRes[exportGlobalGids[i]] = exportToPart[i];
+ }
+ for (int i = 0; i < grid.numCells(); ++i) {
+ gIDtoRank[i] = coarsePartRes[m2c[i]];
+ }
+
+ } else {
+ rc = ZOLTAN_OK + 3; // distinguish Zoltan failures
+ }
+ Zoltan_LB_Free_Part(&exportGlobalGids, &exportLocalGids, &exportProcs, &exportToPart);
+ Zoltan_LB_Free_Part(&importGlobalGids, &importLocalGids, &importProcs, &importToPart);
+ Zoltan_Destroy(&zz);
+ return std::make_tuple(rc, gIDtoRank);
+}
} // end anonymous namespace
std::tuple,
@@ -849,6 +1250,108 @@ zoltanSerialPartitioningWithGraphOfGrid(const Dune::CpGrid& grid,
std::move(myImportList),
std::move(wellConnections));
}
+
+std::tuple,
+ std::vector>,
+ std::vector >,
+ std::vector >,
+ Dune::cpgrid::WellConnections>
+zoltanSerialPartitioningWithCoarseGraph(const Dune::CpGrid& grid,
+ const std::vector * wells,
+ const std::unordered_map>& possibleFutureConnections,
+ const Dune::cpgrid::CpGridDataTraits::Communication& cc,
+ Dune::EdgeWeightMethod edgeWeightMethod,
+ int root,
+ const double zoltanImbalanceTol,
+ bool allowDistributedWells,
+ const std::map& params,
+ const Dune::BCRSMatrix>* transGraph,
+ double coarseThreshold,
+ int coarsePartitionMaxNodeSize)
+{
+ // root process has the whole grid, other ranks nothing
+ bool partitionIsEmpty = cc.rank() != root;
+ int rc = ZOLTAN_OK;
+ std::vector gIDtoRank;
+ using AttributeSet = Dune::cpgrid::CpGridData::AttributeSet;
+ std::vector> myExportList;
+ std::vector> myImportList;
+ std::vector> exportedCells;
+ auto wellConnections = partitionIsEmpty || !wells ? Dune::cpgrid::WellConnections()
+ : Dune::cpgrid::WellConnections(*wells, possibleFutureConnections, grid);
+
+ if (cc.rank() == root) {
+ std::tie(rc, gIDtoRank) = applySerialZoltanCG(grid,
+ wellConnections,
+ cc.size(),
+ edgeWeightMethod,
+ root,
+ zoltanImbalanceTol,
+ allowDistributedWells,
+ params,
+ transGraph,
+ coarseThreshold,
+ coarsePartitionMaxNodeSize);
+ }
+
+ cc.broadcast(&rc, 1, root);
+ if (rc != ZOLTAN_OK) {
+ switch (rc) {
+ case ZOLTAN_OK+1:
+ OPM_THROW(std::runtime_error, "Could not initialize Zoltan!");
+ case ZOLTAN_OK+2:
+ OPM_THROW(std::runtime_error, "Could not create Zoltan!");
+ case ZOLTAN_OK+3:
+ OPM_THROW(std::runtime_error, "Partitioning with Zoltan failed!");
+ default:
+ OPM_THROW(std::runtime_error, "Unknown error reported by Zoltan!");
+ }
+ }
+
+ if (cc.rank() == root) {
+ // prepare exportedCells for communication
+ exportedCells = makeExportListsFromGIDtoRank(gIDtoRank, cc.size());
+ myImportList.reserve(exportedCells[root].size());
+ for (const auto& cell : exportedCells[root]) {
+ myImportList.emplace_back(cell, root, static_cast(AttributeSet::owner), -1);
+ }
+ // exclude root's own cells from communication
+ exportedCells[root].resize(0);
+ }
+ // communicate and create import+export lists
+ auto importedCells = Opm::Impl::communicateExportedCells(exportedCells, cc, root);
+ if (cc.rank() == root) {
+ myExportList.reserve(grid.numCells());
+ for (int i = 0; i < grid.numCells(); ++i) {
+ myExportList.emplace_back(i, gIDtoRank[i], static_cast(AttributeSet::owner));
+ }
+ } else {
+ myImportList.reserve(importedCells.size());
+ for (const auto& cell : importedCells) {
+ myImportList.emplace_back(cell, root, static_cast(AttributeSet::owner), -1);
+ }
+ }
+
+ // get the distribution of wells
+ std::vector> parallel_wells;
+ if (wells) {
+ if (allowDistributedWells) {
+ // wells can be split among several processes
+ auto wellsOnProc = Dune::cpgrid::perforatingWellIndicesOnProc(gIDtoRank, *wells, possibleFutureConnections, grid);
+ parallel_wells = Dune::cpgrid::computeParallelWells(wellsOnProc, *wells, cc, root);
+ } else {
+ // each well is guaranteed to be on a single process
+ auto wellRanks = getWellRanks(gIDtoRank, wellConnections);
+ parallel_wells = wellsOnThisRank(*wells, wellRanks, cc, root);
+ }
+ }
+
+ return std::make_tuple(std::move(gIDtoRank),
+ std::move(parallel_wells),
+ std::move(myExportList),
+ std::move(myImportList),
+ std::move(wellConnections));
+}
#endif // HAVE_MPI
// explicit template instantiations
diff --git a/opm/grid/GraphOfGridWrappers.hpp b/opm/grid/GraphOfGridWrappers.hpp
index 8e4d9c729..13b5116d5 100644
--- a/opm/grid/GraphOfGridWrappers.hpp
+++ b/opm/grid/GraphOfGridWrappers.hpp
@@ -26,9 +26,11 @@
#ifndef GRAPH_OF_GRID_WRAPPERS_HEADER
#define GRAPH_OF_GRID_WRAPPERS_HEADER
+#include
#include
#include
+#include
#include
#include // defines Zoltan and null-callback-functions
@@ -44,6 +46,7 @@ namespace Opm {
/// \brief callback function for ZOLTAN_NUM_OBJ_FN
///
/// returns the number of vertices in the graph
+template >
int getGraphOfGridNumVertices(void* pGraph, int *err);
/// \brief callback function for ZOLTAN_OBJ_LIST_FN
@@ -92,11 +95,58 @@ void getGraphOfGridEdgeList(void *pGraph,
float *edgeWeights,
int *err);
+/// \brief callback function for ZOLTAN_OBJ_LIST_FN
+///
+/// fills the vector gIDs with vertex global IDs
+/// and the vector objWeights with their weights
+void getCoarseGraphVerticesList(void* pGraph,
+ [[maybe_unused]] int dimGlobalID,
+ [[maybe_unused]] int dimLocalID,
+ ZOLTAN_ID_PTR gIDs,
+ [[maybe_unused]] ZOLTAN_ID_PTR lIDs,
+ int weightDim,
+ float *objWeights,
+ int *err);
+
+/// \brief callback function for ZOLTAN_NUM_EDGES_MULTI_FN
+///
+/// takes the list of global IDs (gIDs) and fills (consecutively)
+/// vector numEdges with the number of their edges
+void getCoarseGraphNumEdges(void *pGraph,
+ [[maybe_unused]] int dimGlobalID,
+ [[maybe_unused]] int dimLocalID,
+ int numCells,
+ ZOLTAN_ID_PTR gIDs,
+ [[maybe_unused]] ZOLTAN_ID_PTR lIDs,
+ int *numEdges,
+ int *err);
+
+/// \brief callback function for ZOLTAN_EDGE_LIST_MULTI_FN
+///
+/// takes the list of global IDs (gIDs) and fills (consecutively):
+/// vector nborGIDs with the list of neighbors (all into 1 vector),
+/// vector nborProc with neighbors' process numbers,
+/// vector edgeWeights with edge weights.
+/// The vector numEdges provides the number of edges for each gID
+void getCoarseGraphEdgeList(void *pGraph,
+ [[maybe_unused]] int dimGlobalID,
+ [[maybe_unused]] int dimLocalID,
+ int numCells,
+ ZOLTAN_ID_PTR gIDs,
+ [[maybe_unused]] ZOLTAN_ID_PTR lIDs,
+ int *numEdges,
+ ZOLTAN_ID_PTR nborGIDs,
+ int *nborProc,
+ int weightDim,
+ float *edgeWeights,
+ int *err);
+
/// \brief Register callback functions to Zoltan
-template
+template
void setGraphOfGridZoltanGraphFunctions(Zoltan_Struct *zz,
- GraphOfGrid& gog,
+ GraphType& gog,
bool pretendNull);
+
#endif
/// \brief Adds well to the GraphOfGrid
@@ -146,8 +196,9 @@ void extendAndSortImportList(std::vector>& importLi
/// On root, exportList is extended by well cells that are hidden from the partitioner.
/// These cells are also collected and returned so they can be communicated to other ranks.
/// \return vector[rank][cell] Each entry contains vector of cells exported to that rank.
+template
std::vector>
-extendRootExportList(const GraphOfGrid& gog,
+extendRootExportList(const GOG& gog,
std::vector>& exportList,
int root,
const std::vector& gIDtoRank);
@@ -254,6 +305,45 @@ makeImportAndExportLists(const GraphOfGrid& gog,
const Id* importGlobalGids,
int level);
+/// \brief Transform Zoltan output into tuples
+///
+/// \param gog CoarseGraphOfGrid, has ref. to CpGrid and has a coarsened partitioning graph
+/// \param cc Communication object
+/// \param wells Used to extract well names
+/// \param wellConnections Contains wells' global IDs, ordered as \param wells.
+/// \param possibleFutureConnections parameter needed if allowDistributedWells==true
+/// \param root Rank of the process executing the partitioning (usually 0)
+/// \param numExport Number of cells in the export list
+/// \param numImport Number of cells in the import list
+/// \param exportLocalGids Unused. Partitioning is performed on root
+/// process that has access to all cells.
+/// \param exportGlobalGids Zoltan output: Global IDs of exported cells
+/// \param exportToPart Zoltan output: ranks to which cells are exported
+/// \param importGlobalGids Zoltan output: Global IDs of cells imported to this rank
+/// \param allowDistributedWells should wells be allowed to exist on multiple partitions.
+/// \return gIDtoRank A vector indexed by global ID storing the rank of cell
+/// parallel_wells A vector of pairs wells.name and bool of "Is wells.name on this rank?"
+/// myExportList vector of cells to be moved from this rank
+/// myImportList vector of cells to be moved to this rank
+template
+std::tuple,
+ std::vector>,
+ std::vector >,
+ std::vector > >
+makeImportAndExportLists(const CoarseGraphOfGrid& gog,
+ const Dune::Communication& cc,
+ const std::vector * wells,
+ const Dune::cpgrid::WellConnections& wellConnections,
+ const std::unordered_map>& possibleFutureConnections,
+ int root,
+ int numExport,
+ int numImport,
+ [[maybe_unused]] const Id* exportLocalGids,
+ const Id* exportGlobalGids,
+ const int* exportToPart,
+ const Id* importGlobalGids,
+ bool allowDistributedWells);
+
/// \brief Call Zoltan partitioner on GraphOfGrid
///
/// GraphOfGrid represents a well by one vertex, so wells can not be
@@ -274,6 +364,27 @@ zoltanPartitioningWithGraphOfGrid(const Dune::CpGrid& grid,
const std::map& params,
int level);
+/// \brief Call Zoltan partitioner on Coarsened graph
+///
+/// CoarseGraphOfGrid incorporates transmissibility into the partitioning graph by
+/// coarsening the graph
+std::tuple, std::vector>,
+ std::vector >,
+ std::vector >,
+ Dune::cpgrid::WellConnections>
+zoltanPartitioningWithCoarseGraph(const Dune::CpGrid& grid,
+ const std::vector * wells,
+ const std::unordered_map>& possibleFutureConnections,
+ const Dune::cpgrid::CpGridDataTraits::Communication& cc,
+ Dune::EdgeWeightMethod edgeWeightMethod,
+ int root,
+ const double zoltanImbalanceTol,
+ bool allowDistributedWells,
+ const std::map& params,
+ const Dune::BCRSMatrix>* transGraph,
+ double coarseThreshold,
+ int coarsePartitionMaxNodeSize);
+
/// \brief Make complete export lists from a vector holding destination rank for each global ID
///
/// Intended to use on the root, as other ranks can not construct gIDtoRank.
@@ -302,6 +413,27 @@ zoltanSerialPartitioningWithGraphOfGrid(const Dune::CpGrid& grid,
const double zoltanImbalanceTol,
bool allowDistributedWells,
const std::map& params);
+
+/// \brief Call serial Zoltan partitioner on coarse Graph
+///
+/// CoarseGraphOfGrid incorporates transmissibility into the partitioning graph by
+/// coarsening the graph
+std::tuple, std::vector>,
+ std::vector >,
+ std::vector >,
+ Dune::cpgrid::WellConnections>
+zoltanSerialPartitioningWithCoarseGraph(const Dune::CpGrid& grid,
+ const std::vector * wells,
+ const std::unordered_map>& possibleFutureConnections,
+ const Dune::cpgrid::CpGridDataTraits::Communication& cc,
+ Dune::EdgeWeightMethod edgeWeightMethod,
+ int root,
+ const double zoltanImbalanceTol,
+ bool allowDistributedWells,
+ const std::map& params,
+ const Dune::BCRSMatrix>* transGraph,
+ double coarseThreshold,
+ int coarsePartitionMaxNodeSize);
#endif // HAVE_MPI
} // end namespace Opm
diff --git a/opm/grid/common/GridEnums.hpp b/opm/grid/common/GridEnums.hpp
index c7e097a96..a01bfce5e 100644
--- a/opm/grid/common/GridEnums.hpp
+++ b/opm/grid/common/GridEnums.hpp
@@ -49,7 +49,9 @@ namespace Dune {
/// \brief Use METIS for partitioning
metis=2,
/// \brief use Zoltan on GraphOfGrid for partitioning
- zoltanGoG=3
+ zoltanGoG=3,
+ /// \brief Use Zoltan coarse graph
+ zoltanCG=4
};
}
diff --git a/opm/grid/cpgrid/CpGrid.cpp b/opm/grid/cpgrid/CpGrid.cpp
index a77084222..ac68d3444 100644
--- a/opm/grid/cpgrid/CpGrid.cpp
+++ b/opm/grid/cpgrid/CpGrid.cpp
@@ -218,7 +218,10 @@ CpGrid::scatterGrid(EdgeWeightMethod method,
[[maybe_unused]] bool allowDistributedWells,
[[maybe_unused]] const std::vector& input_cell_part,
int level,
- [[maybe_unused]] bool useTransToFilterOverlap)
+ [[maybe_unused]] bool useTransToFilterOverlap,
+ [[maybe_unused]] const Dune::BCRSMatrix>* transGraph,
+ [[maybe_unused]] double coarseThreshold,
+ [[maybe_unused]] int coarsePartitionMaxNodeSize)
{
// Silence any unused argument warnings that could occur with various configurations.
static_cast(wells);
@@ -227,6 +230,9 @@ CpGrid::scatterGrid(EdgeWeightMethod method,
static_cast(method);
static_cast(imbalanceTol);
static_cast(level);
+ static_cast(transGraph);
+ static_cast(coarseThreshold);
+ static_cast(coarsePartitionMaxNodeSize);
if(!distributed_data_.empty())
{
@@ -381,8 +387,17 @@ CpGrid::scatterGrid(EdgeWeightMethod method,
#else
OPM_THROW(std::runtime_error, "Parallel runs depend on ZOLTAN if useZoltan is true. Please install!");
#endif // HAVE_ZOLTAN
- }
- else
+ } else if (partitionMethod == Dune::PartitionMethod::zoltanCG)
+ {
+#ifdef HAVE_ZOLTAN
+ std::tie(computedCellPart, wells_on_proc, exportList, importList, wellConnections)
+ = serialPartitioning
+ ? Opm::zoltanSerialPartitioningWithCoarseGraph(*this, wells, possibleFutureConnections, cc, method, 0, imbalanceTol, allowDistributedWells, partitioningParams, transGraph, coarseThreshold, coarsePartitionMaxNodeSize)
+ : Opm::zoltanPartitioningWithCoarseGraph(*this, wells, possibleFutureConnections, cc, method, 0, imbalanceTol, allowDistributedWells, partitioningParams, transGraph, coarseThreshold, coarsePartitionMaxNodeSize);
+#else
+ OPM_THROW(std::runtime_error, "Parallel runs depend on ZOLTAN if useZoltan is true. Please install!");
+#endif // HAVE_ZOLTAN
+ } else
{
std::tie(computedCellPart, wells_on_proc, exportList, importList, wellConnections) =
cpgrid::vanillaPartitionGridOnRoot(*this, wells, possibleFutureConnections, transmissibilities, allowDistributedWells);
diff --git a/tests/test_coarsegraphofgrid.cpp b/tests/test_coarsegraphofgrid.cpp
new file mode 100644
index 000000000..1b2cf7914
--- /dev/null
+++ b/tests/test_coarsegraphofgrid.cpp
@@ -0,0 +1,621 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*
+ Copyright 2026 Equinor ASA.
+
+ This file is part of the Open Porous Media project (OPM).
+
+ OPM is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OPM is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OPM. If not, see .
+
+ Consult the COPYING file in the top-level source directory of this
+ module for the precise wording of the license and the list of
+ copyright holders.
+*/
+
+#include
+
+#include
+
+#define BOOST_TEST_MODULE GraphRepresentationOfGrid
+#define BOOST_TEST_NO_MAIN
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+
+#if HAVE_OPM_COMMON
+#include
+#include
+#include
+#include
+#include
+#include
+#endif
+
+#include
+
+// Create "Transmissibility graph" from grid and set all transmissibilities to 1.0
+std::unique_ptr>> createAdjacencyMatrix(const Dune::CpGrid& grid)
+{
+ size_t numCells = grid.numCells();
+ auto gridView = grid.leafGridView();
+
+ using ElementMapper =
+ Dune::MultipleCodimMultipleGeomTypeMapper;
+ const auto elemMapper = ElementMapper { gridView, Dune::mcmgElementLayout() };
+
+ Dune::MatrixIndexSet op;
+ op.resize(numCells, numCells);
+
+ // Create adjacency for transmissibility graph
+ for (const auto& elem : elements(gridView, Dune::Partitions::interiorBorder)) {
+ auto d = elemMapper.index(elem);
+ op.add(d, d);
+ for (const auto& is : intersections(gridView, elem)) {
+ if (!is.neighbor()) {
+ continue;
+ }
+
+ const auto J = static_cast(elemMapper.index(is.outside()));
+ op.add(d, J);
+ }
+ }
+
+ // Allocate the BCRSMatrix wrapped in a unique_ptr
+ auto graph = std::make_unique>>();
+ op.exportIdx(*graph);
+ *graph = 1.0;
+
+ return graph;
+}
+
+// basic test to check if the graph was constructed correctly. Copied from test_graphofgrid.cpp
+BOOST_AUTO_TEST_CASE(SimpleGraph)
+{
+ Dune::CpGrid grid;
+ std::array dims{2,2,2};
+ std::array size{2.,2.,2.};
+ grid.createCartesian(dims,size);
+ if (grid.size(0)==0)
+ return;
+
+ Dune::cpgrid::WellConnections wells;
+
+ auto graph = createAdjacencyMatrix(grid);
+ double strongConnectionThreshold = 1.1; // no nodes will be merged
+ int maxNodeSize = 2;
+ bool allowDistWells = false;
+ Opm::CoarseGraphOfGrid cgog(grid, Dune::EdgeWeightMethod::uniformEdgeWgt,
+ graph.get(), strongConnectionThreshold,
+ maxNodeSize, allowDistWells, 0, wells);
+
+ BOOST_REQUIRE(cgog.size()==8); // number of graph vertices
+ BOOST_REQUIRE(cgog.getCoarseNodes()[0].size()==1); // weight of all nodes is one
+
+ const std::vector< std::map > edges = cgog.getCoarseEdges();
+
+ BOOST_REQUIRE(edges[0].size()==3);
+ auto edgeL = edges[2];
+
+ BOOST_REQUIRE(edgeL.size()==3); // neighbors of vertex 2 are: 0, 3, 6
+ BOOST_REQUIRE(edgeL[0]==1.0);
+ BOOST_REQUIRE(edgeL[3]==1.0);
+ BOOST_REQUIRE(edgeL[6]==1.0);
+ BOOST_REQUIRE_THROW(edgeL.at(4),std::out_of_range); // not a neighbor
+
+ BOOST_REQUIRE_THROW(edges.at(10),std::logic_error);
+}
+
+// basic test to check if nodes get merged.
+//Similar to SimpleGraphWithVertexContraction in test_graphofgrid.cpp
+BOOST_AUTO_TEST_CASE(MergeTwoNodes)
+{
+ Dune::CpGrid grid;
+ std::array dims{2,2,2};
+ std::array size{2.,2.,2.};
+ grid.createCartesian(dims,size);
+ if (grid.size(0)==0)
+ return;
+
+ Dune::cpgrid::WellConnections wells;
+
+ auto graph = createAdjacencyMatrix(grid);
+ (*graph)[0][1] = 2.0; // Set 0->1 connection to 2
+ double strongConnectionThreshold = 1.1; // Node 0 and 1 will be merged
+ int maxNodeSize = 2;
+ bool allowDistWells = false;
+ Opm::CoarseGraphOfGrid cgog(grid, Dune::EdgeWeightMethod::uniformEdgeWgt,
+ graph.get(), strongConnectionThreshold,
+ maxNodeSize, allowDistWells, 0, wells);
+
+ BOOST_REQUIRE(cgog.size()==7); // number of graph vertices
+ BOOST_REQUIRE(cgog.getCoarseNodes()[0].size()==2); // weight of node 0 is two
+
+ const auto edges = cgog.getCoarseEdges();
+
+ auto edge0 = edges[0];
+ BOOST_REQUIRE(edge0.size()==4);
+
+ auto map_to_coarse = cgog.getMapToCoarse();
+
+ BOOST_REQUIRE(map_to_coarse[0]==0);
+ BOOST_REQUIRE(map_to_coarse[1]==0);
+ BOOST_REQUIRE(map_to_coarse[2]==1);
+
+ BOOST_REQUIRE(edge0[map_to_coarse[2]]==1.0);
+ BOOST_REQUIRE(edge0[map_to_coarse[3]]==1.0);
+ BOOST_REQUIRE(edge0[map_to_coarse[4]]==1.0);
+ BOOST_REQUIRE(edge0[map_to_coarse[5]]==1.0);
+}
+
+// basic test to check if the coarse graph gets correct edge weights
+BOOST_AUTO_TEST_CASE(CoarseEdgeWeights)
+{
+ Dune::CpGrid grid;
+ std::array dims{3,3,1};
+ std::array size{2.,2.,2.};
+ grid.createCartesian(dims,size);
+ if (grid.size(0)==0)
+ return;
+
+ Dune::cpgrid::WellConnections wells;
+
+ auto graph = createAdjacencyMatrix(grid);
+ (*graph)[0][1] = 2.0;
+ (*graph)[0][3] = 2.0;
+ (*graph)[7][8] = 0.0; // Fault between cell 7 and 8
+ (*graph)[8][7] = 0.0;
+
+ double strongConnectionThreshold = 1.1;
+ int maxNodeSize = 3;
+ bool allowDistWells = false;
+ Opm::CoarseGraphOfGrid cgog(grid, Dune::EdgeWeightMethod::uniformEdgeWgt,
+ graph.get(), strongConnectionThreshold,
+ maxNodeSize, allowDistWells, 0, wells);
+
+ BOOST_REQUIRE(cgog.size()==7); // number of graph vertices
+ BOOST_REQUIRE(cgog.getCoarseNodes()[0].size()==3); // Node weight of node 0 is 3
+
+ const std::vector< std::map > edges = cgog.getCoarseEdges();
+ auto edge0 = edges[0];
+
+ BOOST_REQUIRE(edge0.size()==3); // The merged node has three edges
+
+ auto map_to_coarse = cgog.getMapToCoarse();
+ BOOST_REQUIRE(edge0[map_to_coarse[2]]==1.0);
+ BOOST_REQUIRE(edge0[map_to_coarse[4]]==2.0); //Node 0 has two connections to the centre node
+ BOOST_REQUIRE(edge0[map_to_coarse[6]]==1.0);
+
+ auto edge7 = edges[map_to_coarse[7]];
+ auto edge8 = edges[map_to_coarse[8]];
+
+ BOOST_REQUIRE(edge7.size()==2); // Only two edges because of fault
+ BOOST_REQUIRE(edge8.size()==1); // Only one edges because of fault
+
+ BOOST_REQUIRE_THROW(edge7.at(map_to_coarse[8]),std::out_of_range); // not a neighbor
+ BOOST_REQUIRE_THROW(edge8.at(map_to_coarse[7]),std::out_of_range); // not a neighbor
+}
+
+
+// basic test to check if the graph uses maxNodeSize correctly
+BOOST_AUTO_TEST_CASE(MaxNodeSize)
+{
+ Dune::CpGrid grid;
+ std::array dims{3,3,1};
+ std::array size{2.,2.,2.};
+ grid.createCartesian(dims,size);
+ if (grid.size(0)==0)
+ return;
+
+ Dune::cpgrid::WellConnections wells;
+
+ auto graph = createAdjacencyMatrix(grid);
+ (*graph)[0][1] = 2.0;
+ (*graph)[0][3] = 3.0; // This is lager, so should be merged
+ (*graph)[6][7] = 2.0;
+ (*graph)[7][8] = 2.0;
+
+ double strongConnectionThreshold = 1.1;
+ int maxNodeSize = 2;
+ bool allowDistWells = false;
+ Opm::CoarseGraphOfGrid cgog(grid, Dune::EdgeWeightMethod::uniformEdgeWgt,
+ graph.get(), strongConnectionThreshold,
+ maxNodeSize, allowDistWells, 0, wells);
+
+ BOOST_REQUIRE(cgog.size()==7); // number of graph vertices
+ BOOST_REQUIRE(cgog.getCoarseNodes()[0].size()==2); // weight of node 0 is two
+
+ auto map_to_coarse = cgog.getMapToCoarse();
+ BOOST_REQUIRE(map_to_coarse[0]==map_to_coarse[3]); // node 0 merged with node 3
+ BOOST_REQUIRE(map_to_coarse[0]!=map_to_coarse[1]); // node 1 not merged with node 0 and 3
+
+ BOOST_REQUIRE(cgog.getCoarseNodes()[map_to_coarse[6]].size()==2); // weight of node 6 is two
+ BOOST_REQUIRE(map_to_coarse[6]==map_to_coarse[7]); // node 6 merged with node 7
+ BOOST_REQUIRE(map_to_coarse[7]!=map_to_coarse[8]); // node 8 not merged with node 6 and 7
+}
+
+#if HAVE_OPM_COMMON
+namespace {
+ // create Wells, we only use well name and cell locations
+ // copied from test_graphofgrid.cpp
+ auto createConnection (int i, int j, int k)
+ {
+ return Opm::Connection(i,j,k,0, 0,Opm::Connection::State::OPEN,
+ Opm::Connection::Direction::Z,
+ Opm::Connection::CTFKind::DeckValue, 0,
+ 5.,Opm::Connection::CTFProperties(),0,false);
+ }
+ auto createWell (const std::string& name)
+ {
+ using namespace Opm;
+ return Dune::cpgrid::OpmWellType(name,name,0,0,0,0,0.,WellType(),
+ Well::ProducerCMode(),Connection::Order(),UnitSystem(),
+ 0.,0.,false,false,0,Well::GasInflowEquation());
+ };
+} // end anonymous namespace
+#endif
+
+#if HAVE_MPI && HAVE_OPM_COMMON
+// Check if merging wells work
+// Similar to addWellConnections in test_graphofgrid.cpp
+BOOST_AUTO_TEST_CASE(MergeWells)
+{
+ Dune::CpGrid grid;
+ std::array dims{2,2,2};
+ std::array size{1.,1.,1.};
+ grid.createCartesian(dims,size);
+ if (grid.size(0)==0)
+ return;
+
+ auto wellCon = std::make_shared(); // do not confuse with Dune::cpgrid::WellConnections
+ wellCon->add(createConnection(0,0,0)); // 0
+ wellCon->add(createConnection(0,1,0)); // 2
+ wellCon->add(createConnection(0,1,1)); // 6
+ std::vector wells;
+ wells.push_back(createWell("first"));
+ wells[0].updateConnections(wellCon,true);
+
+ wellCon = std::make_shared(); //reset
+ wellCon->add(createConnection(0,0,1)); // 4
+ wellCon->add(createConnection(1,1,0)); // 3
+ wells.push_back(createWell("second"));
+ wells[1].updateConnections(wellCon,true);
+
+ wellCon = std::make_shared(); //reset
+ wellCon->add(createConnection(0,0,1)); // 4
+ wellCon->add(createConnection(1,0,1)); // 5
+ wells.push_back(createWell("third")); // intersects with second
+ wells[2].updateConnections(wellCon,true);
+
+ Dune::cpgrid::WellConnections wellConnections(wells,std::unordered_map>(),grid);
+
+ auto graph = createAdjacencyMatrix(grid);
+ (*graph)[0][1] = 2.0; // should not be merged, due to maxNodeSize = 2
+
+ double strongConnectionThreshold = 1.1;
+ int maxNodeSize = 2;
+ bool allowDistWells = false;
+ Opm::CoarseGraphOfGrid cgog(grid, Dune::EdgeWeightMethod::uniformEdgeWgt,graph.get(),
+ strongConnectionThreshold,
+ maxNodeSize, allowDistWells, 0, wellConnections);
+
+ BOOST_REQUIRE(cgog.size()==4);
+
+ int err;
+ int nVer = Opm::getGraphOfGridNumVertices >(&cgog,&err);
+ BOOST_REQUIRE(err==ZOLTAN_OK);
+ BOOST_REQUIRE(nVer == 4);
+ std::vector gIDs(nVer);
+ std::vector objWeights(nVer);
+ getCoarseGraphVerticesList(&cgog, 1, 1, gIDs.data(), nullptr, 1, objWeights.data(), &err);
+ BOOST_REQUIRE(err==ZOLTAN_OK);
+ std::ranges::sort(gIDs);
+ BOOST_REQUIRE(objWeights[0]==3.0 && objWeights[1]==1.0 && objWeights[2]==3.0 && objWeights[3]==1.0);
+ std::vector numEdges(nVer);
+ getCoarseGraphNumEdges(&cgog, 1, 1, nVer, gIDs.data(), nullptr, numEdges.data(), &err);
+ BOOST_REQUIRE(err==ZOLTAN_OK);
+ BOOST_REQUIRE(numEdges[0]==3 && numEdges[1]==2 && numEdges[2]==3 && numEdges[3]==2);
+ int nEdges = 10; // sum of numEdges[i]
+ std::vector nborGIDs(nEdges);
+ std::vector nborProc(nEdges);
+ std::vector edgeWeights(nEdges);
+ getCoarseGraphEdgeList(&cgog, 1, 1, nVer, gIDs.data(), nullptr, numEdges.data(), nborGIDs.data(), nborProc.data(), 1, edgeWeights.data(), &err);
+ BOOST_REQUIRE(err==ZOLTAN_OK);
+
+ // check all edgeWeights. Note that nborGIDs are not sorted
+ for (int i=0; i<3; ++i)
+ {
+ switch (nborGIDs[i])
+ {
+ case 1: BOOST_REQUIRE(edgeWeights[i]==1); break;
+ case 2: BOOST_REQUIRE(edgeWeights[i]==3); break;
+ case 3: BOOST_REQUIRE(edgeWeights[i]==1); break;
+ default: throw("CoarseGraph was constructed badly.");
+ }
+ }
+ for (int i=3; i<5; ++i)
+ {
+ switch (nborGIDs[i])
+ {
+ case 0: BOOST_REQUIRE(edgeWeights[i]==1); break;
+ case 2: BOOST_REQUIRE(edgeWeights[i]==2); break;
+ default: throw("CoarseGraph was constructed badly.");
+ }
+ }
+ for (int i=5; i<8; ++i)
+ {
+ switch (nborGIDs[i])
+ {
+ case 0: BOOST_REQUIRE(edgeWeights[i]==3); break;
+ case 1: BOOST_REQUIRE(edgeWeights[i]==2); break;
+ case 3: BOOST_REQUIRE(edgeWeights[i]==2); break;
+ default: throw("CoarseGraph was constructed badly.");
+ }
+ }
+ for (int i=8; i<10; ++i)
+ {
+ switch (nborGIDs[i])
+ {
+ case 0: BOOST_REQUIRE(edgeWeights[i]==1); break;
+ case 2: BOOST_REQUIRE(edgeWeights[i]==2); break;
+ default: throw("CoarseGraph was constructed badly.");
+ }
+ }
+}
+
+// Check if merging multiple intersecting wells work
+// Similar to IntersectingWells in test_graphofgrid.cpp
+BOOST_AUTO_TEST_CASE(MergeWellsMoreIntersecting)
+{
+ Dune::CpGrid grid;
+ std::array dims{5,4,3};
+ std::array size{1.,1.,1.};
+ grid.createCartesian(dims,size);
+ if (grid.size(0)==0)
+ return;
+
+ auto wellCon = std::make_shared(); // do not confuse with Dune::cpgrid::WellConnections
+ wellCon->add(createConnection(0,0,0)); // 0
+ wellCon->add(createConnection(1,0,0)); // 1
+ wellCon->add(createConnection(2,0,0)); // 2
+ wellCon->add(createConnection(3,0,0)); // 3
+ wellCon->add(createConnection(4,0,0)); // 4
+ std::vector wells;
+ wells.push_back(createWell("first"));
+ wells[0].updateConnections(wellCon,true);
+
+ wellCon = std::make_shared(); //reset
+ wellCon->add(createConnection(2,2,2)); // 52
+ wellCon->add(createConnection(2,2,1)); // 32
+ wellCon->add(createConnection(2,2,0)); // 12
+ wells.push_back(createWell("second"));
+ wells[1].updateConnections(wellCon,true);
+
+ wellCon = std::make_shared(); //reset
+ wellCon->add(createConnection(4,3,2)); // 59
+ wellCon->add(createConnection(3,1,2)); // 48
+ wellCon->add(createConnection(2,3,1)); // 37
+ wells.push_back(createWell("third")); //
+ wells[2].updateConnections(wellCon,true);
+
+ wellCon = std::make_shared(); //reset
+ wellCon->add(createConnection(2,3,1)); // 37
+ wellCon->add(createConnection(3,3,1)); // 38
+ wellCon->add(createConnection(4,3,1)); // 39
+ wellCon->add(createConnection(4,2,1)); // 34
+ wells.push_back(createWell("forth"));
+ wells[3].updateConnections(wellCon,true);
+
+ wellCon = std::make_shared