Add Support for Writing NLDD-Like Partitions in Parallel

This commit adds a new (hidden) debugging option,

    DebugEmitCellPartition (--debug-emit-cell-partition)

which, when set, will cause each rank to write a three-column text
file of the form

    MPI_Rank  Cartesian_Index  NLDD_Domain_ID

into the directory

    partition/CaseName

of the run's output directory.  That file will be named according to
the process' MPI rank, so the first column will be the same as the
file name.

The option is primarily intended for debugging the NLDD partitioning
scheme, so is mostly reserved for runs with low MPI sizes (e.g.,
less than 20).

While here, also make the MPIPartitionFromFile helper class aware of
this format so that we can use concatenated output files as an input
to the MPI partitioning algorithm for repeatability.

ebos/eclcpgridvanguard.hh

@@ -45,6 +45,7 @@
 #include <memory>
 #include <stdexcept>
 #include <tuple>
+#include <unordered_map>
 #include <vector>
 
 #include <fmt/format.h>
@@ -75,16 +76,43 @@ namespace Opm { namespace details {
             std::istream_iterator<int> {}
         };
 
-        if (partition.size() != static_cast<std::vector<int>::size_type>(grid.size(0))) {
-            throw std::invalid_argument {
-                fmt::format("Partition file '{}' with {} values does "
-                            "not match CpGrid instance with {} cells",
-                            this->partitionFile_.generic_string(),
-                            partition.size(), grid.size(0))
-            };
+        const auto nc =
+            static_cast<std::vector<int>::size_type>(grid.size(0));
+
+        if (partition.size() == nc) {
+            // Input is one process ID for each active cell
+            return partition;
         }
 
-        return partition;
+        if (partition.size() == 3 * nc) {
+            // Partition file is of the form
+            //
+            //   Process_ID   Cartesian_Idx   NLDD_Domain
+            //
+            // with one row for each active cell.  Select first column.
+            auto g2l = std::unordered_map<int, int>{};
+            auto locCell = 0;
+            for (const auto& globCell : grid.globalCell()) {
+                g2l.insert_or_assign(globCell, locCell++);
+            }
+
+            auto filtered_partition = std::vector<int>(nc);
+            for (auto c = 0*nc; c < nc; ++c) {
+                auto pos = g2l.find(partition[3*c + 1]);
+                if (pos != g2l.end()) {
+                    filtered_partition[pos->second] = partition[3*c + 0];
+                }
+            }
+
+            return filtered_partition;
+        }
+
+        throw std::invalid_argument {
+            fmt::format("Partition file '{}' with {} values does "
+                        "not match CpGrid instance with {} cells",
+                        this->partitionFile_.generic_string(),
+                        partition.size(), nc)
+        };
     }
 }} // namespace Opm::details
 

opm/simulators/flow/BlackoilModelEbos.hpp

@@ -59,6 +59,8 @@
 #include <algorithm>
 #include <cassert>
 #include <cmath>
+#include <filesystem>
+#include <fstream>
 #include <iomanip>
 #include <ios>
 #include <limits>
@@ -1031,6 +1033,29 @@ namespace Opm {
             return convergence_reports_;
         }
 
+        void writePartitions(const std::filesystem::path& odir) const
+        {
+            if (this->nlddSolver_ != nullptr) {
+                this->nlddSolver_->writePartitions(odir);
+                return;
+            }
+
+            const auto& elementMapper = this->ebosSimulator().model().elementMapper();
+            const auto& cartMapper = this->ebosSimulator().vanguard().cartesianIndexMapper();
+
+            const auto& grid = this->ebosSimulator().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()) };
+
+            for (const auto& cell : elements(grid.leafGridView(), Dune::Partitions::interior)) {
+                pfile << comm.rank() << ' '
+                      << cartMapper.cartesianIndex(elementMapper.index(cell)) << ' '
+                      << comm.rank() << '\n';
+            }
+        }
+
     protected:
         // ---------  Data members  ---------
 

opm/simulators/flow/BlackoilModelEbosNldd.hpp

@@ -54,6 +54,8 @@
 #include <cassert>
 #include <cmath>
 #include <cstddef>
+#include <filesystem>
+#include <fstream>
 #include <functional>
 #include <iomanip>
 #include <ios>
@@ -303,6 +305,26 @@ public:
         return local_reports_accumulated_;
     }
 
+    void writePartitions(const std::filesystem::path& odir) const
+    {
+        const auto& elementMapper = this->model_.ebosSimulator().model().elementMapper();
+        const auto& cartMapper = this->model_.ebosSimulator().vanguard().cartesianIndexMapper();
+
+        const auto& grid = this->model_.ebosSimulator().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()) };
+
+        const auto p = this->reconstitutePartitionVector();
+        auto i = 0;
+        for (const auto& cell : elements(grid.leafGridView(), Dune::Partitions::interior)) {
+            pfile << comm.rank() << ' '
+                  << cartMapper.cartesianIndex(elementMapper.index(cell)) << ' '
+                  << p[i++] << '\n';
+        }
+    }
+
 private:
     //! \brief Solve the equation system for a single domain.
     std::pair<SimulatorReportSingle, ConvergenceReport>
@@ -866,6 +888,34 @@ private:
                                      grid.leafGridView(), wells, zoltan_ctrl);
     }
 
+    std::vector<int> reconstitutePartitionVector() const
+    {
+        const auto& grid = this->model_.ebosSimulator().vanguard().grid();
+
+        auto numD = std::vector<int>(grid.comm().size() + 1, 0);
+        numD[grid.comm().rank() + 1] = static_cast<int>(this->domains_.size());
+        grid.comm().sum(numD.data(), numD.size());
+        std::partial_sum(numD.begin(), numD.end(), numD.begin());
+
+        auto p = std::vector<int>(grid.size(0));
+        auto maxCellIdx = std::numeric_limits<int>::min();
+
+        auto d = numD[grid.comm().rank()];
+        for (const auto& domain : this->domains_) {
+            for (const auto& cell : domain.cells) {
+                p[cell] = d;
+                if (cell > maxCellIdx) {
+                    maxCellIdx = cell;
+                }
+            }
+
+            ++d;
+        }
+
+        p.erase(p.begin() + maxCellIdx + 1, p.end());
+        return p;
+    }
+
     BlackoilModelEbos<TypeTag>& model_; //!< Reference to model
     std::vector<Domain> domains_; //!< Vector of subdomains
     std::vector<std::unique_ptr<Mat>> domain_matrices_; //!< Vector of matrix operator for each subdomain

opm/simulators/flow/BlackoilModelParametersEbos.hpp

@@ -117,6 +117,10 @@ template<class TypeTag, class MyTypeTag>
 struct EnableWellOperabilityCheckIter {
     using type = UndefinedProperty;
 };
+template<class TypeTag, class MyTypeTag>
+struct DebugEmitCellPartition {
+    using type = UndefinedProperty;
+};
 // parameters for multisegment wells
 template<class TypeTag, class MyTypeTag>
 struct TolerancePressureMsWells {
@@ -357,6 +361,10 @@ struct EnableWellOperabilityCheckIter<TypeTag, TTag::FlowModelParameters> {
     static constexpr bool value = false;
 };
 template<class TypeTag>
+struct DebugEmitCellPartition<TypeTag, TTag::FlowModelParameters> {
+    static constexpr bool value = false;
+};
+template<class TypeTag>
 struct RelaxedWellFlowTol<TypeTag, TTag::FlowModelParameters> {
     using type = GetPropType<TypeTag, Scalar>;
     static constexpr type value = 1e-3;
@@ -573,6 +581,8 @@ namespace Opm
         std::string local_domain_partition_method_;
         DomainOrderingMeasure local_domain_ordering_{DomainOrderingMeasure::AveragePressure};
 
+        bool write_partitions_{false};
+
         /// Construct from user parameters or defaults.
         BlackoilModelParametersEbos()
         {
@@ -637,6 +647,8 @@ namespace Opm
             } else {
                 throw std::runtime_error("Invalid domain ordering '" + measure + "' specified.");
             }
+
+            write_partitions_ = EWOMS_GET_PARAM(TypeTag, bool, DebugEmitCellPartition);
         }
 
         static void registerParameters()
@@ -690,6 +702,10 @@ namespace Opm
                                  "Allowed values are 'zoltan', 'simple', and the name of a partition file ending with '.partition'.");
             EWOMS_REGISTER_PARAM(TypeTag, std::string, LocalDomainsOrderingMeasure, "Subdomain ordering measure. "
                                  "Allowed values are 'pressure' and  'residual'.");
+
+            EWOMS_REGISTER_PARAM(TypeTag, bool, DebugEmitCellPartition, "Whether or not to emit cell partitions as a debugging aid.");
+
+            EWOMS_HIDE_PARAM(TypeTag, DebugEmitCellPartition);
         }
     };
 } // namespace Opm

opm/simulators/flow/SimulatorFullyImplicitBlackoilEbos.hpp

@@ -44,6 +44,7 @@
 
 #include <boost/date_time/gregorian/gregorian.hpp>
 
+#include <filesystem>
 #include <memory>
 #include <optional>
 #include <string>
@@ -537,6 +538,24 @@ protected:
                                              wellModel,
                                              terminalOutput_);
 
+        if (this->modelParam_.write_partitions_) {
+            const auto& iocfg = this->eclState().cfg().io();
+
+            const auto odir = iocfg.getOutputDir()
+                / std::filesystem::path { "partition" }
+                / iocfg.getBaseName();
+
+            if (this->grid().comm().rank() == 0) {
+                create_directories(odir);
+            }
+
+            this->grid().comm().barrier();
+
+            model->writePartitions(odir);
+
+            this->modelParam_.write_partitions_ = false;
+        }
+
         return std::make_unique<Solver>(solverParam_, std::move(model));
     }