Reservoir coupling: make a master group's injection target reflect what the slaves are actually producing

opm/simulators/flow/rescoup/ReservoirCouplingSlave.cpp

@@ -460,6 +460,24 @@ markSlaveGroupsInSchedule(Schedule& schedule, const int report_step_idx)
     this->report_step_data_->markSlaveGroupsInSchedule(schedule, report_step_idx);
 }
 
+template <class Scalar>
+void
+ReservoirCouplingSlave<Scalar>::
+setWellsSolvedThisSyncStep(bool value)
+{
+    assert(this->report_step_data_);
+    this->report_step_data_->setWellsSolvedThisSyncStep(value);
+}
+
+template <class Scalar>
+bool
+ReservoirCouplingSlave<Scalar>::
+wellsSolvedThisSyncStep() const
+{
+    assert(this->report_step_data_);
+    return this->report_step_data_->wellsSolvedThisSyncStep();
+}
+
 // ------------------
 // Private methods
 // ------------------

opm/simulators/flow/rescoup/ReservoirCouplingSlave.hpp

@@ -161,6 +161,10 @@ class ReservoirCouplingSlave {
     /// @details Delegates to ReservoirCouplingSlaveReportStep.
     void markSlaveGroupsInSchedule(Schedule& schedule, int report_step_idx);
 
+    /// @brief Record whether the initial well solve for this sync step has run
+    /// @details Delegates to ReservoirCouplingSlaveReportStep
+    void setWellsSolvedThisSyncStep(bool value);
+
     const std::string& slaveGroupIdxToGroupName(std::size_t group_idx) const {
         return this->slave_group_order_.at(group_idx);
     }
@@ -187,10 +191,15 @@ class ReservoirCouplingSlave {
     /// MPI_Comm_disconnect() to complete - it is a collective operation.
     void receiveTerminateAndDisconnect();
 
+    /// @brief True once the initial well solve for this sync step has run
+    /// @details Delegates to ReservoirCouplingSlaveReportStep
+    bool wellsSolvedThisSyncStep() const;
+
 private:
     void checkGrupSlavGroupNames_();
     std::pair<double, bool> getGrupSlavActivationDateAndCheckHistoryMatchingMode_() const;
     bool historyMatchingMode_() const { return this->history_matching_mode_; }
+
     std::size_t numMasterGroups_() const { return this->slave_to_master_group_map_.size(); }
     //! \brief Receive the master's go-ahead after reporting our OK status.
     //!

opm/simulators/flow/rescoup/ReservoirCouplingSlaveReportStep.hpp

@@ -267,6 +267,17 @@ class ReservoirCouplingSlaveReportStep {
     /// @param cmode Production control mode dictated by the master
     void setMasterProductionTarget(const std::string& gname, const Scalar target, const Group::ProductionCMode cmode);
 
+    /// @brief Record whether the initial well solve for this sync step has run
+    /// @param value False before the solve, true after it
+    void setWellsSolvedThisSyncStep(bool value) { wells_solved_this_sync_step_ = value; }
+
+    /// @brief True once the initial well solve for this sync step has run
+    /// @details Before it, a well that opens in this report step still carries the
+    ///   rates SingleWellState::update_producer_targets() derived from its WCONPROD
+    ///   target, which must not be reported to the master as achieved production.
+    ///   See RescoupSendSlaveGroupData::unsolvedNewWellProductionRates_().
+    bool wellsSolvedThisSyncStep() const { return wells_solved_this_sync_step_; }
+
 
 private:
     /// @brief Generic helper method for sending data to the master process via MPI
@@ -301,6 +312,10 @@ class ReservoirCouplingSlaveReportStep {
     // Used to control reservoir coupling synchronization of summary data sent from
     // the slave to the master process.
     bool is_last_substep_of_sync_timestep_{false};
+    // Flag to track whether the initial well solve of this sync timestep has run.
+    // Cleared before the pre-solve send of slave group data and set after the solve;
+    // see wellsSolvedThisSyncStep().
+    bool wells_solved_this_sync_step_{false};
 
     // Master-imposed targets and corresponding control modes, received from the master
     // process at the beginning of each sync timestep. Cleared and repopulated on every

opm/simulators/wells/BlackoilWellModelRescoup.hpp

@@ -140,7 +140,9 @@ class BlackoilWellModelRescoup {
     /// Blocking receive of the per-master-group production targets and
     /// injection limits computed by the master.  The received values are
     /// written into the slave's group state via the receiver helper.
-    /// Called from the slave's beginTimeStep first-substep handshake.
+    /// Called from the slave's beginTimeStep first-substep handshake, and
+    /// once per network iteration from maybeSendSlaveGroupFlowToMaster_(),
+    /// where the master sends injection targets only.
     void receiveGroupConstraintsFromMaster();
 
     /// \brief Receive master-computed network-leaf node pressures and
@@ -242,6 +244,22 @@ class BlackoilWellModelRescoup {
     ///   activated slaves and gates the master's own iteration.
     bool masterNetworkHasMasterGroupLeavesForSlave_(std::size_t slave_idx) const;
 
+    /// \brief Master-side: recompute the master's group state from the slave
+    ///   rates just received and send the resulting injection targets on.
+    ///
+    /// Called immediately after each non-final receiveSlaveGroupData() inside
+    /// the network iteration, so that a target derived from slave production
+    /// (GCONINJE REIN, SALE or VREP) keeps up with the production the slaves
+    /// report as the network iteration proceeds.  Its slave-side counterpart
+    /// is the receiveGroupConstraintsFromMaster() in
+    /// BlackoilWellModel::maybeSendSlaveGroupFlowToMaster_().
+    void refreshAndSendInjectionTargets_();
+
+    /// \brief Send injection targets to each activated slave, replacing the
+    ///   ones the slaves are currently holding.  Production constraints are
+    ///   not resent.  Only called by refreshAndSendInjectionTargets_().
+    void sendMasterGroupInjectionTargetsToSlaves_();
+
     /// \brief Slave-side: true iff this slave's own deck put the named group
     ///   into its own surface network as a fixed-pressure node.
     ///

opm/simulators/wells/BlackoilWellModelRescoup_impl.hpp

@@ -109,6 +109,7 @@ maybeExchangeNetworkOuterIterationWithSlaves(bool more_network_update)
         if (!is_final) {
             // receive slaves' updated network_surface_rates for the next outer iteration.
             this->receiveSlaveGroupData();
+            this->refreshAndSendInjectionTargets_();
         }
     }
 }
@@ -144,6 +145,7 @@ maybeExchangeNetworkSubIterationWithSlaves()
     }
     this->sendMasterGroupNodePressuresToSlaves(/*is_final=*/false);
     this->receiveSlaveGroupData();
+    this->refreshAndSendInjectionTargets_();
 }
 
 template<typename TypeTag>
@@ -413,6 +415,37 @@ masterNetworkHasMasterGroupLeavesForSlave_(std::size_t slave_idx) const
     return false;
 }
 
+template<typename TypeTag>
+void
+BlackoilWellModelRescoup<TypeTag>::
+refreshAndSendInjectionTargets_()
+{
+    // Called right after a receiveSlaveGroupData() inside the network iteration.
+    // Fold the rates just received into the master's group state -- that is what
+    // recomputes the reinjection and voidage rates that a GCONINJE REIN, SALE
+    // or VREP target is built from -- and ship the resulting targets to the
+    // slaves, replacing the ones they are currently holding.  Without this the
+    // targets stay at the values computed in beginTimeStep(), from slave
+    // production of the previous sync step.
+    const int report_step_idx = this->well_model_.simulator().episodeIndex();
+    this->well_model_.updateAndCommunicateGroupData(
+        report_step_idx, /*update_wellgrouptarget=*/false);
+    this->sendMasterGroupInjectionTargetsToSlaves_();
+}
+
+template<typename TypeTag>
+void
+BlackoilWellModelRescoup<TypeTag>::
+sendMasterGroupInjectionTargetsToSlaves_()
+{
+    OPM_TIMEFUNCTION();
+    RescoupConstraintsCalculator<Scalar, IndexTraits> constraints_calculator{
+        this->well_model_.guideRateHandler(),
+        this->groupStateHelper()
+    };
+    constraints_calculator.recalculateInjectionTargetsAndSendToSlaves();
+}
+
 template<typename TypeTag>
 bool
 BlackoilWellModelRescoup<TypeTag>::

opm/simulators/wells/BlackoilWellModel_impl.hpp

@@ -476,6 +476,10 @@ namespace Opm {
 #ifdef RESERVOIR_COUPLING_ENABLED
         if (this->isReservoirCouplingSlave()) {
             if (this->reservoirCouplingSlave().isFirstSubstepOfSyncTimestep()) {
+                // The wells have not been solved for this sync step yet, so a well that
+                // opens in this report step still carries its WCONPROD target as its rate.
+                // See RescoupSendSlaveGroupData::collectSlaveGroupSurfaceProductionRates_().
+                this->reservoirCouplingSlave().setWellsSolvedThisSyncStep(false);
                 this->rescoupHelper_.sendSlaveGroupDataToMaster();
                 this->rescoupHelper_.receiveGroupConstraintsFromMaster();
                 this->rescoupHelper_.receiveCoupledNetworkActiveStatus();
@@ -538,6 +542,8 @@ namespace Opm {
 
 #ifdef RESERVOIR_COUPLING_ENABLED
         if (slave_needs_well_solution) {  // isReservoirCouplingSlave()
+            // The initial well solve above has run, so the well states now hold solved rates.
+            this->reservoirCouplingSlave().setWellsSolvedThisSyncStep(true);
             // Need to update group data based on new well solution.
             this->updateAndCommunicateGroupData(reportStepIdx, /*update_wellgrouptarget*/ false);
             this->rescoupHelper_.sendSlaveGroupDataToMaster();
@@ -2362,6 +2368,11 @@ namespace Opm {
         if (!is_final) {
             this->updateAndCommunicateGroupData(reportStepIdx, /*update_wellgrouptarget=*/false);
             this->rescoupHelper_.sendSlaveGroupDataToMaster();
+            // The master turns the rates just sent into fresh injection targets
+            // for the groups it controls through a derived GCONINJE mode, and
+            // sends them straight back.  See
+            // BlackoilWellModelRescoup::refreshAndSendInjectionTargets_().
+            this->rescoupHelper_.receiveGroupConstraintsFromMaster();
             return /*more_network_update=*/true;
         }
         return /*more_network_update=*/false;

opm/simulators/wells/GroupStateHelper.cpp

@@ -1155,39 +1155,8 @@ GroupStateHelper<Scalar, IndexTraits>::sumWellPhaseRates(bool res_rates,
     }
 
     for (const std::string& well_name : group.wells()) {
-        const auto well_index = this->wellState().index(well_name);
-        if (!well_index.has_value())
-            continue;
-
-        if (!this->wellState().wellIsOwned(well_index.value(), well_name)) // Only sum once
-        {
-            continue;
-        }
-
-        const auto& well_ecl = this->schedule_.getWell(well_name, this->report_step_);
-        // only count producers or injectors
-        if ((well_ecl.isProducer() && is_injector) || (well_ecl.isInjector() && !is_injector))
-            continue;
-
-        const auto& ws = this->wellState().well(well_index.value());
-        if (ws.status == Opm::Well::Status::SHUT)
-            continue;
-
-        const Scalar factor = well_ecl.getEfficiencyFactor(network)
-            * this->wellState().well(well_index.value()).efficiency_scaling_factor;
-        if (res_rates) {
-            const auto& well_rates = ws.reservoir_rates;
-            if (is_injector)
-                rate += factor * well_rates[phase_pos];
-            else
-                rate -= factor * well_rates[phase_pos];
-        } else {
-            const auto& well_rates = ws.surface_rates;
-            if (is_injector)
-                rate += factor * well_rates[phase_pos];
-            else
-                rate -= factor * well_rates[phase_pos];
-        }
+        rate += this->wellRateContributionToGroup(
+            well_name, phase_pos, res_rates, is_injector, network);
     }
     return rate;
 }
@@ -1615,6 +1584,41 @@ GroupStateHelper<Scalar, IndexTraits>::updateWellRatesFromGroupTargetScale(
     }
 }
 
+template <typename Scalar, typename IndexTraits>
+Scalar
+GroupStateHelper<Scalar, IndexTraits>::wellRateContributionToGroup(const std::string& well_name,
+                                                                  const int phase_pos,
+                                                                  const bool res_rates,
+                                                                  const bool is_injector,
+                                                                  const bool network) const
+{
+    const auto well_index = this->wellState().index(well_name);
+    if (!well_index.has_value())
+        return 0.0;
+
+    if (!this->wellState().wellIsOwned(well_index.value(), well_name)) // Only sum once
+    {
+        return 0.0;
+    }
+
+    const auto& well_ecl = this->schedule_.getWell(well_name, this->report_step_);
+    // only count producers or injectors
+    if ((well_ecl.isProducer() && is_injector) || (well_ecl.isInjector() && !is_injector))
+        return 0.0;
+
+    const auto& ws = this->wellState().well(well_index.value());
+    if (ws.status == Opm::Well::Status::SHUT)
+        return 0.0;
+
+    const Scalar factor = well_ecl.getEfficiencyFactor(network)
+        * ws.efficiency_scaling_factor;
+    const auto& well_rates = res_rates ? ws.reservoir_rates : ws.surface_rates;
+    // Production rates are negative in the well state; a group rate sum is a positive
+    // magnitude for both directions.
+    return is_injector ? factor * well_rates[phase_pos]
+                       : -factor * well_rates[phase_pos];
+}
+
 template <typename Scalar, typename IndexTraits>
 std::pair<std::optional<std::string>, Scalar>
 GroupStateHelper<Scalar, IndexTraits>::worstOffendingWell(const Group& group,

opm/simulators/wells/GroupStateHelper.hpp

@@ -517,6 +517,28 @@ class GroupStateHelper
                                              bool is_injector,
                                              WellState<Scalar, IndexTraits>& well_state) const;
 
+    /// \brief The contribution of a single well to a group rate sum.
+    ///
+    /// Applies the filters, the efficiency factor and the sign convention that
+    /// sumWellPhaseRates() uses, so a caller that needs to add or remove one well's
+    /// share of such a sum stays consistent with it.  Returns zero when the well does
+    /// not contribute at all: not present in this rank's well state, not owned by this
+    /// rank (each well is counted once), of the wrong type for the sum, or shut.
+    ///
+    /// Being owner-filtered, the result is rank-local, exactly like sumWellPhaseRates();
+    /// a caller that compares it against an already-reduced quantity must reduce it too.
+    ///
+    /// \param well_name Name of the well
+    /// \param phase_pos Active phase index
+    /// \param res_rates Use reservoir rates instead of surface rates
+    /// \param is_injector Sum injectors rather than producers
+    /// \param network Use the network efficiency factors (GEFAC/WEFAC item 3)
+    Scalar wellRateContributionToGroup(const std::string& well_name,
+                                       const int phase_pos,
+                                       const bool res_rates,
+                                       const bool is_injector,
+                                       const bool network = false) const;
+
     /// Returns the name of the worst offending well and its fraction (i.e. violated_phase / preferred_phase)
     std::pair<std::optional<std::string>, Scalar>
     worstOffendingWell(const Group& group,