Merge pull request #3796 from bska/accumulate-iregflow-values

Accumulate Inter-Region Phase Flow Rate Contributions

CMakeLists_files.cmake

@@ -31,6 +31,7 @@ list (APPEND MAIN_SOURCE_FILES
   ebos/eclgenerictracermodel.cc
   ebos/eclgenericvanguard.cc
   ebos/eclgenericwriter.cc
+  ebos/eclinterregflows.cc
   ebos/ecltransmissibility.cc
   opm/core/props/BlackoilPhases.cpp
   opm/core/props/phaseUsageFromDeck.cpp
@@ -135,30 +136,31 @@ endif()
 # originally generated with the command:
 # find tests -name '*.cpp' -a ! -wholename '*/not-unit/*' -printf '\t%p\n' | sort
 list (APPEND TEST_SOURCE_FILES
-  tests/test_equil.cc
-  tests/test_ecl_output.cc
+  tests/test_ALQState.cpp
   tests/test_blackoil_amg.cpp
   tests/test_convergencereport.cpp
+  tests/test_deferredlogger.cpp
+  tests/test_ecl_output.cc
+  tests/test_eclinterregflows.cpp
+  tests/test_equil.cc
   tests/test_flexiblesolver.cpp
-  tests/test_preconditionerfactory.cpp
+  tests/test_glift1.cpp
   tests/test_graphcoloring.cpp
-  tests/test_vfpproperties.cpp
+  tests/test_GroupState.cpp
+  tests/test_invert.cpp
+  tests/test_keyword_validator.cpp
   tests/test_milu.cpp
   tests/test_multmatrixtransposed.cpp
-  tests/test_wellmodel.cpp
-  tests/test_deferredlogger.cpp
-  tests/test_timer.cpp
-  tests/test_invert.cpp
-  tests/test_stoppedwells.cpp
-  tests/test_relpermdiagnostics.cpp
   tests/test_norne_pvt.cpp
+  tests/test_parallelwellinfo.cpp
+  tests/test_preconditionerfactory.cpp
+  tests/test_relpermdiagnostics.cpp
+  tests/test_stoppedwells.cpp
+  tests/test_timer.cpp
+  tests/test_vfpproperties.cpp
+  tests/test_wellmodel.cpp
   tests/test_wellprodindexcalculator.cpp
   tests/test_wellstate.cpp
-  tests/test_parallelwellinfo.cpp
-  tests/test_glift1.cpp
-  tests/test_keyword_validator.cpp
-  tests/test_GroupState.cpp
-  tests/test_ALQState.cpp
   )
 
 if(MPI_FOUND)
@@ -229,6 +231,7 @@ list (APPEND TEST_DATA_FILES
 # originally generated with the command:
 # find opm -name '*.h*' -a ! -name '*-pch.hpp' -printf '\t%p\n' | sort
 list (APPEND PUBLIC_HEADER_FILES
+  ebos/eclinterregflows.hh
   opm/simulators/flow/countGlobalCells.hpp
   opm/simulators/flow/BlackoilModelEbos.hpp
   opm/simulators/flow/BlackoilModelParametersEbos.hpp

ebos/eclgenericoutputblackoilmodule.cc

@@ -18,6 +18,7 @@
 */
 
 #include <config.h>
+
 #include <ebos/eclgenericoutputblackoilmodule.hh>
 
 #include <opm/common/OpmLog/OpmLog.hpp>
@@ -33,11 +34,16 @@
 #include <opm/input/eclipse/Units/Units.hpp>
 
 #include <cassert>
-#include <fmt/format.h>
+#include <cstddef>
+#include <functional>
 #include <initializer_list>
 #include <iomanip>
 #include <sstream>
 #include <stdexcept>
+#include <string>
+#include <vector>
+
+#include <fmt/format.h>
 
 namespace {
 
@@ -59,6 +65,24 @@ std::string EclString(Opm::Inplace::Phase phase) {
     }
 }
 
+    std::size_t numCells(const Opm::EclipseState& eclState)
+    {
+        return eclState.fieldProps().get_int("FIPNUM").size();
+    }
+
+    std::vector<Opm::EclInterRegFlowMap::SingleRegion>
+    defineInterRegionFlowArrays(const Opm::EclipseState&  eclState,
+                                const Opm::SummaryConfig& summaryConfig)
+    {
+        auto regions = std::vector<Opm::EclInterRegFlowMap::SingleRegion>{};
+
+        const auto& fprops = eclState.fieldProps();
+        for (const auto& arrayName : summaryConfig.fip_regions_interreg_flow()) {
+            regions.push_back({ arrayName, std::cref(fprops.get_int(arrayName)) });
+        }
+
+        return regions;
+    }
 }
 
 namespace Opm {
@@ -82,6 +106,7 @@ EclGenericOutputBlackoilModule(const EclipseState& eclState,
     , schedule_(schedule)
     , summaryConfig_(summaryConfig)
     , summaryState_(summaryState)
+    , interRegionFlows_(numCells(eclState), defineInterRegionFlowArrays(eclState, summaryConfig))
     , enableEnergy_(enableEnergy)
     , enableTemperature_(enableTemperature)
     , enableSolvent_(enableSolvent)

ebos/eclgenericoutputblackoilmodule.hh

@@ -43,6 +43,8 @@
 #endif
 #include <dune/common/parallel/mpihelper.hh>
 
+#include <ebos/eclinterregflows.hh>
+
 namespace Opm {
 
 namespace data { class Solution; }
@@ -370,6 +372,9 @@ protected:
     const Schedule& schedule_;
     const SummaryConfig& summaryConfig_;
     const SummaryState& summaryState_;
+
+    EclInterRegFlowMap interRegionFlows_;
+
     bool enableEnergy_;
     bool enableTemperature_;
 

ebos/eclinterregflows.cc

@@ -0,0 +1,212 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*
+  Copyright 2022 Equinor AS
+
+  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 <http://www.gnu.org/licenses/>.
+*/
+
+#include <ebos/eclinterregflows.hh>
+
+#include <algorithm>
+#include <stdexcept>
+
+Opm::EclInterRegFlowMapSingleFIP::
+EclInterRegFlowMapSingleFIP(const std::vector<int>& region)
+    : region_(region.size(), 0)
+{
+    if (! region.empty()) {
+        this->maxLocalRegionID_= this->maxGlobalRegionID_ =
+            *std::max_element(region.begin(), region.end());
+    }
+
+    std::transform(region.begin(), region.end(),
+                   this->region_.begin(),
+                   [](const int regID) { return regID - 1; });
+}
+
+void
+Opm::EclInterRegFlowMapSingleFIP::
+addConnection(const Cell& source,
+              const Cell& destination,
+              const data::InterRegFlowMap::FlowRates& rates)
+{
+    if (this->isReadFromStream_) {
+        throw std::logic_error {
+            "Cannot add new connection to deserialised object"
+        };
+    }
+
+    if (! source.isInterior ||
+        (source.cartesianIndex > destination.cartesianIndex))
+    {
+        // Connection handled in different call.  Don't double-count
+        // contributions.
+        return;
+    }
+
+    const auto r1 = this->region_[ source.activeIndex ];
+    const auto r2 = this->region_[ destination.activeIndex ];
+
+    if (r1 == r2) {
+        // Connection is internal to a region.  Nothing to do.
+        return;
+    }
+
+    // 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()
+{
+    this->iregFlow_.compress(this->maxGlobalRegionID_);
+}
+
+void Opm::EclInterRegFlowMapSingleFIP::clear()
+{
+    this->iregFlow_.clear();
+    this->isReadFromStream_ = false;
+}
+
+const Opm::data::InterRegFlowMap&
+Opm::EclInterRegFlowMapSingleFIP::getInterRegFlows() const
+{
+    return this->iregFlow_;
+}
+
+std::size_t Opm::EclInterRegFlowMapSingleFIP::getLocalMaxRegionID() const
+{
+    return this->maxLocalRegionID_;
+}
+
+bool
+Opm::EclInterRegFlowMapSingleFIP::
+assignGlobalMaxRegionID(const std::size_t regID)
+{
+    if (regID < this->maxLocalRegionID_) {
+        return false;
+    }
+
+    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

@@ -0,0 +1,362 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*
+  Copyright 2022 Equinor AS
+
+  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 <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef ECL_INTERREG_FLOWS_MODULE_HH
+#define ECL_INTERREG_FLOWS_MODULE_HH
+
+#include <opm/output/data/InterRegFlowMap.hpp>
+
+#include <algorithm>
+#include <cstddef>
+#include <string>
+#include <utility>
+#include <vector>
+
+/// \file
+///
+/// MPI-aware facility for converting collection of tuples of region ID
+/// pairs and associate flow rates into a sparse (CSR) adjacency matrix
+/// representation of a graph.  Supports O(nnz) compression.
+
+namespace Opm {
+
+    /// Form CSR adjacency matrix representation of inter-region flow rate
+    /// graph provided as a list of connections between regions on local MPI
+    /// rank.  Pertains to a single FIP definition array (e.g., FIPNUM).
+    class EclInterRegFlowMapSingleFIP
+    {
+    public:
+        /// Minimal characteristics of a cell from a simulation grid.
+        struct Cell {
+            /// 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};
+        };
+
+        /// Constructor
+        ///
+        /// \param[in] region Local rank's FIP region definition array.
+        explicit EclInterRegFlowMapSingleFIP(const std::vector<int>& region);
+
+        /// Add flow rate connection between regions.
+        ///
+        /// \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();
+
+        /// Get read-only access to the underlying CSR representation.
+        ///
+        /// Mostly intended for summary output purposes.
+        const data::InterRegFlowMap& getInterRegFlows() const;
+
+        /// Retrieve maximum FIP region ID on local MPI rank.
+        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::size_t regID);
+
+        /// 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
+        {
+            // Note: this->region_ is *intentionally* omitted here.
+            buffer.write(this->maxGlobalRegionID_);
+            this->iregFlow_.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)
+        {
+            auto otherMaxRegionID = 0 * this->maxGlobalRegionID_;
+            buffer.read(otherMaxRegionID);
+            this->maxGlobalRegionID_ = std::max(this->maxGlobalRegionID_, otherMaxRegionID);
+
+            this->iregFlow_.read(buffer);
+            this->isReadFromStream_ = true;
+        }
+
+    private:
+        /// Zero-based FIP region IDs on local MPI rank.
+        std::vector<int> region_{};
+
+        /// Maximum one-based FIP region ID on local MPI rank.
+        std::size_t maxLocalRegionID_{0};
+
+        /// Maximum one-based FIP region ID for this FIP region definition
+        /// array across all MPI ranks.
+        std::size_t maxGlobalRegionID_{0};
+
+        /// Rank-local inter-regional flow map.
+        data::InterRegFlowMap iregFlow_{};
+
+        /// Whether or not this object contains contributions deserialised
+        /// 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

ebos/ecloutputblackoilmodule.hh

@@ -31,6 +31,8 @@
 #include <numeric>
 #include <optional>
 #include <stdexcept>
+#include <type_traits>
+#include <utility>
 
 #include <opm/models/blackoil/blackoilproperties.hh>
 
@@ -576,6 +578,94 @@ public:
         }
     }
 
+    /*!
+     * \brief Capture connection fluxes, particularly to account for inter-region flows.
+     *
+     * \tparam ActiveIndex Callable type, typically a lambda, that enables
+     *    retrieving the active index, on the local MPI rank, of a
+     *    particular cell/element.  Must support a function call operator of
+     *    the form
+     \code
+        int operator()(const Element& elem) const
+     \endcode
+     *
+     * \tparam CartesianIndex Callable type, typically a lambda, that
+     *    enables retrieving the globally unique Cartesian index of a
+     *    particular cell/element given its active index on the local MPI
+     *    rank.  Must support a function call operator of the form
+     \code
+        int operator()(const int activeIndex) const
+     \endcode
+     *
+     * \param[in] elemCtx Primary lookup structure for per-cell/element
+     *    dynamic information.
+     *
+     * \param[in] activeIndex Mapping from cell/elements to linear indices
+     *    on local MPI rank.
+     *
+     * \param[in] cartesianIndex Mapping from active index on local MPI rank
+     *    to globally unique Cartesian cell/element index.
+     */
+    template <class ActiveIndex, class CartesianIndex>
+    void processFluxes(const ElementContext& elemCtx,
+                       ActiveIndex&&         activeIndex,
+                       CartesianIndex&&      cartesianIndex)
+    {
+        const auto identifyCell = [&activeIndex, &cartesianIndex](const Element& elem)
+            -> EclInterRegFlowMap::Cell
+        {
+            const auto cellIndex = activeIndex(elem);
+
+            return {
+                cellIndex,
+                cartesianIndex(cellIndex),
+                elem.partitionType() == Dune::InteriorEntity
+            };
+        };
+
+        const auto timeIdx = 0u;
+        const auto& stencil = elemCtx.stencil(timeIdx);
+        const auto numInteriorFaces = elemCtx.numInteriorFaces(timeIdx);
+
+        for (auto scvfIdx = 0*numInteriorFaces; scvfIdx < numInteriorFaces; ++scvfIdx) {
+            const auto& face = stencil.interiorFace(scvfIdx);
+            const auto left  = identifyCell(stencil.element(face.interiorIndex()));
+            const auto right = identifyCell(stencil.element(face.exteriorIndex()));
+
+            const auto rates = this->
+                getComponentSurfaceRates(elemCtx, face.area(), scvfIdx, timeIdx);
+
+            this->interRegionFlows_.addConnection(left, right, rates);
+        }
+    }
+
+    /*!
+     * \brief Prepare for capturing connection fluxes, particularly to
+     *    account for inter-region flows.
+     */
+    void initializeFluxData()
+    {
+        // Inter-region flow rates.  Note: ".clear()" prepares to accumulate
+        // contributions per bulk connection between FIP regions.
+        this->interRegionFlows_.clear();
+    }
+
+    /*!
+     * \brief Finalize capturing connection fluxes.
+     */
+    void finalizeFluxData()
+    {
+        this->interRegionFlows_.compress();
+    }
+
+    /*!
+     * \brief Get read-only access to collection of inter-region flows.
+     */
+    const EclInterRegFlowMap& getInterRegFlows() const
+    {
+        return this->interRegionFlows_;
+    }
+
     template <class FluidState>
     void assignToFluidState(FluidState& fs, unsigned elemIdx) const
     {
@@ -772,6 +862,106 @@ private:
         }
     }
 
+    /*!
+     * \brief Compute surface level component flow rates across a single
+     *   intersection.
+     *
+     * \param[in] elemCtx Primary lookup structure for per-cell/element
+     *    dynamic information.
+     *
+     * \param[in] scvfIdx Linear index of current interior bulk connection.
+     *
+     * \param[in] timeIdx Historical time-point at which to evaluate dynamic
+     *    quantities (e.g., reciprocal FVF or dissolved gas concentration).
+     *    Zero for the current time.
+     *
+     * \return Surface level component flow rates.
+     */
+    data::InterRegFlowMap::FlowRates
+    getComponentSurfaceRates(const ElementContext& elemCtx,
+                             const Scalar          faceArea,
+                             const std::size_t     scvfIdx,
+                             const std::size_t     timeIdx) const
+    {
+        using Component = data::InterRegFlowMap::Component;
+
+        auto rates = data::InterRegFlowMap::FlowRates{};
+
+        const auto& extQuant = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
+
+        const auto alpha = getValue(extQuant.extrusionFactor()) * faceArea;
+
+        if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
+            const auto& up = elemCtx
+                .intensiveQuantities(extQuant.upstreamIndex(oilPhaseIdx), timeIdx);
+
+            using FluidState = std::remove_cv_t<std::remove_reference_t<
+                decltype(up.fluidState())>>;
+
+            const auto pvtReg = up.pvtRegionIndex();
+
+            const auto bO = getValue(getInvB_<FluidSystem, FluidState, Scalar>
+                                     (up.fluidState(), oilPhaseIdx, pvtReg));
+
+            const auto qO = alpha * bO * getValue(extQuant.volumeFlux(oilPhaseIdx));
+
+            rates[Component::Oil] += qO;
+
+            if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
+                const auto Rs = getValue(
+                    BlackOil::getRs_<FluidSystem, FluidState, Scalar>
+                    (up.fluidState(), pvtReg));
+
+                rates[Component::Gas]    += qO * Rs;
+                rates[Component::Disgas] += qO * Rs;
+            }
+        }
+
+        if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
+            const auto& up = elemCtx
+                .intensiveQuantities(extQuant.upstreamIndex(gasPhaseIdx), timeIdx);
+
+            using FluidState = std::remove_cv_t<std::remove_reference_t<
+                decltype(up.fluidState())>>;
+
+            const auto pvtReg = up.pvtRegionIndex();
+
+            const auto bG = getValue(getInvB_<FluidSystem, FluidState, Scalar>
+                                     (up.fluidState(), gasPhaseIdx, pvtReg));
+
+            const auto qG = alpha * bG * getValue(extQuant.volumeFlux(gasPhaseIdx));
+
+            rates[Component::Gas] += qG;
+
+            if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
+                const auto Rv = getValue(
+                    BlackOil::getRv_<FluidSystem, FluidState, Scalar>
+                    (up.fluidState(), pvtReg));
+
+                rates[Component::Oil]    += qG * Rv;
+                rates[Component::Vapoil] += qG * Rv;
+            }
+        }
+
+        if (FluidSystem::phaseIsActive(waterPhaseIdx)) {
+            const auto& up = elemCtx
+                .intensiveQuantities(extQuant.upstreamIndex(waterPhaseIdx), timeIdx);
+
+            using FluidState = std::remove_cv_t<std::remove_reference_t<
+                decltype(up.fluidState())>>;
+
+            const auto pvtReg = up.pvtRegionIndex();
+
+            const auto bW = getValue(getInvB_<FluidSystem, FluidState, Scalar>
+                                     (up.fluidState(), waterPhaseIdx, pvtReg));
+
+            rates[Component::Water] +=
+                alpha * bW * getValue(extQuant.volumeFlux(waterPhaseIdx));
+        }
+
+        return rates;
+    }
+
     const Simulator& simulator_;
 };
 

ebos/eclwriter.hh

@@ -38,8 +38,11 @@
 
 #include <ebos/eclgenericwriter.hh>
 
-#include <string>
+#include <dune/grid/common/partitionset.hh>
+
+#include <functional>
 #include <limits>
+#include <string>
 
 namespace Opm::Properties {
 
@@ -181,11 +184,12 @@ public:
         const auto localGroupAndNetworkData = simulator_.problem().wellModel()
             .groupAndNetworkData(reportStepNum);
 
-
         const auto localAquiferData = simulator_.problem().aquiferModel().aquiferData();
         const auto localWellTestState = simulator_.problem().wellModel().wellTestState();
         this->prepareLocalCellData(isSubStep, reportStepNum);
 
+        this->captureLocalFluxData();
+
         if (this->collectToIORank_.isParallel())
             this->collectToIORank_.collect({},
                                            eclOutputModule_->getBlockData(),
@@ -435,6 +439,42 @@ private:
         OPM_END_PARALLEL_TRY_CATCH("EclWriter::prepareLocalCellData() failed: ", simulator_.vanguard().grid().comm())
     }
 
+    void captureLocalFluxData()
+    {
+        const auto& gridView = this->simulator_.vanguard().gridView();
+        const auto timeIdx = 0u;
+
+        auto elemCtx = ElementContext { this->simulator_ };
+
+        const auto elemMapper = ElementMapper { gridView, Dune::mcmgElementLayout() };
+        const auto activeIndex = [&elemMapper](const Element& e)
+        {
+            return elemMapper.index(e);
+        };
+
+        const auto cartesianIndex = [this](const int elemIndex)
+        {
+            return this->cartMapper_.cartesianIndex(elemIndex);
+        };
+
+        this->eclOutputModule_->initializeFluxData();
+
+        OPM_BEGIN_PARALLEL_TRY_CATCH();
+
+        for (const auto& elem : elements(gridView, Dune::Partitions::interiorBorder)) {
+            elemCtx.updateStencil(elem);
+            elemCtx.updateIntensiveQuantities(timeIdx);
+            elemCtx.updateExtensiveQuantities(timeIdx);
+
+            this->eclOutputModule_->processFluxes(elemCtx, activeIndex, cartesianIndex);
+        }
+
+        OPM_END_PARALLEL_TRY_CATCH("EclWriter::captureLocalFluxData() failed: ",
+                                   this->simulator_.vanguard().grid().comm())
+
+        this->eclOutputModule_->finalizeFluxData();
+    }
+
     Simulator& simulator_;
     std::unique_ptr<EclOutputBlackOilModule<TypeTag>> eclOutputModule_;
     Scalar restartTimeStepSize_;