Add Multi-Array Aware Accumulator for Inter-Region Flow Rates

This commit introduces a new helper class

    Opm::EclInterRegFlowMapSingleFIP

that wraps a vector of

    Opm::EclInterRegFlowMapSingleFIP

along with the associate array names (e.g., "FIPNUM" and any
additional "FIP*" arrays).  We implement the same operations as the
*SingleFIP type and simply loop over the internal accumulators to
affect the operations.

Add unit tests to exercise the new class, including simulating
multiple MPI ranks that are communicated to a single I/O rank for
summary output purposes.

ebos/eclinterregflows.cc

@@ -50,9 +50,11 @@ addConnection(const Cell& source,
         };
     }
 
-    if (! (source.isInterior || destination.isInterior)) {
-        // Connection between two cells not on this process.  Unlikely, but
-        // nothing to do here.
+    if (! source.isInterior ||
+        (source.cartesianIndex > destination.cartesianIndex))
+    {
+        // Connection handled in different call.  Don't double-count
+        // contributions.
         return;
     }
 
@@ -64,14 +66,9 @@ addConnection(const Cell& source,
         return;
     }
 
-    if ((source.isInterior && destination.isInterior) ||
-        (source.isInterior && (r1 < r2)) ||
-        (destination.isInterior && (r2 < r1)))
-    {
-        // Inter-region connection internal to an MPI rank or this rank owns
-        // the flow rate across this connection.
-        this->iregFlow_.addConnection(r1, r2, rates);
-    }
+    // Inter-region connection internal to an MPI rank or this rank owns
+    // the flow rate across this connection.
+    this->iregFlow_.addConnection(r1, r2, rates);
 }
 
 void Opm::EclInterRegFlowMapSingleFIP::compress()
@@ -107,3 +104,109 @@ assignGlobalMaxRegionID(const std::size_t regID)
     this->maxGlobalRegionID_ = regID;
     return true;
 }
+
+// =====================================================================
+//
+// Implementation of EclInterRegFlowMap (wrapper for multiple arrays)
+//
+// =====================================================================
+
+Opm::EclInterRegFlowMap::
+EclInterRegFlowMap(const std::size_t                numCells,
+                   const std::vector<SingleRegion>& regions)
+{
+    this->regionMaps_.reserve(regions.size());
+    this->names_.reserve(regions.size());
+
+    this->numCells_ = numCells;
+
+    for (const auto& region : regions) {
+        this->regionMaps_.emplace_back(region.definition);
+        this->names_.push_back(region.name);
+    }
+}
+
+void
+Opm::EclInterRegFlowMap::
+addConnection(const Cell& source,
+              const Cell& destination,
+              const data::InterRegFlowMap::FlowRates& rates)
+{
+    for (auto& regionMap : this->regionMaps_) {
+        regionMap.addConnection(source, destination, rates);
+    }
+}
+
+void Opm::EclInterRegFlowMap::compress()
+{
+    for (auto& regionMap : this->regionMaps_) {
+        regionMap.compress();
+    }
+}
+
+void Opm::EclInterRegFlowMap::clear()
+{
+    for (auto& regionMap : this->regionMaps_) {
+        regionMap.clear();
+    }
+
+    this->readIsConsistent_ = true;
+}
+
+const std::vector<std::string>&
+Opm::EclInterRegFlowMap::names() const
+{
+    return this->names_;
+}
+
+std::vector<Opm::data::InterRegFlowMap>
+Opm::EclInterRegFlowMap::getInterRegFlows() const
+{
+    auto maps = std::vector<data::InterRegFlowMap>{};
+    maps.reserve(this->regionMaps_.size());
+
+    for (auto& regionMap : this->regionMaps_) {
+        maps.push_back(regionMap.getInterRegFlows());
+    }
+
+    return maps;
+}
+
+std::vector<std::size_t>
+Opm::EclInterRegFlowMap::getLocalMaxRegionID() const
+{
+    auto maxLocalRegionID = std::vector<std::size_t>{};
+    maxLocalRegionID.reserve(this->regionMaps_.size());
+
+    for (auto& regionMap : this->regionMaps_) {
+        maxLocalRegionID.push_back(regionMap.getLocalMaxRegionID());
+    }
+
+    return maxLocalRegionID;
+}
+
+bool
+Opm::EclInterRegFlowMap::
+assignGlobalMaxRegionID(const std::vector<std::size_t>& regID)
+{
+    if (regID.size() != this->regionMaps_.size()) {
+        return false;
+    }
+
+    auto assignmentOK = true;
+
+    const auto numReg = regID.size();
+    for (auto region = 0*numReg; region < numReg; ++region) {
+        const auto ok = this->regionMaps_[region]
+            .assignGlobalMaxRegionID(regID[region]);
+
+        assignmentOK = assignmentOK && ok;
+    }
+
+    return assignmentOK;
+}
+
+bool Opm::EclInterRegFlowMap::readIsConsistent() const
+{
+    return this->readIsConsistent_;
+}

ebos/eclinterregflows.hh

@@ -49,6 +49,9 @@ namespace Opm {
             /// Cell's active index on local rank.
             int activeIndex{-1};
 
+            /// Cell's global cell ID.
+            int cartesianIndex{-1};
+
             /// Whether or not cell is interior to local rank.
             bool isInterior{true};
         };
@@ -163,6 +166,197 @@ namespace Opm {
         /// from a stream.  For error detection.
         bool isReadFromStream_{false};
     };
+
+    /// Inter-region flow accumulation maps for all region definition arrays
+    class EclInterRegFlowMap
+    {
+    public:
+        /// Minimal representation of a single named region defintion.
+        ///
+        /// Name is typically "FIPNUM" or any additional "FIP*" array names.
+        struct SingleRegion {
+            /// Region definition array name
+            std::string name;
+
+            /// Region definition array.
+            std::reference_wrapper<const std::vector<int>> definition;
+        };
+
+        /// Characteristics of a cell from a simulation grid.
+        using Cell = EclInterRegFlowMapSingleFIP::Cell;
+
+        /// Default constructor.
+        EclInterRegFlowMap() = default;
+
+        /// Constructor.
+        ///
+        /// \param[in] numCells Number of cells on local MPI rank, including
+        ///    overlap cells if applicable.
+        ///
+        /// \param[in] regions All applicable region definition arrays.
+        explicit EclInterRegFlowMap(const std::size_t                numCells,
+                                    const std::vector<SingleRegion>& regions);
+
+        EclInterRegFlowMap(const EclInterRegFlowMap& rhs) = default;
+        EclInterRegFlowMap(EclInterRegFlowMap&& rhs) noexcept = default;
+
+        EclInterRegFlowMap& operator=(const EclInterRegFlowMap& rhs) = default;
+        EclInterRegFlowMap& operator=(EclInterRegFlowMap&& rhs) noexcept = default;
+
+        /// Add flow rate connection between regions for all region
+        /// definitions.
+        ///
+        /// \param[in] source Cell from which the flow nominally originates.
+        ///
+        /// \param[in] destination Cell into which flow nominally goes.
+        ///
+        /// \param[in] rates Flow rates associated to single connection.
+        ///
+        /// If both cells are in the same region, or if neither cell is
+        /// interior to this MPI rank, then this function does nothing.  If
+        /// one cell is interior to this MPI rank and the other isn't, then
+        /// this function will include the flow rate contribution if and
+        /// only if the cell with the smallest associate region ID is
+        /// interior to this MPI rank.
+        void addConnection(const Cell& source,
+                           const Cell& destination,
+                           const data::InterRegFlowMap::FlowRates& rates);
+
+        /// Form CSR adjacency matrix representation of input graph from
+        /// connections established in previous calls to addConnection().
+        ///
+        /// Number of rows in the CSR representation is the maximum FIP
+        /// region ID.
+        void compress();
+
+        /// Clear all internal buffers, but preserve allocated capacity.
+        void clear();
+
+        /// Names of all applicable region definition arrays.
+        ///
+        /// Mostly intended for summary output purposes.
+        const std::vector<std::string>& names() const;
+
+        /// Get read-only access to the underlying CSR representation.
+        ///
+        /// Mostly intended for summary output purposes.
+        std::vector<data::InterRegFlowMap> getInterRegFlows() const;
+
+        /// Retrieve maximum FIP region ID on local MPI rank.
+        std::vector<std::size_t> getLocalMaxRegionID() const;
+
+        /// Assign maximum FIP region ID across all MPI ranks.
+        ///
+        /// Fails if global maximum is smaller than local maximum region ID.
+        ///
+        /// \param[in] regID Global maximum FIP region ID for this FIP
+        ///    definition array across all MPI ranks.
+        ///
+        /// \return Whether or not assignment succeeded.
+        bool assignGlobalMaxRegionID(const std::vector<std::size_t>& regID);
+
+        /// Whether or not previous read() operation succeeded.
+        bool readIsConsistent() const;
+
+        /// Serialise internal representation to MPI message buffer
+        ///
+        /// \tparam MessageBufferType Linear MPI message buffer.  API should
+        ///    be similar to Dune::MessageBufferIF
+        ///
+        /// \param[in,out] buffer Linear MPI message buffer instance.
+        ///    Function appends a partially linearised representation of
+        ///    \code *this \endcode to the buffer contents.
+        template <class MessageBufferType>
+        void write(MessageBufferType& buffer) const
+        {
+            buffer.write(this->names_.size());
+            for (const auto& name : this->names_) {
+                buffer.write(name);
+            }
+
+            for (const auto& map : this->regionMaps_) {
+                map.write(buffer);
+            }
+        }
+
+        /// Reconstitute internal object representation from MPI message
+        /// buffer
+        ///
+        /// This object (\code *this \endcode) is not usable in subsequent
+        /// calls to \code addConnection() \endcode following a call to
+        /// member function \code read() \endcode.
+        ///
+        /// \tparam MessageBufferType Linear MPI message buffer.  API should
+        ///    be similar to Dune::MessageBufferIF
+        ///
+        /// \param[in,out] buffer Linear MPI message buffer instance.
+        ///    Function reads a partially linearised representation of \code
+        ///    *this \endcode from the buffer contents and advances the
+        ///    buffer's read position.
+        template <class MessageBufferType>
+        void read(MessageBufferType& buffer)
+        {
+            const auto names = this->readNames(buffer);
+
+            if (names == this->names_) {
+                // Input stream holds the maps in expected order (common
+                // case).  Read the maps and merge with internal values.
+                for (auto& map : this->regionMaps_) {
+                    map.read(buffer);
+                }
+            }
+            else {
+                // Input stream does not hold the maps in expected order (or
+                // different number of maps).  Unexpected.  Read the values
+                // from the input stream, but do not merge with internal
+                // values.
+                auto map = EclInterRegFlowMapSingleFIP {
+                    std::vector<int>(this->numCells_, 1)
+                };
+
+                const auto numMaps = names.size();
+                for (auto n = 0*numMaps; n < numMaps; ++n) {
+                    map.read(buffer);
+                }
+
+                this->readIsConsistent_ = false;
+            }
+        }
+
+    private:
+        /// Inter-region flow accumulators.  One accumulator map for each
+        /// region definition array.
+        std::vector<EclInterRegFlowMapSingleFIP> regionMaps_{};
+
+        /// Names of region definition arrays.  Typically "FIPNUM" and other
+        /// "FIPXYZ" array names.
+        std::vector<std::string> names_;
+
+        /// Number of cells, including overlap, reachable from this MPI
+        /// rank.
+        std::size_t numCells_{0};
+
+        /// Wheter or not read() successfully updated this object from an
+        /// input stream.
+        bool readIsConsistent_{true};
+
+        /// Retrieve array names from an input stream.
+        ///
+        /// Needed for consistency checks.
+        template <class MessageBufferType>
+        std::vector<std::string> readNames(MessageBufferType& buffer) const
+        {
+            auto numNames = 0 * this->names_.size();
+            buffer.read(numNames);
+
+            auto names = std::vector<std::string>(numNames);
+            for (auto name = 0*numNames; name < numNames; ++name) {
+                buffer.read(names[name]);
+            }
+
+            return names;
+        }
+    };
 } // namespace Opm
 
 #endif // ECL_INTERREG_FLOWS_MODULE_HH