Reservoir coupling: Populate satellite data for summary output

opm/simulators/flow/EclGenericWriter.hpp

@@ -30,6 +30,8 @@
 
 #include <opm/models/parallel/tasklets.hpp>
 
+#include <opm/output/data/Groups.hpp>
+
 #include <opm/simulators/flow/CollectDataOnIORank.hpp>
 #include <opm/simulators/flow/Transmissibility.hpp>
 #include <opm/simulators/timestepping/SimulatorReport.hpp>
@@ -149,7 +151,8 @@ class EclGenericWriter
                      const Inplace*                                       initialInPlace,
                      const InterRegFlowMap&                               interRegFlows,
                      SummaryState&                                        summaryState,
-                     UDQState&                                            udqState);
+                     UDQState&                                            udqState,
+                     const data::ReservoirCouplingGroupRates*             rcGroupRates = nullptr);
 
     CollectDataOnIORankType collectOnIORank_;
     const Grid& grid_;
@@ -180,7 +183,7 @@ class EclGenericWriter
     mutable std::vector<std::vector<NNCdata>> outputNncGlobalLocal_; // here GLOBAL refers to level 0 grid, local to any LGR (refined grid)
 
     // Amalgamated NNCs: nncs between different LGRs. For example, nested refinement or neighboring LGRs.
-    // The cells belong to different refined level grids. 
+    // The cells belong to different refined level grids.
     // nnc.cell1 (NNA1 in *.EGRID) Level/Local Cartesian index of cell1 (in its level grid: level1)
     // nnc.cell2 (NNA2 in *.EGRID) Level/Local Cartesian index of cell2 (in its level grid: level2, with level2 > level1).
     // Equivalent to TRANLL in *.INIT
@@ -209,16 +212,16 @@ class EclGenericWriter
     bool isCartesianNeighbour_(const std::array<int,3>& levelCartDims,
                                const std::size_t levelCartIdx1,
                                const std::size_t levelCartIdx2) const;
-    
+
     bool isDirectNeighbours_(const std::unordered_map<int,int>& levelCartesianToActive,
                              const std::array<int,3>& levelCartDims,
                              const std::size_t levelCartIdx1,
                              const std::size_t levelCartIdx2) const;
-    
+
     auto activeCell_(const std::unordered_map<int,int>& levelCartToLevelCompressed,
                      const std::size_t levelCartIdx) const;
-    
-    
+
+
 
     /// Returns true if the given Cartesian cell index belongs to a numerical aquifer.
     /// If no numerical aquifer exists, always returns false.
@@ -285,7 +288,7 @@ class EclGenericWriter
     /// Allocates NNCdata  vectors for
     /// - TRANNNC for level grids,
     /// - TRANGL for NNCs between level zero grid and an LGR (refined level grid)
-    /// - TRANLL for NNCs between two different refined level grids. 
+    /// - TRANLL for NNCs between two different refined level grids.
     /// Only for CpGrid. Other grid types do not support refinement yet.
     void allocateAllNncs_(int maxLevel) const;
 };

opm/simulators/flow/EclGenericWriter_impl.hpp

@@ -536,7 +536,7 @@ computeTrans_(const std::vector<std::unordered_map<int,int>>&  levelCartToLevelC
 
             // For level-zero grid, level Cartesian indices coincide with the grid Cartesian indices.
             if (isNumAquCell_(originCartIdxIn) || isNumAquCell_(originCartIdxOut)) {
-                // Check there are no refined aquifer cells. 
+                // Check there are no refined aquifer cells.
                 assert(level == 0);
                 // Connections involving numerical aquifers are always NNCs
                 // for the purpose of file output.  This holds even for
@@ -597,7 +597,7 @@ EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>::
 isDirectNeighbours_(const std::unordered_map<int,int>& levelCartesianToActive,
                     const std::array<int,3>& levelCartDims,
                     const std::size_t levelCartIdx1,
-                    const std::size_t levelCartIdx2) const 
+                    const std::size_t levelCartIdx2) const
 {
     return isCartesianNeighbour_(levelCartDims, levelCartIdx1, levelCartIdx2)
         || directVerticalNeighbors(levelCartDims, levelCartesianToActive, levelCartIdx1, levelCartIdx2);
@@ -607,7 +607,7 @@ template<class Grid, class EquilGrid, class GridView, class ElementMapper, class
 auto
 EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>::
 activeCell_(const std::unordered_map<int,int>& levelCartToLevelCompressed,
-            const std::size_t levelCartIdx) const 
+            const std::size_t levelCartIdx) const
 {
     auto pos = levelCartToLevelCompressed.find(levelCartIdx);
     return (pos == levelCartToLevelCompressed.end()) ? -1 : pos->second;
@@ -634,7 +634,7 @@ allocateAllNncs_(int maxLevel) const
         //          outputAmalgamatedNnc_[0][0] -> NNCs between level 1 and level 2
         //          outputAmalgamatedNnc_[0][1] -> NNCs between level 1 and level 3
         //          outputAmalgamatedNnc_[1][2] -> NNCs between level 2 and level 3
-        this->outputAmalgamatedNnc_.resize(maxLevel-1); 
+        this->outputAmalgamatedNnc_.resize(maxLevel-1);
         for (int i = 0; i < maxLevel-1; ++i) {
             this->outputAmalgamatedNnc_[i].resize(maxLevel-1-i);
         }
@@ -660,7 +660,7 @@ exportNncStructure_(const std::vector<std::unordered_map<int,int>>& levelCartToL
 
     // Cartesian index mapper for the serial I/O grid
     const auto& equilCartMapper = *equilCartMapper_;
-    
+
     const auto& level0CartDims = equilCartMapper.cartesianDimensions();
 
     int maxLevel = this->equilGrid_->maxLevel();
@@ -702,7 +702,7 @@ exportNncStructure_(const std::vector<std::unordered_map<int,int>>& levelCartToL
 
                 const auto& [smallerLevel, smallerLevelCartIdx] = smallerPair;
                 const auto& [largerLevel, largerLevelCartIdx] = largerPair;
-                
+
                 auto t = this->globalTrans().transmissibility(c1, c2);
 
                 // ECLIPSE ignores NNCs with zero transmissibility
@@ -765,13 +765,13 @@ exportNncStructure_(const std::vector<std::unordered_map<int,int>>& levelCartToL
                     // We need to check whether an NNC for this face was also
                     // specified via the NNC keyword in the deck.
                     auto t = this->globalTrans().transmissibility(c1, c2);
-                    
+
                     if (level == 0) {
                         auto candidate = std::lower_bound(nncData.begin(), nncData.end(),
                                                           NNCdata { originCartIdxIn, originCartIdxOut, 0.0 });
                         const auto transMlt = transMult.getRegionMultiplierNNC(originCartIdxIn, originCartIdxOut);
                         bool foundNncEditr = false;
-                    
+
                         while ((candidate != nncData.end()) &&
                                (candidate->cell1 == originCartIdxIn) &&
                                (candidate->cell2 == originCartIdxOut))
@@ -802,7 +802,7 @@ exportNncStructure_(const std::vector<std::unordered_map<int,int>>& levelCartToL
                             continue;
                         }
                     }
-                    
+
                     // ECLIPSE ignores NNCs with zero transmissibility
                     // (different threshold than for NNC with corresponding
                     // EDITNNC above).  In addition we do set small
@@ -1012,7 +1012,8 @@ evalSummary(const int                                            reportStepNum,
             const Inplace*                                       initialInPlace,
             const InterRegFlowMap&                               interRegFlows,
             SummaryState&                                        summaryState,
-            UDQState&                                            udqState)
+            UDQState&                                            udqState,
+            const data::ReservoirCouplingGroupRates*             rcGroupRates)
 {
     if (collectOnIORank_.isIORank()) {
         const auto& summary = eclIO_->summary();
@@ -1044,6 +1045,7 @@ evalSummary(const int                                            reportStepNum,
             /* block_values            = */ &blockData,
             /* aquifer_values          = */ &aquiferData,
             /* interreg_flows          = */ &interreg_flows,
+            /* rc_group_rates          = */ rcGroupRates,
             /* inplace                 = */ {
                 /* current = */ &inplace,
                 /* initial = */ initialInPlace

opm/simulators/flow/EclWriter.hpp

@@ -53,6 +53,11 @@
 #include <opm/simulators/utils/ParallelRestart.hpp>
 #include <opm/simulators/utils/ParallelSerialization.hpp>
 
+#include <opm/simulators/flow/rescoup/ReservoirCouplingEnabled.hpp>
+#ifdef RESERVOIR_COUPLING_ENABLED
+#include <opm/simulators/flow/rescoup/ReservoirCouplingMaster.hpp>
+#endif
+
 #include <boost/date_time/posix_time/posix_time.hpp>
 
 #include <limits>
@@ -128,7 +133,7 @@ class EclWriter : public EclGenericWriter<GetPropType<TypeTag, Properties::Grid>
 
     typedef Dune::MultipleCodimMultipleGeomTypeMapper< GridView > VertexMapper;
 
-    enum { enableEnergy = getPropValue<TypeTag, Properties::EnergyModuleType>() == EnergyModules::FullyImplicitThermal || 
+    enum { enableEnergy = getPropValue<TypeTag, Properties::EnergyModuleType>() == EnergyModules::FullyImplicitThermal ||
            getPropValue<TypeTag, Properties::EnergyModuleType>() == EnergyModules::SequentialImplicitThermal };
     enum { enableMech = getPropValue<TypeTag, Properties::EnableMech>() };
     enum { enableSolvent = getPropValue<TypeTag, Properties::EnableSolvent>() };
@@ -320,6 +325,10 @@ class EclWriter : public EclGenericWriter<GetPropType<TypeTag, Properties::Grid>
             miscSummaryData["MSUMNEWT"] = this->simulation_report_.success.total_newton_iterations;
         }
 
+        // For reservoir coupling master: collect slave production/injection
+        // rates to pass through to Summary::eval() via DynamicSimulatorState.
+        const auto rcGroupRates = this->collectReservoirCouplingGroupRates_();
+
         {
             OPM_TIMEBLOCK(evalSummary);
 
@@ -344,7 +353,8 @@ class EclWriter : public EclGenericWriter<GetPropType<TypeTag, Properties::Grid>
                               this->outputModule_->initialInplace(),
                               interRegFlows,
                               this->summaryState(),
-                              this->udqState());
+                              this->udqState(),
+                              rcGroupRates ? &(*rcGroupRates) : nullptr);
         }
     }
 
@@ -752,6 +762,29 @@ class EclWriter : public EclGenericWriter<GetPropType<TypeTag, Properties::Grid>
     const Schedule& schedule() const
     { return simulator_.vanguard().schedule(); }
 
+    /// Collect reservoir coupling master group rates for Summary::eval().
+    /// Returns nullopt for non-RC simulations or non-master processes.
+    std::optional<data::ReservoirCouplingGroupRates> collectReservoirCouplingGroupRates_()
+    {
+#ifdef RESERVOIR_COUPLING_ENABLED
+        // Guard: only BlackoilWellModel has reservoir coupling support.
+        // CompWellModel (compositional) does not, so we use if constexpr
+        // to avoid compilation errors when EclWriter is instantiated with
+        // a compositional TypeTag.
+        using WellModelType = std::remove_cvref_t<
+            decltype(simulator_.problem().wellModel())>;
+        if constexpr (requires(WellModelType& wm) { wm.isReservoirCouplingMaster(); }) {
+            auto& wellModel = simulator_.problem().wellModel();
+            if (!wellModel.isReservoirCouplingMaster()) {
+                return std::nullopt;
+            }
+            return wellModel.reservoirCouplingMaster()
+                .collectGroupRatesForSummary();
+        }
+#endif
+        return std::nullopt;
+    }
+
     void prepareLocalCellData(const bool isSubStep,
                               const int  reportStepNum)
     {

opm/simulators/flow/rescoup/ReservoirCouplingMaster.cpp

@@ -497,6 +497,15 @@ updateMasterGroupNameOrderMap(
     this->master_group_name_order_[slave_name] = master_group_name_map;
 }
 
+template <class Scalar>
+data::ReservoirCouplingGroupRates
+ReservoirCouplingMaster<Scalar>::
+collectGroupRatesForSummary() const
+{
+    assert(this->report_step_data_);
+    return this->report_step_data_->collectGroupRatesForSummary();
+}
+
 // ------------------
 // Private methods
 // ------------------

opm/simulators/flow/rescoup/ReservoirCouplingMaster.hpp

@@ -154,6 +154,11 @@ class ReservoirCouplingMaster {
          this->logger_.setDeferredLogger(deferred_logger);
     }
     void setFirstSubstepOfSyncTimestep(bool value);
+
+    /// @brief Collect production/injection rates for all master groups.
+    /// @details Returns a ReservoirCouplingGroupRates struct that can be
+    ///   passed through DynamicSimulatorState to Summary::eval().
+    data::ReservoirCouplingGroupRates collectGroupRatesForSummary() const;
     // These are currently only used for unit testing
     void setSlaveActivationDate(int index, double date) { this->slave_activation_dates_[index] = date; }
     void setSlaveNextReportTimeOffset(int index, double offset);

opm/simulators/flow/rescoup/ReservoirCouplingMasterReportStep.cpp

@@ -320,6 +320,47 @@ setReportStepIdx(int report_step_idx)
     this->report_step_idx_ = report_step_idx;
 }
 
+template <class Scalar>
+data::ReservoirCouplingGroupRates
+ReservoirCouplingMasterReportStep<Scalar>::
+collectGroupRatesForSummary() const
+{
+    using RcPhase = ReservoirCoupling::Phase;
+    using RateKind = ReservoirCoupling::RateKind;
+
+    data::ReservoirCouplingGroupRates result;
+
+    const auto numSlaves = this->numSlaves();
+    for (unsigned int i = 0; i < numSlaves; ++i) {
+        const auto& masterGroupNames = this->getMasterGroupNamesForSlave(i);
+        for (const auto& groupName : masterGroupNames) {
+            // Production rates (positive values, SI units)
+            data::ReservoirCouplingGroupRates::ProductionRates prod;
+            prod.oil = this->getMasterGroupRate_(groupName, RcPhase::Oil, RateKind::ProductionSurface);
+            prod.gas = this->getMasterGroupRate_(groupName, RcPhase::Gas, RateKind::ProductionSurface);
+            prod.water = this->getMasterGroupRate_(groupName, RcPhase::Water, RateKind::ProductionSurface);
+            prod.resv =
+                this->getMasterGroupRate_(groupName, RcPhase::Oil, RateKind::ProductionReservoir)
+                + this->getMasterGroupRate_(groupName, RcPhase::Gas, RateKind::ProductionReservoir)
+                + this->getMasterGroupRate_(groupName, RcPhase::Water, RateKind::ProductionReservoir);
+            result.production[groupName] = prod;
+
+            // Injection rates (positive values, SI units)
+            for (const auto phase : {RcPhase::Oil, RcPhase::Gas, RcPhase::Water}) {
+                const double surfRate = this->getMasterGroupRate_(
+                    groupName, phase, RateKind::InjectionSurface);
+                const double resRate = this->getMasterGroupRate_(
+                    groupName, phase, RateKind::InjectionReservoir);
+                if (surfRate != 0.0 || resRate != 0.0) {
+                    const auto opmPhase = convertToOpmPhase(phase);
+                    result.injection[groupName][opmPhase] = {surfRate, resRate};
+                }
+            }
+        }
+    }
+    return result;
+}
+
 
 // ------------------
 // Private methods

opm/simulators/flow/rescoup/ReservoirCouplingMasterReportStep.hpp

@@ -22,6 +22,7 @@
 
 #include <opm/simulators/flow/rescoup/ReservoirCoupling.hpp>
 #include <opm/simulators/flow/rescoup/ReservoirCouplingMpiTraits.hpp>
+#include <opm/output/data/Groups.hpp>
 #include <opm/input/eclipse/Schedule/Schedule.hpp>
 #include <opm/simulators/utils/ParallelCommunication.hpp>
 #include <opm/common/OpmLog/OpmLog.hpp>
@@ -221,6 +222,10 @@ class ReservoirCouplingMasterReportStep {
     /// being processed.
     void setReportStepIdx(int report_step_idx);
 
+    /// @brief Collect production/injection rates for all master groups.
+    /// @return ReservoirCouplingGroupRates struct with per-group rates.
+    data::ReservoirCouplingGroupRates collectGroupRatesForSummary() const;
+
     /// @brief Check if a specific slave process has been activated
     /// @param index Index of the slave process
     /// @return true if the slave is activated, false otherwise