Skip to content

Commit

Permalink
Merge pull request #809 from michal-toth/feature/graph-of-grid-second…
Browse files Browse the repository at this point in the history
…-try

Use transmissibilities as edge weights
  • Loading branch information
blattms authored Dec 3, 2024
2 parents 31946cc + 29072ca commit 73423a7
Show file tree
Hide file tree
Showing 5 changed files with 201 additions and 48 deletions.
45 changes: 35 additions & 10 deletions opm/grid/GraphOfGrid.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,7 @@
namespace Opm{

template<typename Grid>
void GraphOfGrid<Grid>::createGraph ()
void GraphOfGrid<Grid>::createGraph (const double* transmissibilities)
{
const auto& rank = grid.comm().rank();
// load vertices (grid cells) into graph
Expand Down Expand Up @@ -57,11 +57,11 @@ void GraphOfGrid<Grid>::createGraph ()
{
continue;
}
WeightType weight = 1; // default edge weight
vertex.edges.try_emplace(otherCell,weight);
WeightType weight = transmissibilities ? transmissibilities[face] : 1; // default edge weight is 1
vertex.edges.try_emplace(otherCell, weight);
}

graph.try_emplace(gID,vertex);
graph.try_emplace(gID, vertex);
}

}
Expand All @@ -87,7 +87,7 @@ int GraphOfGrid<Grid>::contractVertices (int gID1, int gID2)
}
if (gID2<gID1)
{
std::swap(gID1,gID2);
std::swap(gID1, gID2);
}

// add up vertex weights
Expand All @@ -107,7 +107,7 @@ int GraphOfGrid<Grid>::contractVertices (int gID1, int gID2)
v1e.insert(edge);
// remap neighbor's edge
graph[edge.first].edges.erase(gID2);
graph[edge.first].edges.emplace(gID1,edge.second);
graph[edge.first].edges.emplace(gID1, edge.second);
}
}
else
Expand Down Expand Up @@ -164,26 +164,51 @@ void GraphOfGrid<Grid>::addWell (const std::set<int>& well, bool checkIntersecti
wID = *(w->begin());
}
gID = *(w->begin());
newWell.insert(w->begin(),w->end());
newWell.insert(w->begin(), w->end());
wells.erase(w);
break; // GraphOfGrid::wells are constructed to be disjoint, each gID has max 1 match
}
}
wID = contractVertices(wID,gID);
wID = contractVertices(wID, gID);
assert(wID!=-1 && "Added well vertex was not found in the grid (or its wells).");
}
newWell.insert(well.begin(),well.end());
newWell.insert(well.begin(), well.end());
wells.push_front(newWell);
}
else
{
for (int gID : well)
{
wID = contractVertices(wID,gID);
wID = contractVertices(wID, gID);
}
wells.emplace_front(well);
}
}

template<typename Grid>
void GraphOfGrid<Grid>::addNeighboringCellsToWells ()
{
// mark all cells that will be added to wells (addding them one
// by one would require recursive checks for neighboring wells)
std::vector<std::set<int>> buffer(wells.size());
int i=0;
for (auto& w : wells)
{
buffer[i].insert(*w.begin()); // intersects with its well
for (const auto& v : edgeList(*w.begin()))
{
buffer[i].insert(v.first);
}
++i;
}

// intersecting wells and buffers will be merged
for (const auto& b : buffer)
{
addWell(b);
}
}


template class GraphOfGrid<Dune::CpGrid>;
} // namespace Opm
21 changes: 18 additions & 3 deletions opm/grid/GraphOfGrid.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -52,10 +52,10 @@ class GraphOfGrid{
};

public:
explicit GraphOfGrid (const Grid& grid_)
explicit GraphOfGrid (const Grid& grid_, const double* transmissibilities=nullptr)
: grid(grid_)
{
createGraph();
createGraph(transmissibilities);
}

const Grid& getGrid() const
Expand Down Expand Up @@ -162,9 +162,24 @@ class GraphOfGrid{
return wells;
}

/// \brief Contract a layer of verices around each well into it
///
/// Representing a well by one node guarantees that the well won't
/// be split over several processes. Giving the well an extra layer
/// of cells distances that well from the subdomain boundary.
void addNeighboringCellsToWells ();
void addNeighboringCellsToWells (int layers)
{
for (int i=0; i<layers; ++i)
{
addNeighboringCellsToWells();
}
}

private:
/// \brief Create a graph representation of the grid
void createGraph (); // edge weight=1
/// If transmissibilities are not supplied, edge weight=1
void createGraph (const double* transmissibilities=nullptr);

/// \brief Identify the well containing the cell with this global ID
///
Expand Down
64 changes: 32 additions & 32 deletions opm/grid/GraphOfGridWrappers.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -174,7 +174,7 @@ void addFutureConnectionWells(GraphOfGrid<Dune::CpGrid>& gog,
assert(gID!=-1); // well should be an active cell
wellsgID.insert(gID);
}
gog.addWell(wellsgID,checkWellIntersections);
gog.addWell(wellsgID, checkWellIntersections);
}
}

Expand All @@ -184,7 +184,7 @@ void addWellConnections(GraphOfGrid<Dune::CpGrid>& gog,
{
for (const auto& w : wells)
{
gog.addWell(w,checkWellIntersections);
gog.addWell(w, checkWellIntersections);
}
}

Expand Down Expand Up @@ -225,20 +225,20 @@ extendRootExportList(const GraphOfGrid<Dune::CpGrid>& gog,
using ExportList = std::vector<std::tuple<int,int,char>>;
// make a list of wells for easy identification. Contains ID, begin, end
using iter = std::set<int>::const_iterator;
std::unordered_map<int,std::tuple<iter,iter,int>> wellMap;
std::unordered_map<int, std::tuple<iter,iter,int>> wellMap;
for (const auto& well : gog.getWells())
{
if (gIDtoRank.size()>0)
{
auto wellID = *well.begin();
if (gIDtoRank[wellID]!=root)
{
wellMap[wellID] = std::make_tuple(well.begin(),well.end(),well.size());
wellMap[wellID] = std::make_tuple(well.begin(), well.end(), well.size());
}
}
else
{
wellMap[*well.begin()] = std::make_tuple(well.begin(),well.end(),well.size());
wellMap[*well.begin()] = std::make_tuple(well.begin(), well.end(), well.size());
}
}

Expand All @@ -253,7 +253,7 @@ extendRootExportList(const GraphOfGrid<Dune::CpGrid>& gog,
int rankToExport = std::get<1>(cellProperties);
if (rankToExport!=root)
{
const auto& [begin,end,wSize] = pWell->second;
const auto& [begin, end, wSize] = pWell->second;
std::vector<int> wellToExport;
wellToExport.reserve(wSize);
wellToExport.push_back(*begin);
Expand All @@ -279,12 +279,12 @@ extendRootExportList(const GraphOfGrid<Dune::CpGrid>& gog,
}

// add new cells to the exportList and sort it. It is assumed that exportList starts sorted.
std::sort(addToList.begin(),addToList.end());
std::sort(addToList.begin(), addToList.end());
auto origSize = exportList.size();
auto totsize = origSize+addToList.size();
exportList.reserve(totsize);
exportList.insert(exportList.end(),addToList.begin(),addToList.end());
std::inplace_merge(exportList.begin(),exportList.begin()+origSize,exportList.end());
exportList.insert(exportList.end(), addToList.begin(), addToList.end());
std::inplace_merge(exportList.begin(), exportList.begin()+origSize, exportList.end());

return exportedWells;
}
Expand Down Expand Up @@ -343,7 +343,7 @@ std::vector<std::vector<int>> communicateExportedWells(
int wellSize = receivedData[index++];
assert(index+wellSize<=totsize);
const auto dataBegin = receivedData.begin()+index;
result[i] = std::vector<int>(dataBegin,dataBegin+wellSize);
result[i] = std::vector<int>(dataBegin, dataBegin+wellSize);
index+=wellSize;
}
}
Expand All @@ -355,7 +355,7 @@ void extendImportList(std::vector<std::tuple<int,int,char,int>>& importList,
{
using ImportList = std::vector<std::tuple<int,int,char,int>>;
// make a list of wells for easy identification
std::unordered_map<int,std::size_t> wellMap;
std::unordered_map<int, std::size_t> wellMap;
for (std::size_t i=0; i<extraWells.size(); ++i)
{
if (extraWells[i].size()>1)
Expand Down Expand Up @@ -392,12 +392,12 @@ void extendImportList(std::vector<std::tuple<int,int,char,int>>& importList,
}

// add new cells to the importList and sort it. It is assumed that importList starts sorted.
std::sort(addToList.begin(),addToList.end());
std::sort(addToList.begin(), addToList.end());
auto origSize = importList.size();
auto totsize = origSize+addToList.size();
importList.reserve(totsize);
importList.insert(importList.end(),addToList.begin(),addToList.end());
std::inplace_merge(importList.begin(),importList.begin()+origSize,importList.end());
importList.insert(importList.end(), addToList.begin(), addToList.end());
std::inplace_merge(importList.begin(), importList.begin()+origSize, importList.end());
}

} // end namespace Impl
Expand All @@ -410,13 +410,13 @@ void extendExportAndImportLists(const GraphOfGrid<Dune::CpGrid>& gog,
const std::vector<int>& gIDtoRank)
{
// extend root's export list and get sets of well cells for other ranks
auto expListToComm = Impl::extendRootExportList(gog,exportList,root,gIDtoRank);
auto expListToComm = Impl::extendRootExportList(gog, exportList, root, gIDtoRank);
// obtain wells on this rank from root
auto extraWells = Impl::communicateExportedWells(expListToComm,cc,root);
auto extraWells = Impl::communicateExportedWells(expListToComm, cc, root);
if (cc.rank()!=root)
{
std::sort(importList.begin(),importList.end());
Impl::extendImportList(importList,extraWells);
std::sort(importList.begin(), importList.end());
Impl::extendImportList(importList, extraWells);
}
}
#endif // HAVE_MPI
Expand All @@ -434,7 +434,7 @@ std::vector<int> getWellRanks(const std::vector<int>& gIDtoRank,
}

#if HAVE_MPI
std::vector<std::pair<std::string,bool>>
std::vector<std::pair<std::string, bool>>
wellsOnThisRank(const std::vector<Dune::cpgrid::OpmWellType>& wells,
const std::vector<int>& wellRanks,
const Dune::cpgrid::CpGridDataTraits::Communication& cc,
Expand All @@ -451,7 +451,7 @@ wellsOnThisRank(const std::vector<Dune::cpgrid::OpmWellType>& wells,

template<class Id>
std::tuple<std::vector<int>,
std::vector<std::pair<std::string,bool>>,
std::vector<std::pair<std::string, bool>>,
std::vector<std::tuple<int,int,char> >,
std::vector<std::tuple<int,int,char,int> > >
makeImportAndExportLists(const GraphOfGrid<Dune::CpGrid>& gog,
Expand Down Expand Up @@ -488,7 +488,7 @@ makeImportAndExportLists(const GraphOfGrid<Dune::CpGrid>& gog,

for ( int i=0; i < numImport; ++i )
{
myImportList.emplace_back(importGlobalGids[i], root, static_cast<char>(AttributeSet::owner),-1);
myImportList.emplace_back(importGlobalGids[i], root, static_cast<char>(AttributeSet::owner), -1);
}
assert(rank==root || numExport==0);
if (rank==root)
Expand All @@ -498,9 +498,9 @@ makeImportAndExportLists(const GraphOfGrid<Dune::CpGrid>& gog,
gIDtoRank[exportGlobalGids[i]] = exportToPart[i];
myExportList.emplace_back(exportGlobalGids[i], exportToPart[i], static_cast<char>(AttributeSet::owner));
}
std::sort(myExportList.begin(),myExportList.end());
std::sort(myExportList.begin(), myExportList.end());
// partitioner sees only one cell per well, modify remaining
extendGIDtoRank(gog,gIDtoRank,rank);
extendGIDtoRank(gog, gIDtoRank, rank);

// Add cells that stay here to the lists. Somehow I could not persuade Zoltan to do this.
// This also adds all well cells that were missing in the importGlobalIDs.
Expand All @@ -517,11 +517,11 @@ makeImportAndExportLists(const GraphOfGrid<Dune::CpGrid>& gog,
}


std::vector<std::pair<std::string,bool>> parallel_wells;
std::vector<std::pair<std::string, bool>> parallel_wells;
if( wells )
{
// complete root's export and other's import list by adding remaining well cells
extendExportAndImportLists(gog,cc,root,myExportList,myImportList,gIDtoRank);
extendExportAndImportLists(gog, cc, root, myExportList, myImportList, gIDtoRank);

auto wellRanks = getWellRanks(gIDtoRank, wellConnections);
parallel_wells = wellsOnThisRank(*wells, wellRanks, cc, root);
Expand All @@ -540,7 +540,7 @@ void setDefaultZoltanParameters(Zoltan_Struct* zz)
Zoltan_Set_Param(zz, "NUM_GID_ENTRIES", "1");
Zoltan_Set_Param(zz, "NUM_LID_ENTRIES", "0");
Zoltan_Set_Param(zz, "RETURN_LISTS", "ALL");
Zoltan_Set_Param(zz, "EDGE_WEIGHT_DIM","1");
Zoltan_Set_Param(zz, "EDGE_WEIGHT_DIM", "1");
Zoltan_Set_Param(zz, "OBJ_WEIGHT_DIM", "1");
Zoltan_Set_Param(zz, "PHG_EDGE_SIZE_THRESHOLD", ".35"); /* 0-remove all, 1-remove none */
Zoltan_Set_Param(zz, "DEBUG_LEVEL", "0");
Expand All @@ -553,7 +553,7 @@ void setDefaultZoltanParameters(Zoltan_Struct* zz)

} // anon namespace

std::tuple<std::vector<int>, std::vector<std::pair<std::string,bool>>,
std::tuple<std::vector<int>, std::vector<std::pair<std::string, bool>>,
std::vector<std::tuple<int,int,char> >,
std::vector<std::tuple<int,int,char,int> >,
Dune::cpgrid::WellConnections>
Expand All @@ -565,7 +565,7 @@ zoltanPartitioningWithGraphOfGrid(const Dune::CpGrid& grid,
[[maybe_unused]] Dune::EdgeWeightMethod edgeWeightsMethod,
int root,
const double zoltanImbalanceTol,
const std::map<std::string,std::string>& params)
const std::map<std::string, std::string>& params)
{
int rc = ZOLTAN_OK - 1;
float ver = 0;
Expand All @@ -591,14 +591,14 @@ zoltanPartitioningWithGraphOfGrid(const Dune::CpGrid& grid,

// prepare graph and contract well cells
// non-root processes have empty grid and no wells
GraphOfGrid gog(grid);
GraphOfGrid gog(grid, transmissibilities);
assert(gog.size()==0 || !partitionIsEmpty);
auto wellConnections=partitionIsEmpty ? Dune::cpgrid::WellConnections()
: Dune::cpgrid::WellConnections(*wells,possibleFutureConnections,grid);
addWellConnections(gog,wellConnections);
: Dune::cpgrid::WellConnections(*wells, possibleFutureConnections, grid);
addWellConnections(gog, wellConnections);

// call partitioner
setGraphOfGridZoltanGraphFunctions(zz,gog,partitionIsEmpty);
setGraphOfGridZoltanGraphFunctions(zz, gog, 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 */
Expand Down
6 changes: 3 additions & 3 deletions opm/grid/GraphOfGridWrappers.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -209,7 +209,7 @@ std::vector<int> getWellRanks(const std::vector<int>& gIDtoRank,
///
/// This function only gets the information from wellRanks into proper
/// format to call computeParallelWells.
std::vector<std::pair<std::string,bool>>
std::vector<std::pair<std::string, bool>>
wellsOnThisRank(const std::vector<Dune::cpgrid::OpmWellType>& wells,
const std::vector<int>& wellRanks,
const Dune::cpgrid::CpGridDataTraits::Communication& cc,
Expand All @@ -235,7 +235,7 @@ wellsOnThisRank(const std::vector<Dune::cpgrid::OpmWellType>& wells,
/// myImportList vector of cells to be moved to this rank
template<class Id>
std::tuple<std::vector<int>,
std::vector<std::pair<std::string,bool>>,
std::vector<std::pair<std::string, bool>>,
std::vector<std::tuple<int,int,char> >,
std::vector<std::tuple<int,int,char,int> > >
makeImportAndExportLists(const GraphOfGrid<Dune::CpGrid>& gog,
Expand All @@ -255,7 +255,7 @@ makeImportAndExportLists(const GraphOfGrid<Dune::CpGrid>& gog,
/// GraphOfGrid represents a well by one vertex, so wells can not be
/// spread over several processes.
/// transmissiblities are currently not supported, but are queued
std::tuple<std::vector<int>, std::vector<std::pair<std::string,bool>>,
std::tuple<std::vector<int>, std::vector<std::pair<std::string, bool>>,
std::vector<std::tuple<int,int,char> >,
std::vector<std::tuple<int,int,char,int> >,
Dune::cpgrid::WellConnections>
Expand Down
Loading

0 comments on commit 73423a7

Please sign in to comment.