OilfieldToolsNet/opm-simulators
4
0 Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2026-10-26 15:40:09 -05:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #6103 from jakobtorben/add_support_for_merging_well_cells_for_nldd_partitioning

Browse Source 
Add support for merging well cells for NLDD partitioning
This commit is contained in:
Bård Skaflestad authored and GitHub committed 2025-03-28 09:59:24 +01:00
2 parents 43697e9a7c 36cda5abd6
commit 97175336ab
8 files changed
+860 -94
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

opm/simulators/flow/BlackoilModelNldd.hpp
+66 -11
View File
@@ -107,7 +107,26 @@ public:
, rank_(model_.simulator().vanguard().grid().comm().rank())
{
// Create partitions.
const auto& [partition_vector, num_domains] = this->partitionCells();
const auto& [partition_vector_initial, num_domains_initial] = this->partitionCells();
int num_domains = num_domains_initial;
std::vector<int> partition_vector = partition_vector_initial;
// Fix-up for an extreme case: Interior cells who do not have any on-rank
// neighbours. Move all such cells into a single domain on this rank,
// and mark the domain for skipping. For what it's worth, we've seen this
// case occur in practice when testing on field cases.
bool isolated_cells = false;
for (auto& domainId : partition_vector) {
if (domainId < 0) {
domainId = num_domains;
domains_to_skip_.push_back(domainId);
isolated_cells = true;
}
}
if (isolated_cells) {
num_domains++;
}
// Set nldd handler in main well model
model.wellModel().setNlddAdapter(&wellModel_);
@@ -222,6 +241,9 @@ public:
for (const int domain_index : domain_order) {
const auto& domain = domains_[domain_index];
SimulatorReportSingle local_report;
if (std::find(domains_to_skip_.begin(), domains_to_skip_.end(), domain.index) != domains_to_skip_.end()) {
continue;
}
try {
switch (model_.param().local_solve_approach_) {
case DomainSolveApproach::Jacobi:
@@ -340,23 +362,55 @@ public:
return local_reports_accumulated_;
}
/// Write the partition vector to a file in ResInsight compatible format for inspection
/// and a partition file for each rank, that can be used as input for OPM.
void writePartitions(const std::filesystem::path& odir) const
{
const auto& elementMapper = this->model_.simulator().model().elementMapper();
const auto& cartMapper = this->model_.simulator().vanguard().cartesianIndexMapper();
const auto& dims = cartMapper.cartesianDimensions();
const auto total_size = dims[0] * dims[1] * dims[2];
const auto& grid = this->model_.simulator().vanguard().grid();
const auto& comm = grid.comm();
const auto nDigit = 1 + static_cast<int>(std::floor(std::log10(comm.size())));
std::ofstream pfile { odir / fmt::format("{1:0>{0}}", nDigit, comm.rank()) };
// Create a full-sized vector initialized with -1 (indicating inactive cells)
auto full_partition = std::vector<int>(total_size, -1);
// Fill active cell values for this rank
const auto p = this->reconstitutePartitionVector();
auto i = 0;
for (const auto& cell : elements(grid.leafGridView(), Dune::Partitions::interior)) {
full_partition[cartMapper.cartesianIndex(elementMapper.index(cell))] = p[i++];
}
// Gather all partitions using max operation
comm.max(full_partition.data(), full_partition.size());
// Only rank 0 writes the file
if (this->rank_ == 0) {
auto fname = odir / "ResInsight_compatible_partition.txt";
std::ofstream resInsightFile { fname };
// Write header
resInsightFile << "NLDD_DOM" << '\n';
// Write all cells, including inactive ones
for (const auto& val : full_partition) {
resInsightFile << val << '\n';
}
resInsightFile << "/" << '\n';
}
const auto nDigit = 1 + static_cast<int>(std::floor(std::log10(comm.size())));
auto partition_fname = odir / fmt::format("{1:0>{0}}", nDigit, comm.rank());
std::ofstream pfile { partition_fname };
auto cell_index = 0;
for (const auto& cell : elements(grid.leafGridView(), Dune::Partitions::interior)) {
pfile << comm.rank() << ' '
<< cartMapper.cartesianIndex(elementMapper.index(cell)) << ' '
<< p[i++] << '\n';
<< p[cell_index++] << '\n';
}
}
@@ -763,7 +817,7 @@ private:
} else if (model_.param().local_solve_approach_ == DomainSolveApproach::GaussSeidel) {
// Calculate the measure used to order the domains.
std::vector<Scalar> measure_per_domain(domains_.size());
switch (model_.param().local_domain_ordering_) {
switch (model_.param().local_domains_ordering_) {
case DomainOrderingMeasure::AveragePressure: {
// Use average pressures to order domains.
for (const auto& domain : domains_) {
@@ -802,7 +856,7 @@ private:
}
break;
}
} // end of switch (model_.param().local_domain_ordering_)
} // end of switch (model_.param().local_domains_ordering_)
// Sort by largest measure, keeping index order if equal.
const auto& m = measure_per_domain;
@@ -937,7 +991,7 @@ private:
auto zoltan_ctrl = ZoltanPartitioningControl<Element>{};
zoltan_ctrl.domain_imbalance = param.local_domain_partition_imbalance_;
zoltan_ctrl.domain_imbalance = param.local_domains_partition_imbalance_;
zoltan_ctrl.index =
[elementMapper = &this->model_.simulator().model().elementMapper()]
@@ -954,7 +1008,7 @@ private:
};
// Forming the list of wells is expensive, so do this only if needed.
const auto need_wells = param.local_domain_partition_method_ == "zoltan";
const auto need_wells = param.local_domains_partition_method_ == "zoltan";
const auto wells = need_wells
? this->model_.simulator().vanguard().schedule().getWellsatEnd()
@@ -969,10 +1023,10 @@ private:
const int num_domains = (param.num_local_domains_ > 0)
? param.num_local_domains_
: detail::countGlobalCells(grid) / default_cells_per_domain;
return ::Opm::partitionCells(param.local_domain_partition_method_,
return ::Opm::partitionCells(param.local_domains_partition_method_,
num_domains, grid.leafGridView(), wells,
possibleFutureConnectionSet, zoltan_ctrl);
possibleFutureConnectionSet, zoltan_ctrl,
param.local_domains_partition_well_neighbor_levels_);
}
std::vector<int> reconstitutePartitionVector() const
@@ -1008,6 +1062,7 @@ private:
std::vector<Domain> domains_; //!< Vector of subdomains
std::vector<std::unique_ptr<Mat>> domain_matrices_; //!< Vector of matrix operator for each subdomain
std::vector<ISTLSolverType> domain_linsolvers_; //!< Vector of linear solvers for each domain
std::vector<int> domains_to_skip_; //!< Vector of domains to skip
SimulatorReportSingle local_reports_accumulated_; //!< Accumulated convergence report for subdomain solvers
int rank_ = 0; //!< MPI rank of this process
};
opm/simulators/flow/BlackoilModelParameters.cpp
+6 -3
View File
@@ -93,15 +93,16 @@ BlackoilModelParameters<Scalar>::BlackoilModelParameters()
local_tolerance_scaling_cnv_ = Parameters::Get<Parameters::LocalToleranceScalingCnv<Scalar>>();
nldd_num_initial_newton_iter_ = Parameters::Get<Parameters::NlddNumInitialNewtonIter>();
num_local_domains_ = Parameters::Get<Parameters::NumLocalDomains>();
local_domain_partition_imbalance_ = std::max(Scalar{1.0}, Parameters::Get<Parameters::LocalDomainsPartitioningImbalance<Scalar>>());
local_domain_partition_method_ = Parameters::Get<Parameters::LocalDomainsPartitioningMethod>();
local_domains_partition_imbalance_ = std::max(Scalar{1.0}, Parameters::Get<Parameters::LocalDomainsPartitioningImbalance<Scalar>>());
local_domains_partition_method_ = Parameters::Get<Parameters::LocalDomainsPartitioningMethod>();
local_domains_partition_well_neighbor_levels_ = Parameters::Get<Parameters::LocalDomainsPartitionWellNeighborLevels>();
deck_file_name_ = Parameters::Get<Parameters::EclDeckFileName>();
network_max_strict_outer_iterations_ = Parameters::Get<Parameters::NetworkMaxStrictOuterIterations>();
network_max_outer_iterations_ = Parameters::Get<Parameters::NetworkMaxOuterIterations>();
network_max_sub_iterations_ = Parameters::Get<Parameters::NetworkMaxSubIterations>();
network_pressure_update_damping_factor_ = Parameters::Get<Parameters::NetworkPressureUpdateDampingFactor<Scalar>>();
network_max_pressure_update_in_bars_ = Parameters::Get<Parameters::NetworkMaxPressureUpdateInBars<Scalar>>();
local_domain_ordering_ = domainOrderingMeasureFromString(Parameters::Get<Parameters::LocalDomainsOrderingMeasure>());
local_domains_ordering_ = domainOrderingMeasureFromString(Parameters::Get<Parameters::LocalDomainsOrderingMeasure>());
write_partitions_ = Parameters::Get<Parameters::DebugEmitCellPartition>();
monitor_params_.enabled_ = Parameters::Get<Parameters::ConvergenceMonitoring>();
@@ -253,6 +254,8 @@ void BlackoilModelParameters<Scalar>::registerParameters()
"'zoltan', "
"'simple', "
"and the name of a partition file ending with '.partition'.");
Parameters::Register<Parameters::LocalDomainsPartitionWellNeighborLevels>
("Number of neighbor levels around wells to include in the same domain during NLDD partitioning");
Parameters::Register<Parameters::LocalDomainsOrderingMeasure>
("Subdomain ordering measure. Allowed values are "
"'maxpressure', "
opm/simulators/flow/BlackoilModelParameters.hpp
+6 -3
View File
@@ -152,6 +152,7 @@ template<class Scalar>
struct LocalDomainsPartitioningImbalance { static constexpr Scalar value = 1.03; };
struct LocalDomainsPartitioningMethod { static constexpr auto value = "zoltan"; };
struct LocalDomainsPartitionWellNeighborLevels { static constexpr int value = 1; };
struct LocalDomainsOrderingMeasure { static constexpr auto value = "maxpressure"; };
struct ConvergenceMonitoring { static constexpr bool value = false; };
@@ -159,6 +160,7 @@ struct ConvergenceMonitoringCutOff { static constexpr int value = 6; };
template<class Scalar>
struct ConvergenceMonitoringDecayFactor { static constexpr Scalar value = 0.75; };
template<class Scalar>
struct NupcolGroupRateTolerance { static constexpr Scalar value = 0.001; };
@@ -327,9 +329,10 @@ public:
int nldd_num_initial_newton_iter_{1};
int num_local_domains_{0};
Scalar local_domain_partition_imbalance_{1.03};
std::string local_domain_partition_method_;
DomainOrderingMeasure local_domain_ordering_{DomainOrderingMeasure::MaxPressure};
Scalar local_domains_partition_imbalance_{1.03};
std::string local_domains_partition_method_;
int local_domains_partition_well_neighbor_levels_{1};
DomainOrderingMeasure local_domains_ordering_{DomainOrderingMeasure::MaxPressure};
bool write_partitions_{false};
opm/simulators/flow/partitionCells.cpp
+104 -37
View File
@@ -109,11 +109,24 @@ public:
///
/// \param[in] zoltan_ctrl Control parameters for on-rank subdomain
/// partitioning.
///
/// \param[in] num_neighbor_levels Number of neighbor levels to include when connecting well cells.
/// Default is 0, which means only direct well connections are considered.
template <class GridView, class Element>
void buildLocalGraph(const GridView& grid_view,
const std::vector<Opm::Well>& wells,
const std::unordered_map<std::string, std::set<int>>& possibleFutureConnections,
const Opm::ZoltanPartitioningControl<Element>& zoltan_ctrl);
const Opm::ZoltanPartitioningControl<Element>& zoltan_ctrl,
const int num_neighbor_levels);
/// Connect neighbors of cells up to a specified level.
///
/// \param[in,out] cells Initial cells to find neighbors of. Updated to include all neighbors
/// up to specified level.
///
/// \param[in] num_levels Number of neighbor levels to include.
void connectNeighbors(std::vector<int>& cells,
const int num_levels) const;
/// Partition rank's interior cells into non-overlapping domains using
/// the Zoltan graph partitioning software package.
@@ -148,6 +161,9 @@ private:
/// Zoltan partitioning backend.
Opm::ParallelNLDDPartitioningZoltan partitioner_;
/// Map of cell indices to their neighbors for fast lookup
std::unordered_map<int, std::unordered_set<int>> neighbor_map_;
/// Form internal connectivity graph between rank's reachable, interior
/// cells.
///
@@ -166,7 +182,8 @@ private:
template <class GridView, class Element>
std::unordered_map<int, int>
connectElements(const GridView& grid_view,
const Opm::ZoltanPartitioningControl<Element>& zoltan_ctrl);
const Opm::ZoltanPartitioningControl<Element>& zoltan_ctrl,
const bool build_neighbor_map);
/// Adapt internal connectivity graph to account for connections induced
/// by well reservoir connections.
@@ -182,11 +199,15 @@ private:
///
/// \param[in] g2l Mapping from globally unique cell IDs to local,
/// on-rank active cell IDs. Return value from \c connectElements().
///
/// \param[in] num_neighbor_levels Number of neighbor levels to include when connecting well cells.
/// Default is 0, which means only direct well connections are considered.
template <typename Comm>
void connectWells(const Comm comm,
const std::vector<Opm::Well>& wells,
std::unordered_map<std::string, std::set<int>> possibleFutureConnections,
const std::unordered_map<int, int>& g2l);
const std::unordered_map<int, int>& g2l,
const int num_neighbor_levels);
};
// Note: "grid_view.size(0)" is intentional here. It is not an error. The
@@ -203,9 +224,43 @@ template <class GridView, class Element>
void ZoltanPartitioner::buildLocalGraph(const GridView& grid_view,
const std::vector<Opm::Well>& wells,
const std::unordered_map<std::string, std::set<int>>& possibleFutureConnections,
const Opm::ZoltanPartitioningControl<Element>& zoltan_ctrl)
const Opm::ZoltanPartitioningControl<Element>& zoltan_ctrl,
const int num_neighbor_levels)
{
this->connectWells(grid_view.comm(), wells, possibleFutureConnections, this->connectElements(grid_view, zoltan_ctrl));
const bool build_neighbor_map = num_neighbor_levels > 0;
auto g2l = this->connectElements(grid_view, zoltan_ctrl, build_neighbor_map);
this->connectWells(grid_view.comm(), wells, possibleFutureConnections, g2l, num_neighbor_levels);
}
void ZoltanPartitioner::connectNeighbors(std::vector<int>& cells,
const int num_levels) const
{
// Initialize search structures
std::unordered_set<int> visited(cells.begin(), cells.end());
std::vector<int> frontier;
std::vector<int> new_frontier;
frontier.swap(cells);
// Breadth-first search for neighbor cells up to num_levels away
for (int level = 0; level < num_levels && !frontier.empty(); ++level) {
new_frontier.clear();
new_frontier.reserve(frontier.size() * 2);
// Expand current frontier using precomputed neighbor map
for (const int cell : frontier) {
for (const int neighbor : neighbor_map_.at(cell)) {
if (visited.insert(neighbor).second) { // New discovery
new_frontier.push_back(neighbor);
}
}
}
frontier.swap(new_frontier); // Prepare next level
}
// Sort final result for faster processing in later steps
cells.assign(visited.begin(), visited.end());
std::sort(cells.begin(), cells.end());
}
template <class GridView, class Element>
@@ -233,35 +288,30 @@ ZoltanPartitioner::partition(const int num_
num_domains = domainId + 1;
}
}
// Fix-up for an extreme case: Interior cells whose neighbours are all
// on a different rank--i.e., interior cells which do not have on-rank
// neighbours. Make each such cell be a separate domain on this rank.
// Don't remove this code unless you've taken remedial action elsewhere
// to ensure the situation doesn't happen. For what it's worth, we've
// seen this case occur in practice when testing on field cases.
for (auto& domainId : domains) {
if (domainId < 0) {
domainId = num_domains++;
}
}
return partition;
}
template <class GridView, class Element>
std::unordered_map<int, int>
ZoltanPartitioner::connectElements(const GridView& grid_view,
const Opm::ZoltanPartitioningControl<Element>& zoltan_ctrl)
const Opm::ZoltanPartitioningControl<Element>& zoltan_ctrl,
const bool build_neighbor_map)
{
auto g2l = std::unordered_map<int, int>{};
if (build_neighbor_map) {
neighbor_map_.clear();
}
for (const auto& element : elements(grid_view, Dune::Partitions::interior)) {
{
const auto c = zoltan_ctrl.index(element);
g2l.insert_or_assign(zoltan_ctrl.local_to_global(c), c);
}
if (build_neighbor_map) {
// Initialize neighbor set for this cell
neighbor_map_[c];
}
}
for (const auto& is : intersections(grid_view, element)) {
if (! is.neighbor()) {
continue;
@@ -282,6 +332,12 @@ ZoltanPartitioner::connectElements(const GridView&
this->partitioner_.registerConnection(c1, c2);
this->partitioner_.registerConnection(c2, c1);
if (build_neighbor_map) {
// Store neighbors in both directions
neighbor_map_[c1].insert(c2);
neighbor_map_[c2].insert(c1);
}
}
}
@@ -292,7 +348,8 @@ template <typename Comm>
void ZoltanPartitioner::connectWells(const Comm comm,
const std::vector<Opm::Well>& wells,
std::unordered_map<std::string, std::set<int>> possibleFutureConnections,
const std::unordered_map<int, int>& g2l)
const std::unordered_map<int, int>& g2l,
const int num_neighbor_levels)
{
auto distributedWells = 0;
@@ -310,7 +367,6 @@ void ZoltanPartitioner::connectWells(const Comm
++otherProc;
continue;
}
cellIx.push_back(locPos->second);
}
const auto possibleFutureConnectionSetIt = possibleFutureConnections.find(well.name());
@@ -330,17 +386,23 @@ void ZoltanPartitioner::connectWells(const Comm
continue;
}
const auto nc = cellIx.size();
for (auto ic1 = 0*nc; ic1 < nc; ++ic1) {
for (auto ic2 = ic1 + 1; ic2 < nc; ++ic2) {
this->partitioner_.registerConnection(cellIx[ic1], cellIx[ic2]);
this->partitioner_.registerConnection(cellIx[ic2], cellIx[ic1]);
}
// Process neighbors and get deduplicated result
if (num_neighbor_levels > 0 && !cellIx.empty()) {
connectNeighbors(cellIx, num_neighbor_levels);
}
if (! cellIx.empty()) {
// All cells intersected by a single well go in the same domain.
this->partitioner_.forceSameDomain(std::move(cellIx));
const auto nc = cellIx.size();
if (nc == 0) {
continue;
}
for (auto ic1 = 0*nc; ic1 < nc-1; ++ic1) {
this->partitioner_.registerConnection(cellIx[ic1], cellIx[ic1+1]);
this->partitioner_.registerConnection(cellIx[ic1+1], cellIx[ic1]);
}
// Collect all cells for the well into a vertex group
if (!cellIx.empty()) {
this->partitioner_.addVertexGroup(cellIx);
}
}
@@ -356,7 +418,8 @@ This is not supported in the current NLDD implementation.)",
};
}
OPM_END_PARALLEL_TRY_CATCH(std::string { "ZoltanPartitioner::connectWells()" }, comm)
OPM_END_PARALLEL_TRY_CATCH(std::string{"ZoltanPartitioner::connectWells(): Distributed well detected. "
"This is not supported in the current NLDD implementation"}, comm)
}
template <class GridView, class Element>
@@ -365,7 +428,8 @@ partitionCellsZoltan(const int num_d
const GridView& grid_view,
const std::vector<Opm::Well>& wells,
const std::unordered_map<std::string, std::set<int>>& possibleFutureConnections,
const Opm::ZoltanPartitioningControl<Element>& zoltan_ctrl)
const Opm::ZoltanPartitioningControl<Element>& zoltan_ctrl,
const int num_neighbor_levels)
{
if (num_domains <= 1) { // No partitioning => every cell in domain zero.
const auto num_interior_cells =
@@ -379,7 +443,8 @@ partitionCellsZoltan(const int num_d
}
auto partitioner = ZoltanPartitioner { grid_view, zoltan_ctrl.local_to_global };
partitioner.buildLocalGraph(grid_view, wells, possibleFutureConnections, zoltan_ctrl);
partitioner.buildLocalGraph(grid_view, wells, possibleFutureConnections, zoltan_ctrl, num_neighbor_levels);
return partitioner.partition(num_domains, grid_view, zoltan_ctrl);
}
@@ -600,12 +665,13 @@ Opm::partitionCells(const std::string& method,
const GridView& grid_view,
[[maybe_unused]] const std::vector<Well>& wells,
[[maybe_unused]] const std::unordered_map<std::string, std::set<int>>& possibleFutureConnections,
[[maybe_unused]] const ZoltanPartitioningControl<Element>& zoltan_ctrl)
[[maybe_unused]] const ZoltanPartitioningControl<Element>& zoltan_ctrl,
[[maybe_unused]] const int num_neighbor_levels)
{
if (method == "zoltan") {
#if HAVE_MPI && HAVE_ZOLTAN
return partitionCellsZoltan(num_local_domains, grid_view, wells, possibleFutureConnections, zoltan_ctrl);
return partitionCellsZoltan(num_local_domains, grid_view, wells, possibleFutureConnections, zoltan_ctrl, num_neighbor_levels);
#else // !HAVE_MPI || !HAVE_ZOLTAN
@@ -700,7 +766,8 @@ Opm::partitionCellsSimple(const int num_cells, const int num_domains)
const std::unordered_map<std::string, std::set<int>>&, \
const Opm::ZoltanPartitioningControl< \
typename std::remove_cv_t<std::remove_reference_t< \
decltype(std::declval<Grid>().leafGridView())>>::template Codim<0>::Entity>&)
decltype(std::declval<Grid>().leafGridView())>>::template Codim<0>::Entity>&, \
const int)
// ---------------------------------------------------------------------------
// Grid types built into Flow.
opm/simulators/flow/partitionCells.hpp
+4 -1
View File
@@ -84,6 +84,8 @@ struct ZoltanPartitioningControl
/// \param[in] zoltan_ctrl Control parameters for local Zoltan-based
/// partitioning. Not used unless \code method == "zoltan" \endcode.
///
/// \param[in] num_neighbor_levels Number of neighbor levels to consider for partitioning.
///
/// \return Partition vector--subdomain ID for each cell in \p grid_view
/// traversal order for its interior cells--and the number of subdomains
/// on current rank.
@@ -94,7 +96,8 @@ partitionCells(const std::string& method,
const GridView& grid_view,
const std::vector<Well>& wells,
const std::unordered_map<std::string, std::set<int>>& possibleFutureConnections,
const ZoltanPartitioningControl<Element>& zoltan_ctrl);
const ZoltanPartitioningControl<Element>& zoltan_ctrl,
const int num_neighbor_levels = 0);
/// Read a partitioning from file, assumed to contain one number per cell, its partition number.
/// \return pair containing a partition vector (partition number for each cell), and the number of partitions.
opm/simulators/utils/ParallelNLDDPartitioningZoltan.cpp
+46 -18
View File
@@ -42,6 +42,8 @@
#include <stdexcept>
#include <utility>
#include <vector>
#include <set>
#include <map>
namespace {
@@ -155,6 +157,8 @@ namespace {
///
/// \param[in] vertexId Enumeration of reachable vertices.
///
/// \param[in] vertices_to_merge Map of vertices to be merged.
///
/// \param[in] globalCell Callback for mapping (local) vertex IDs to
/// globally unique vertex IDs.
template <typename Edge, typename GlobalCellID>
@@ -162,6 +166,7 @@ namespace {
const std::size_t numVertices,
const std::vector<Edge>& edges,
const EnumerateSeenVertices& vertexId,
const std::vector<std::vector<int>>& vertexGroups,
GlobalCellID&& globalCell)
: myRank_ { myRank }
{
@@ -171,7 +176,14 @@ namespace {
this->graph_.addConnection(vertexId[v2], vertexId[v1]);
}
this->graph_.compress(vertexId.numVertices());
// Apply any vertex merges that were collected
for (const auto& group : vertexGroups) {
this->graph_.addVertexGroup(group);
}
auto numVerticesAfterMerge = this->graph_.applyVertexMerges();
this->graph_.compress(numVerticesAfterMerge);
// Form local-to-global vertex ID mapping for reachable vertices.
this->globalCell_.resize(vertexId.numVertices());
@@ -220,6 +232,15 @@ namespace {
return this->graph_.columnIndices();
}
/// Get the final vertex ID after merging and renumbering for a given original vertex ID.
///
/// \param[in] originalVertexID The original vertex ID to look up
/// \return The final vertex ID after merges and renumbering
int getFinalVertexID(const int originalVertexID) const
{
return this->graph_.getFinalVertexID(originalVertexID);
}
private:
// VertexID = int, TrackCompressedIdx = false
using Backend = Opm::utility::CSRGraphFromCoordinates<>;
@@ -536,27 +557,36 @@ extern "C" {
return blocks;
}
void forceSameDomain(const std::vector<int>& cells,
std::vector<int>& parts)
{
if (cells.empty()) { return; }
const auto first = parts[cells.front()];
for (const auto& cell : cells) {
parts[cell] = first;
}
}
} // Anonymous namespace
void Opm::ParallelNLDDPartitioningZoltan::addVertexGroup(const std::vector<int>& vertices)
{
vertexGroups_.push_back(vertices);
}
std::vector<int>
Opm::ParallelNLDDPartitioningZoltan::partitionElements(const ZoltanParamMap& params) const
{
const auto vertexId = EnumerateSeenVertices { this->numElements_, this->conns_ };
// The vertexGraph is built using the reachable indices, so for the vertexGroups_
// we need to map the original local indices to the reachable indices
std::vector<std::vector<int>> reachableVertexGroups;
for (const auto& group : vertexGroups_) {
std::vector<int> reachableGroup;
for (const int originalLocalIdx : group) {
if (const int reachableIdx = vertexId[originalLocalIdx]; reachableIdx != -1) {
reachableGroup.push_back(reachableIdx);
}
}
if (!reachableGroup.empty()) {
reachableVertexGroups.push_back(reachableGroup);
}
}
auto graph = VertexGraph {
this->comm_.rank(), this->numElements_,
this->conns_, vertexId, this->globalCell_
this->conns_, vertexId, reachableVertexGroups, this->globalCell_
};
const auto partsForReachableCells = Partitioner {
@@ -567,13 +597,11 @@ Opm::ParallelNLDDPartitioningZoltan::partitionElements(const ZoltanParamMap& par
auto parts = std::vector<int>(this->numElements_, -1);
for (auto elm = 0*this->numElements_; elm < this->numElements_; ++elm) {
if (const auto reachableElmIx = vertexId[elm]; reachableElmIx >= 0) {
parts[elm] = partsForReachableCells[reachableElmIx];
// Get the final vertex ID after merging and renumbering
const auto finalVertexId = graph.getFinalVertexID(reachableElmIx);
parts[elm] = partsForReachableCells[finalVertexId];
}
}
for (const auto& cells : this->sameDomain_) {
::forceSameDomain(cells, parts);
}
return util::compressPartitionIDs(std::move(parts));
}
opm/simulators/utils/ParallelNLDDPartitioningZoltan.hpp
+7 -18
View File
@@ -72,18 +72,6 @@ namespace Opm {
this->conns_.emplace_back(c1, c2);
}
/// Force collection of cells to be in same result domain.
///
/// Mostly as a means to ensuring wells do not intersect multiple
/// domains/blocks.
///
/// \param[in] cells Cell collection. Typically those cells which
/// are intersected by a single well.
void forceSameDomain(std::vector<int>&& cells)
{
this->sameDomain_.emplace_back(std::move(cells));
}
/// Partition connectivity graph using Zoltan graph partitioning
/// package.
///
@@ -99,9 +87,16 @@ namespace Opm {
std::vector<int>
partitionElements(const ZoltanParamMap& params) const;
/// Add a group of vertices that should be merged together.
/// Must be called before compress().
///
/// \param[in] vertices Vector of vertex IDs to merge
void addVertexGroup(const std::vector<int>& vertices);
private:
/// Connection/graph edge.
using Connection = std::pair<std::size_t, std::size_t>;
std::vector<std::vector<int>> vertexGroups_{};
/// MPI communication object. Needed by Zoltan.
Parallel::Communication comm_{};
@@ -116,12 +111,6 @@ namespace Opm {
/// Connectivity graph edges.
std::vector<Connection> conns_{};
/// Collections of vertices/cells which must be coalesced to the
/// same domain/block. All vertices within a single collection will
/// be placed on the same domain, but cells from different
/// collections may be placed on different domains.
std::vector<std::vector<int>> sameDomain_{};
};
} // namespace Opm
tests/test_partitionCells.cpp
+621 -3
View File
@@ -1,5 +1,6 @@
/*
Copyright 2021 SINTEF Digital, Mathematics and Cybernetics.
Copyright 2021,2025 SINTEF Digital, Mathematics and Cybernetics.
Copyright 2025 Equinor ASA
This file is part of the Open Porous Media project (OPM).
@@ -19,11 +20,84 @@
#include <config.h>
#define BOOST_TEST_MODULE OPM_test_aspinPartition
#define BOOST_TEST_MODULE OPM_test_partitionCells
#define BOOST_TEST_NO_MAIN
#include <boost/test/unit_test.hpp>
#include <dune/grid/common/gridview.hh>
#include <opm/grid/common/WellConnections.hpp>
#include <opm/grid/CpGrid.hpp>
#include <opm/grid/utility/OpmWellType.hpp>
#include <opm/input/eclipse/Schedule/Well/Connection.hpp>
#include <opm/input/eclipse/Schedule/Well/Well.hpp>
#include <opm/input/eclipse/Schedule/Well/WellConnections.hpp>
#include <opm/input/eclipse/Parser/Parser.hpp>
#include <opm/input/eclipse/Deck/Deck.hpp>
#include <opm/input/eclipse/EclipseState/Grid/EclipseGrid.hpp>
#include <opm/simulators/flow/partitionCells.hpp>
// Helper functions
namespace {
Opm::Connection createConnection(int i, int j, int k, int global_index) {
return Opm::Connection(i, j, k, // i,j,k
global_index, // global_index
0, // complnum
Opm::Connection::State::OPEN, // state
Opm::Connection::Direction::Z, // direction
Opm::Connection::CTFKind::DeckValue, // kind
0, // satTableId
0.0, // depth
Opm::Connection::CTFProperties(), // properties
0, // sort_value
false); // defaultSatTabId
}
Dune::cpgrid::OpmWellType createWell(const std::string& name) {
using namespace Opm;
return Dune::cpgrid::OpmWellType(name, name, 0, 0, 0, 0, 0.0, WellType(),
Well::ProducerCMode(), Connection::Order::TRACK,
UnitSystem::newMETRIC(),
0.0, 0.0, false, false, 0, Well::GasInflowEquation());
}
std::vector<Dune::cpgrid::OpmWellType> createWellsWithConnections(
const std::vector<std::pair<std::string, std::vector<int>>>& well_specs) {
std::vector<Dune::cpgrid::OpmWellType> wells;
for (const auto& [well_name, cell_indices] : well_specs) {
auto well_conn = std::make_shared<Opm::WellConnections>();
for (int idx : cell_indices) {
well_conn->add(createConnection(idx, 0, 0, idx));
}
auto well = createWell(well_name);
well.updateConnections(well_conn, true);
wells.push_back(well);
}
return wells;
}
Dune::CpGrid createTestGrid(const std::array<int, 3>& dims,
const std::array<double, 3>& size) {
Dune::CpGrid grid;
grid.createCartesian(dims, size);
return grid;
}
template<typename Entity>
Opm::ZoltanPartitioningControl<Entity> createZoltanControl(const Dune::CpGrid& grid) {
Opm::ZoltanPartitioningControl<Entity> zoltan_ctrl;
zoltan_ctrl.domain_imbalance = 1.1;
zoltan_ctrl.index = [](const auto& element) {
return element.index();
};
zoltan_ctrl.local_to_global = [&grid](const int local_idx) {
return grid.globalCell()[local_idx];
};
return zoltan_ctrl;
}
}
BOOST_AUTO_TEST_CASE(FileBased)
{
@@ -54,6 +128,550 @@ BOOST_AUTO_TEST_CASE(Simple2)
BOOST_CHECK_EQUAL(num_part, 7);
std::vector<int> expected = { 0, 0, 1, 1, 2, 2, 3, 4, 5, 6 };
BOOST_CHECK_EQUAL_COLLECTIONS(expected.begin(), expected.end(), part.begin(), part.end());
}
BOOST_AUTO_TEST_CASE(PartitionCellsTest)
{
auto grid = createTestGrid({3, 4, 1}, {3.0, 4.0, 1.0});
using Entity = typename Dune::CpGrid::LeafGridView::template Codim<0>::Entity;
auto zoltan_ctrl = createZoltanControl<Entity>(grid);
// Test simple partitioning
{
auto [part, num_part] = Opm::partitionCells("simple", 3, grid.leafGridView(),
std::vector<Opm::Well>{},
std::unordered_map<std::string, std::set<int>>{},
zoltan_ctrl);
BOOST_CHECK_EQUAL(num_part, 3);
BOOST_CHECK_EQUAL(part.size(), 12);
// Check that all partition numbers are valid
for (const auto& p : part) {
BOOST_CHECK(p >= 0 && p < 3);
}
}
#if HAVE_MPI && HAVE_ZOLTAN
// Test Zoltan partitioning
{
auto [part, num_part] = Opm::partitionCells("zoltan", 3, grid.leafGridView(),
std::vector<Opm::Well>{},
std::unordered_map<std::string, std::set<int>>{},
zoltan_ctrl);
BOOST_CHECK_EQUAL(num_part, 3);
BOOST_CHECK_EQUAL(part.size(), 12);
// Check that all partition numbers are valid
for (const auto& p : part) {
BOOST_CHECK(p >= 0 && p < 3);
}
}
#endif
}
BOOST_AUTO_TEST_CASE(PartitionCellsWithWellMergeTest)
{
// Create a 1D grid with 7 cells
// Visual representation of the grid and wells:
// Cell indices: 0 1 2 3 4 5 6
// Well layout: W1---W1 | W2---W2 | |
// Expected: Cell 0 and 1 should be merged
// Cell 3 and 4 should be merged
// Cells 2, 5, and 6 can be in any partition
auto grid = createTestGrid({7, 1, 1}, {7.0, 1.0, 1.0});
// Setup Zoltan control parameters
using Entity = typename Dune::CpGrid::LeafGridView::template Codim<0>::Entity;
auto zoltan_ctrl = createZoltanControl<Entity>(grid);
auto wells = createWellsWithConnections({
{"TESTW1", {0, 1}}, // Well 1 connects cells 0,1
{"TESTW2", {3, 4}} // Well 2 connects cells 3,4
});
// Create well connections object
Dune::cpgrid::WellConnections wellConnections(wells, std::unordered_map<std::string, std::set<int>>{}, grid);
#if HAVE_MPI && HAVE_ZOLTAN
// Test Zoltan partitioning with well cells
auto [part, num_part] = Opm::partitionCells("zoltan", 4, grid.leafGridView(),
wells,
std::unordered_map<std::string, std::set<int>>{},
zoltan_ctrl);
// Verify number of partitions
BOOST_CHECK_EQUAL(num_part, 4);
BOOST_CHECK_EQUAL(part.size(), 7);
// The key tests:
// 1. Cells 0 and 1 (well1) should be in the same partition
BOOST_CHECK_EQUAL(part[0], part[1]);
// 2. Cells 3 and 4 (well2) should be in the same partition
BOOST_CHECK_EQUAL(part[3], part[4]);
// 3. Cells 2, 5, and 6 can be in any partition, but well1 and well2 cells
// don't need to be in the same partition
#endif
}
BOOST_AUTO_TEST_CASE(PartitionCellsWithOverlappingWellsTest)
{
// Create a 1D grid with 7 cells
// Visual representation of the grid and wells:
// Cell indices: 0 1 2 3 4 5 6
// Well layout: W1---W1---W1 | | |
// W2---W2
// Expected: Cells 1,2,3 should be in same partition (Well 1)
// Cells 3,4 should be in same partition (Well 2)
// Therefore cells 1,2,3,4 should all end up in same partition
auto grid = createTestGrid({7, 1, 1}, {7.0, 1.0, 1.0});
using Entity = typename Dune::CpGrid::LeafGridView::template Codim<0>::Entity;
auto zoltan_ctrl = createZoltanControl<Entity>(grid);
auto wells = createWellsWithConnections({
{"TESTW1", {1, 2, 3}}, // Well 1 connects cells 1,2,3
{"TESTW2", {3, 4}} // Well 2 connects cells 3,4
});
#if HAVE_MPI && HAVE_ZOLTAN
// Test Zoltan partitioning with overlapping well cells
auto [part, num_part] = Opm::partitionCells("zoltan", 4, grid.leafGridView(),
wells,
std::unordered_map<std::string, std::set<int>>{},
zoltan_ctrl);
// Verify number of partitions
BOOST_CHECK_EQUAL(num_part, 4);
BOOST_CHECK_EQUAL(part.size(), 7);
// The key tests:
// 1. All cells connected by Well 1 should be in the same partition
BOOST_CHECK_EQUAL(part[1], part[2]);
BOOST_CHECK_EQUAL(part[2], part[3]);
// 2. All cells connected by Well 2 should be in the same partition
BOOST_CHECK_EQUAL(part[3], part[4]);
// 3. Due to the overlap at cell 3, all well-connected cells should be in the same partition
BOOST_CHECK_EQUAL(part[1], part[4]);
// 4. Cells 0, 5, and 6 can be in any partition
#endif
}
BOOST_AUTO_TEST_CASE(PartitionCellsWithOverlappingWells3DTest)
{
// Create a 3x3x3 grid with three wells:
// - Well1: Vertical well through (1,1,0) -> (1,1,1) -> (1,1,2)
// - Well2: Diagonal well through (0,0,0) -> (1,1,1) -> (2,2,2) [overlaps with Well1 at (1,1,1)]
// - Well3: Horizontal well through (0,2,0) -> (1,2,0) -> (2,2,0) [no overlap]
auto grid = createTestGrid({3, 3, 3}, {3.0, 3.0, 3.0});
using Entity = typename Dune::CpGrid::LeafGridView::template Codim<0>::Entity;
auto zoltan_ctrl = createZoltanControl<Entity>(grid);
// Helper to convert i,j,k coordinates to global index for 3x3x3 grid
auto ijkToGlobal = [](int i, int j, int k) { return i + (3 * j) + (9 * k); };
// Create wells with proper 3D connections
auto wells = std::vector<Dune::cpgrid::OpmWellType>();
// Well 1 - vertical well
{
auto well_conn = std::make_shared<Opm::WellConnections>();
well_conn->add(createConnection(1, 1, 0, ijkToGlobal(1, 1, 0))); // Global index: 4
well_conn->add(createConnection(1, 1, 1, ijkToGlobal(1, 1, 1))); // Global index: 13
well_conn->add(createConnection(1, 1, 2, ijkToGlobal(1, 1, 2))); // Global index: 22
auto well = createWell("VERTICAL");
well.updateConnections(well_conn, true);
wells.push_back(well);
}
// Well 2 - diagonal well (overlaps with Well 1 at (1,1,1))
{
auto well_conn = std::make_shared<Opm::WellConnections>();
well_conn->add(createConnection(0, 0, 0, ijkToGlobal(0, 0, 0))); // Global index: 0
well_conn->add(createConnection(1, 1, 1, ijkToGlobal(1, 1, 1))); // Global index: 13
well_conn->add(createConnection(2, 2, 2, ijkToGlobal(2, 2, 2))); // Global index: 26
auto well = createWell("DIAGONAL");
well.updateConnections(well_conn, true);
wells.push_back(well);
}
// Well 3 - horizontal well (no overlap)
{
auto well_conn = std::make_shared<Opm::WellConnections>();
well_conn->add(createConnection(0, 2, 0, ijkToGlobal(0, 2, 0))); // Global index: 6
well_conn->add(createConnection(1, 2, 0, ijkToGlobal(1, 2, 0))); // Global index: 7
well_conn->add(createConnection(2, 2, 0, ijkToGlobal(2, 2, 0))); // Global index: 8
auto well = createWell("HORIZONTAL");
well.updateConnections(well_conn, true);
wells.push_back(well);
}
// Create well connections object
Dune::cpgrid::WellConnections wellConnections(wells, std::unordered_map<std::string, std::set<int>>{}, grid);
#if HAVE_MPI && HAVE_ZOLTAN
// Test Zoltan partitioning with overlapping well cells
auto [part, num_part] = Opm::partitionCells("zoltan", 10, grid.leafGridView(),
wells,
std::unordered_map<std::string, std::set<int>>{},
zoltan_ctrl);
// Verify number of partitions
BOOST_CHECK_EQUAL(num_part, 10);
BOOST_CHECK_EQUAL(part.size(), 27); // 3x3x3 grid
// The key tests:
// 1. All cells in Well 1 (vertical) should be in the same partition
BOOST_CHECK_EQUAL(part[ijkToGlobal(1,1,0)], part[ijkToGlobal(1,1,1)]);
BOOST_CHECK_EQUAL(part[ijkToGlobal(1,1,1)], part[ijkToGlobal(1,1,2)]);
// 2. All cells in Well 2 (diagonal) should be in the same partition
BOOST_CHECK_EQUAL(part[ijkToGlobal(0,0,0)], part[ijkToGlobal(1,1,1)]);
BOOST_CHECK_EQUAL(part[ijkToGlobal(1,1,1)], part[ijkToGlobal(2,2,2)]);
// 3. Due to overlap at (1,1,1), all cells from Well 1 and Well 2 should be in same partition
BOOST_CHECK_EQUAL(part[ijkToGlobal(1,1,0)], part[ijkToGlobal(0,0,0)]);
BOOST_CHECK_EQUAL(part[ijkToGlobal(1,1,2)], part[ijkToGlobal(2,2,2)]);
// 4. All cells in Well 3 (horizontal) should be in the same partition
BOOST_CHECK_EQUAL(part[ijkToGlobal(0,2,0)], part[ijkToGlobal(1,2,0)]);
BOOST_CHECK_EQUAL(part[ijkToGlobal(1,2,0)], part[ijkToGlobal(2,2,0)]);
// 5. Well 3 can be in a different partition from Well 1 and 2 since it doesn't overlap
// No test needed as it can be in any partition
#endif
}
BOOST_AUTO_TEST_CASE(PartitionCellsComplexWellNetworkTest)
{
// Create a 5x5x4 grid to accommodate a complex well network
auto grid = createTestGrid({5, 5, 4}, {5.0, 5.0, 4.0});
using Entity = typename Dune::CpGrid::LeafGridView::template Codim<0>::Entity;
auto zoltan_ctrl = createZoltanControl<Entity>(grid);
// Helper to convert i,j,k coordinates to global index for 5x5x4 grid
auto ijkToGlobal = [](int i, int j, int k) { return i + (5 * j) + (25 * k); };
// Create wells with complex connections
auto wells = std::vector<Dune::cpgrid::OpmWellType>();
// Well 1 - Main vertical producer with two branches
{
auto well_conn = std::make_shared<Opm::WellConnections>();
// Main wellbore - vertical
well_conn->add(createConnection(2, 2, 0, ijkToGlobal(2, 2, 0))); // Bottom
well_conn->add(createConnection(2, 2, 1, ijkToGlobal(2, 2, 1)));
well_conn->add(createConnection(2, 2, 2, ijkToGlobal(2, 2, 2)));
well_conn->add(createConnection(2, 2, 3, ijkToGlobal(2, 2, 3))); // Top
// Branch 1 - horizontal in x direction
well_conn->add(createConnection(3, 2, 2, ijkToGlobal(3, 2, 2)));
well_conn->add(createConnection(4, 2, 2, ijkToGlobal(4, 2, 2)));
// Branch 2 - horizontal in y direction
well_conn->add(createConnection(2, 3, 1, ijkToGlobal(2, 3, 1)));
well_conn->add(createConnection(2, 4, 1, ijkToGlobal(2, 4, 1)));
auto well = createWell("PRODUCER1");
well.updateConnections(well_conn, true);
wells.push_back(well);
}
// Well 2 - Diagonal injector crossing the producer
{
auto well_conn = std::make_shared<Opm::WellConnections>();
well_conn->add(createConnection(0, 0, 0, ijkToGlobal(0, 0, 0)));
well_conn->add(createConnection(1, 1, 1, ijkToGlobal(1, 1, 1)));
well_conn->add(createConnection(2, 2, 2, ijkToGlobal(2, 2, 2))); // Intersects with Well 1
well_conn->add(createConnection(3, 3, 3, ijkToGlobal(3, 3, 3)));
auto well = createWell("INJECTOR1");
well.updateConnections(well_conn, true);
wells.push_back(well);
}
// Well 3 - Horizontal well that crosses both previous wells
{
auto well_conn = std::make_shared<Opm::WellConnections>();
well_conn->add(createConnection(0, 2, 2, ijkToGlobal(0, 2, 2)));
well_conn->add(createConnection(1, 2, 2, ijkToGlobal(1, 2, 2)));
well_conn->add(createConnection(2, 2, 2, ijkToGlobal(2, 2, 2))); // Intersects with Well 1 and 2
well_conn->add(createConnection(3, 2, 2, ijkToGlobal(3, 2, 2))); // Intersects with Well 1's branch
well_conn->add(createConnection(4, 2, 2, ijkToGlobal(4, 2, 2)));
auto well = createWell("PRODUCER2");
well.updateConnections(well_conn, true);
wells.push_back(well);
}
// Well 4 - L-shaped well connecting to the network
{
auto well_conn = std::make_shared<Opm::WellConnections>();
well_conn->add(createConnection(2, 4, 1, ijkToGlobal(2, 4, 1))); // Connects to Well 1's branch
well_conn->add(createConnection(2, 4, 2, ijkToGlobal(2, 4, 2)));
well_conn->add(createConnection(2, 4, 3, ijkToGlobal(2, 4, 3)));
well_conn->add(createConnection(3, 4, 3, ijkToGlobal(3, 4, 3)));
well_conn->add(createConnection(4, 4, 3, ijkToGlobal(4, 4, 3)));
auto well = createWell("PRODUCER3");
well.updateConnections(well_conn, true);
wells.push_back(well);
}
// Well 5 - Isolated well (not connected to the network)
{
auto well_conn = std::make_shared<Opm::WellConnections>();
well_conn->add(createConnection(0, 4, 0, ijkToGlobal(0, 4, 0)));
well_conn->add(createConnection(1, 4, 0, ijkToGlobal(1, 4, 0)));
auto well = createWell("ISOLATED");
well.updateConnections(well_conn, true);
wells.push_back(well);
}
#if HAVE_MPI && HAVE_ZOLTAN
// Test Zoltan partitioning with complex well network
auto [part, num_part] = Opm::partitionCells("zoltan", 15, grid.leafGridView(),
wells,
std::unordered_map<std::string, std::set<int>>{},
zoltan_ctrl);
// Verify number of partitions and grid size
BOOST_CHECK_EQUAL(num_part, 15);
BOOST_CHECK_EQUAL(part.size(), 100); // 5x5x4 grid
// Helper to check if two cells are in the same partition
auto inSamePartition = [&part](int idx1, int idx2) {
return part[idx1] == part[idx2];
};
// Test 1: All cells in Well 1 (main bore and branches) should be in the same partition
// Main wellbore
BOOST_CHECK(inSamePartition(ijkToGlobal(2,2,0), ijkToGlobal(2,2,1)));
BOOST_CHECK(inSamePartition(ijkToGlobal(2,2,1), ijkToGlobal(2,2,2)));
BOOST_CHECK(inSamePartition(ijkToGlobal(2,2,2), ijkToGlobal(2,2,3)));
// Branch 1
BOOST_CHECK(inSamePartition(ijkToGlobal(2,2,2), ijkToGlobal(3,2,2)));
BOOST_CHECK(inSamePartition(ijkToGlobal(3,2,2), ijkToGlobal(4,2,2)));
// Branch 2
BOOST_CHECK(inSamePartition(ijkToGlobal(2,2,1), ijkToGlobal(2,3,1)));
BOOST_CHECK(inSamePartition(ijkToGlobal(2,3,1), ijkToGlobal(2,4,1)));
// Test 2: Well 2 (diagonal) should be in same partition as Well 1 due to intersection
BOOST_CHECK(inSamePartition(ijkToGlobal(0,0,0), ijkToGlobal(2,2,2)));
BOOST_CHECK(inSamePartition(ijkToGlobal(3,3,3), ijkToGlobal(2,2,2)));
// Test 3: Well 3 (horizontal) should be in same partition as Well 1 and 2
BOOST_CHECK(inSamePartition(ijkToGlobal(0,2,2), ijkToGlobal(2,2,2)));
BOOST_CHECK(inSamePartition(ijkToGlobal(4,2,2), ijkToGlobal(2,2,2)));
// Test 4: Well 4 (L-shaped) should be in same partition as Well 1 due to connection
BOOST_CHECK(inSamePartition(ijkToGlobal(2,4,1), ijkToGlobal(2,2,1)));
BOOST_CHECK(inSamePartition(ijkToGlobal(4,4,3), ijkToGlobal(2,2,1)));
// Test 5: Well 5 (isolated) should maintain its own cells in same partition
BOOST_CHECK(inSamePartition(ijkToGlobal(0,4,0), ijkToGlobal(1,4,0)));
// But should be in different partition from the connected network
BOOST_CHECK(part[ijkToGlobal(0,4,0)] != part[ijkToGlobal(2,2,2)]);
// Test 6: Verify the connected network (Wells 1-4) forms a single partition
std::set<int> network_partition_ids;
network_partition_ids.insert(part[ijkToGlobal(2,2,2)]); // Well 1
network_partition_ids.insert(part[ijkToGlobal(1,1,1)]); // Well 2
network_partition_ids.insert(part[ijkToGlobal(2,2,2)]); // Well 3
network_partition_ids.insert(part[ijkToGlobal(2,4,1)]); // Well 4
BOOST_CHECK_EQUAL(network_partition_ids.size(), 1);
#endif
}
BOOST_AUTO_TEST_CASE(PartitionCellsWithNonReachableCellsTest)
{
// Create a 3x3x3 grid where:
// - cell (0,0,0) is isolated (surrounded by inactive cells)
// - cell (0,2,0) is isolated (surrounded by inactive cells)
// - cells needed for wells are active
// Grid structure (top view):
// Layer k=0: Layer k=1: Layer k=2: W = well cells
// [I][X][X ] [X][ ][ ] [X][ ][ ] X = inactive cells
// [X][X][W2] [X][W1][W2] [X][W1][W2] I = isolated active cell
// [I][X][ ] [X][ ][ ] [X][ ][ ] [ ] = active cells
const std::string deckString =
R"(RUNSPEC
DIMENS
3 3 3 /
GRID
COORD
0 0 0 0 0 1
1 0 0 1 0 1
2 0 0 2 0 1
3 0 0 3 0 1
0 1 0 0 1 1
1 1 0 1 1 1
2 1 0 2 1 1
3 1 0 3 1 1
0 2 0 0 2 1
1 2 0 1 2 1
2 2 0 2 2 1
3 2 0 3 2 1
0 3 0 0 3 1
1 3 0 1 3 1
2 3 0 2 3 1
3 3 0 3 3 1
/
ZCORN
36*0
36*1
36*1
36*2
36*2
36*3
/
ACTNUM
-- First layer (3x3)
1 0 0 -- cell (0,0,0) is isolated
0 0 1
1 0 1 -- cell (1,1,1) is isolated
-- Middle layer (3x3)
0 1 1
0 1 1
0 1 1
-- Top layer (3x3)
0 1 1
0 1 1
0 1 1
/
END
)";
// Parse and create grid
Opm::Parser parser;
const auto deck = parser.parseString(deckString);
Opm::EclipseGrid ecl_grid(deck);
Dune::CpGrid grid;
grid.processEclipseFormat(&ecl_grid, nullptr, false, false, false);
// Verify grid properties
const auto& grid_view = grid.leafGridView();
const std::size_t num_active_cells = grid_view.size(0);
BOOST_CHECK_EQUAL(num_active_cells, 16); // 27 total cells - 11 inactive = 16 active
// Helper functions
using Entity = typename Dune::CpGrid::LeafGridView::template Codim<0>::Entity;
auto zoltan_ctrl = createZoltanControl<Entity>(grid);
auto ijkToGlobal = [](int i, int j, int k) { return i + (3 * j) + (9 * k); };
auto wells = std::vector<Dune::cpgrid::OpmWellType>();
// Well 1 - Middle column well
{
auto well_conn = std::make_shared<Opm::WellConnections>();
well_conn->add(createConnection(1, 1, 1, ijkToGlobal(1, 1, 1)));
well_conn->add(createConnection(1, 1, 2, ijkToGlobal(1, 1, 2)));
auto well = createWell("WELL1");
well.updateConnections(well_conn, true);
wells.push_back(well);
}
// Well 2 - Right column well
{
auto well_conn = std::make_shared<Opm::WellConnections>();
well_conn->add(createConnection(2, 1, 0, ijkToGlobal(2, 1, 0)));
well_conn->add(createConnection(2, 1, 1, ijkToGlobal(2, 1, 1)));
well_conn->add(createConnection(2, 1, 2, ijkToGlobal(2, 1, 2)));
auto well = createWell("WELL2");
well.updateConnections(well_conn, true);
wells.push_back(well);
}
#if HAVE_MPI && HAVE_ZOLTAN
auto [part, num_part] = Opm::partitionCells("zoltan", 5, grid.leafGridView(),
wells,
std::unordered_map<std::string, std::set<int>>{},
zoltan_ctrl);
BOOST_CHECK_EQUAL(num_part, 5);
BOOST_CHECK_EQUAL(part.size(), 16);
// For this test global indices != local indices, so we need to
// create a map from global index to local index
std::map<int, int> g2l;
for (const auto& element : elements(grid.leafGridView(), Dune::Partitions::interior)) {
const auto globalIndex = grid.globalCell()[zoltan_ctrl.index(element)];
g2l[globalIndex] = zoltan_ctrl.index(element);
}
// Verify that the two isolated cells are assigned to the -1 partition
BOOST_CHECK_EQUAL(part[g2l[ijkToGlobal(0,0,0)]], -1);
BOOST_CHECK_EQUAL(part[g2l[ijkToGlobal(0,2,0)]], -1);
// Well 1: All cells should be in the same partition
BOOST_CHECK(part[g2l[ijkToGlobal(1,1,1)]] == part[g2l[ijkToGlobal(1,1,2)]]);
// Well 2: All cells should be in the same partition
BOOST_CHECK(part[g2l[ijkToGlobal(2,1,0)]] == part[g2l[ijkToGlobal(2,1,1)]]);
BOOST_CHECK(part[g2l[ijkToGlobal(2,1,1)]] == part[g2l[ijkToGlobal(2,1,2)]]);
// Wells can be in different partitions as they don't connect
BOOST_CHECK(part[g2l[ijkToGlobal(1,1,1)]] != part[g2l[ijkToGlobal(2,1,2)]]);
#endif
}
BOOST_AUTO_TEST_CASE(PartitionCellsWithNeighborConnectivityTest)
{
// Create a 1D grid with 7 cells and two wells that are not directly connected
// but will be connected through their neighboring cells when using neighbor_levels=1
// Visual representation:
// Cell indices: 0 1 2 3 4 5 6
// Well layout: W1 W2
// Connected: * * * * *
// The wells will be connected because cell 2 is a neighbor of W1 and W2
auto grid = createTestGrid({7, 1, 1}, {7.0, 1.0, 1.0});
using Entity = typename Dune::CpGrid::LeafGridView::template Codim<0>::Entity;
auto zoltan_ctrl = createZoltanControl<Entity>(grid);
// Create two wells that are not directly connected
auto wells = createWellsWithConnections({
{"WELL1", {1}}, // Well 1 in cell 1
{"WELL2", {3}} // Well 2 in cell 3 (not adjacent to Well 1)
});
#if HAVE_MPI && HAVE_ZOLTAN
// Test Zoltan partitioning with neighbor connectivity
auto [part, num_part] = Opm::partitionCells("zoltan", 3, grid.leafGridView(),
wells,
std::unordered_map<std::string, std::set<int>>{},
zoltan_ctrl,
1);
// Verify number of partitions
BOOST_CHECK_EQUAL(num_part, 3);
BOOST_CHECK_EQUAL(part.size(), 7); // 7x1x1 grid
// Check that the wells are connected through their neighbors
BOOST_CHECK(part[0] == part[1]);
BOOST_CHECK(part[1] == part[2]);
BOOST_CHECK(part[2] == part[3]);
BOOST_CHECK(part[3] == part[4]);
// Check that the remaining cells are in different partitions
BOOST_CHECK(part[0] != part[5]);
BOOST_CHECK(part[5] != part[6]);
#endif
}
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);
}