diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 0d31c27e3..df18b82b0 100644 --- a/CMakeLists_files.cmake +++ b/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 diff --git a/opm/grid/GraphOfGrid.cpp b/opm/grid/GraphOfGrid.cpp index 8b49cd821..3c3d5ce98 100644 --- a/opm/grid/GraphOfGrid.cpp +++ b/opm/grid/GraphOfGrid.cpp @@ -356,6 +356,37 @@ void GraphOfGrid::addNeighboringCellsToWells () } } +template +void GraphOfGrid::multiplyWellConnectivity (const std::set& 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; } // namespace Opm diff --git a/opm/grid/GraphOfGrid.hpp b/opm/grid/GraphOfGrid.hpp index 4784c1e02..ec3618ae5 100644 --- a/opm/grid/GraphOfGrid.hpp +++ b/opm/grid/GraphOfGrid.hpp @@ -179,6 +179,9 @@ class GraphOfGrid{ } } + /// \brief Multiply edges between well cells by a factor + void multiplyWellConnectivity (const std::set& well, const WeightType& factor); + private: /// \brief Create a graph representation of the grid /// diff --git a/opm/grid/GraphOfGridWrappers.cpp b/opm/grid/GraphOfGridWrappers.cpp index 34e540df6..197b9ae47 100644 --- a/opm/grid/GraphOfGridWrappers.cpp +++ b/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 diff --git a/opm/grid/common/MetisPartition.cpp b/opm/grid/common/MetisPartition.cpp index a89021d86..a24709548 100644 --- a/opm/grid/common/MetisPartition.cpp +++ b/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(cpgrid, *gridAndWells); + idx_t weWeight = calculateWellEdgeWeight(cpgrid, *gridAndWells); int neighborCounter = 0; for( int cell = 0; cell < n; cell++ ) diff --git a/opm/grid/common/ZoltanGraphFunctions.cpp b/opm/grid/common/ZoltanGraphFunctions.cpp index 63005ab05..62e73d90b 100644 --- a/opm/grid/common/ZoltanGraphFunctions.cpp +++ b/opm/grid/common/ZoltanGraphFunctions.cpp @@ -180,13 +180,29 @@ void getNullEdgeList(void *cpGridPointer, int sizeGID, int sizeLID, } template -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::infinity()) { + total /= grid.numFaces(); + } else { + // grid is too big, use maximum instead of the average + total = 0; + for (int edge=0; edge(std::numeric_limits::max()); if (total > maxVal) { return std::numeric_limits::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(grid, graph); + float weWeight = calculateWellEdgeWeight(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 diff --git a/opm/grid/common/ZoltanGraphFunctions.hpp b/opm/grid/common/ZoltanGraphFunctions.hpp index c1dd3a075..3add0e9a4 100644 --- a/opm/grid/common/ZoltanGraphFunctions.hpp +++ b/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 -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 diff --git a/opm/grid/common/ZoltanPartition.cpp b/opm/grid/common/ZoltanPartition.cpp index 815a11edf..f64dc0ce6 100644 --- a/opm/grid/common/ZoltanPartition.cpp +++ b/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. diff --git a/tests/cpgrid/combinedgridwellgraph_test.cpp b/tests/cpgrid/combinedgridwellgraph_test.cpp new file mode 100644 index 000000000..54f6f01a4 --- /dev/null +++ b/tests/cpgrid/combinedgridwellgraph_test.cpp @@ -0,0 +1,174 @@ +// -*- 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 CombinedGridWellGraph +#define BOOST_TEST_NO_MAIN +#include +#include +#include +#include +#include + +#if HAVE_OPM_COMMON +#include +#include +#include +#include +#include +#include +#endif + +#include +#include + +#if HAVE_OPM_COMMON +namespace { + // create Wells, we only use well name and cell locations + 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 + +BOOST_AUTO_TEST_CASE(CombinedGridWellGraph) +{ + ///! beware, test functionality, not implementation + + /// construct CombinedGridWellGraph + Dune::CpGrid grid; + std::array dims { 3, 3, 1 }; + std::array size { 1., 1., 1. }; + grid.createCartesian(dims, size); + + auto wellCon = std::make_shared(); // do not confuse with Dune::cpgrid::WellConnections + std::vector wells; + wellCon->add(createConnection(0,0,0)); + wellCon->add(createConnection(1,0,0)); + wellCon->add(createConnection(2,0,0)); + wells.push_back(createWell("first_row")); + wells[0].updateConnections(wellCon,true); + + wellCon = std::make_shared(); // reset + wellCon->add(createConnection(0,0,0)); + wellCon->add(createConnection(1,1,0)); + wellCon->add(createConnection(2,2,0)); + wells.push_back(createWell("diag")); + wells[1].updateConnections(wellCon,true); + + // std::vector wells; + std::unordered_map> possibleFutureConnections; + std::vector transmissibilities{0., 1., 2., 0., 0., 3., 4., 0., 0., 5., 6., 0., // X + 0., 0., 0.,11.,12.,13.,14.,15.,16., 0., 0., 0., // Y + 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}; // Z + Dune::cpgrid::CombinedGridWellGraph gridWellGraph(grid, &wells, possibleFutureConnections, + transmissibilities.data(), false, Dune::EdgeWeightMethod::defaultTransEdgeWgt); + const auto& wellEdges = gridWellGraph.getWellsGraph(); + + /// log transmissibilities do not affect well connections + { + Dune::cpgrid::CombinedGridWellGraph logGridWellGraph(grid, &wells, possibleFutureConnections, + transmissibilities.data(), false, Dune::EdgeWeightMethod::logTransEdgeWgt); + const auto& logWellEdges = logGridWellGraph.getWellsGraph(); + BOOST_REQUIRE(wellEdges == logWellEdges); + } + + /// connections of every well are interconnected + BOOST_REQUIRE(wellEdges.size() == 9); // all vertices + std::set edges{1, 2, 4, 8}; // both wells have the cell 0 + BOOST_REQUIRE(wellEdges[0].size() == 4); + BOOST_REQUIRE(wellEdges[0] == edges); + edges = std::set{0,2}; + BOOST_REQUIRE(wellEdges[1].size() == 2); + BOOST_REQUIRE(wellEdges[1] == edges); + edges = std::set{0,1}; + BOOST_REQUIRE(wellEdges[2].size() == 2); + BOOST_REQUIRE(wellEdges[2] == edges); + BOOST_REQUIRE(wellEdges[3].size() == 0); + edges = std::set{0,8}; + BOOST_REQUIRE(wellEdges[4].size() == 2); + BOOST_REQUIRE(wellEdges[4] == edges); + BOOST_REQUIRE(wellEdges[5].size() == 0); + BOOST_REQUIRE(wellEdges[6].size() == 0); + BOOST_REQUIRE(wellEdges[7].size() == 0); + edges = std::set{0,4}; + BOOST_REQUIRE(wellEdges[8].size() == 2); + BOOST_REQUIRE(wellEdges[8] == edges); + + /// well edge weight is by defaul equal to the sum of all grid edges + auto wellEdgeWeight = calculateWellEdgeWeight(grid, gridWellGraph); + // the sum (102) gets multipied by 1e18 to accomodate partitioners that use integral weights (Metis) + BOOST_REQUIRE(std::abs(wellEdgeWeight / 1e18 - 102.) < 1e-5); + + // a positive multiplier changes the weight from sum_of_faces to multiplier_\times_average_face + gridWellGraph.setMultiplyWellConnectivities(21.); // there are 42 grid edges in total + wellEdgeWeight = calculateWellEdgeWeight(grid, gridWellGraph); + BOOST_REQUIRE(std::abs(wellEdgeWeight / 1e18 - 51.) < 1e-5); + + /// edges from well connections overwrite grid edges + int neighborCounter{0}; + std::vector gID(9), nborGID(32, 0); // 32=2x(12 inner faces + 4 new connections due to wells) + ZOLTAN_ID_PTR nborGIDData = nborGID.data(), gIDData = gID.data(); + std::vector ewgts(32, 0.); + for (int i=0; i<9; ++i) { + gID[i] = i; + fillNBORGIDAndWeightsForSpecificCellAndIncrementNeighborCounterForGridWithWells(gridWellGraph, i, gIDData, neighborCounter, nborGIDData, ewgts.data(), wellEdgeWeight); + } + + BOOST_REQUIRE(neighborCounter == 32); + // check if the well edge weights are at correct places + for (const auto i : std::vector{0,1,2,3, 5,6, 8,9, 14,15, 28,29}) { + BOOST_REQUIRE(ewgts[i] == wellEdgeWeight); + } + const auto weightSum = std::accumulate(ewgts.begin(), ewgts.end(), 0.); + // sum of edge weights is 102, there are 6 well connections in total and grid edges {1,2} were overwritten + BOOST_REQUIRE(std::abs(weightSum/1e18 - 2*(102 + 6*wellEdgeWeight/1e18 - 1 - 2)) < 1e-3); +} + +bool init_unit_test_func() +{ + return true; +} + +int main(int argc, char** argv) +{ + Dune::MPIHelper::instance(argc, argv); + boost::unit_test::unit_test_main(&init_unit_test_func, + argc, argv); +} diff --git a/tests/test_graphofgrid.cpp b/tests/test_graphofgrid.cpp index a18f21462..17402ac20 100644 --- a/tests/test_graphofgrid.cpp +++ b/tests/test_graphofgrid.cpp @@ -892,6 +892,61 @@ BOOST_AUTO_TEST_CASE(test_getWellRanks) } #endif +BOOST_AUTO_TEST_CASE(MultiplyWellConnectivities) +{ + Dune::CpGrid grid; + std::array dims { 3, 3, 1 }; + std::array size { 1., 1., 1. }; + grid.createCartesian(dims, size); + Opm::GraphOfGrid gog(grid); + if (grid.size(0) == 0) + return; + + std::unordered_map> wells { + { "first row", { 0, 1, 2 } }, + { "first column", { 0, 3, 6 } }, + { "intersecting first row", { 1, 2, 4, 5 } } + }; + float factor = 5; + for (const auto& well : wells) { + // multiplies weights of well edges, not between two wells + gog.multiplyWellConnectivity(well.second, factor); + } + BOOST_REQUIRE(gog.size() == 9); + BOOST_REQUIRE(gog.getWells().size() == 0); + int err; + int nVer = getGraphOfGridNumVertices(&gog, &err); + BOOST_REQUIRE(err == ZOLTAN_OK); + BOOST_REQUIRE(nVer == 9); + + auto checkEdge = [&gog](int from, int to, float weight) { + const auto& edgesFrom = gog.edgeList(from); + BOOST_REQUIRE(edgesFrom.at(to) == weight); + const auto& edgesTo = gog.edgeList(to); + BOOST_REQUIRE(edgesTo.at(from) == weight); + }; + // left-to-right edges + checkEdge(0, 1, 5.); + checkEdge(1, 2, 25.); + checkEdge(3, 4, 1.); + checkEdge(4, 5, 5.); + checkEdge(6, 7, 1.); + checkEdge(7, 8, 1.); + // front-to-back edges + checkEdge(0, 3, 5.); + checkEdge(1, 4, 5.); + checkEdge(2, 5, 5.); + checkEdge(3, 6, 5.); + checkEdge(4, 7, 1.); + checkEdge(5, 8, 1.); + + int nrEdges = 0; + for (int i=0; i